%Project 2, Part C - HEY

clear all;
close all;

A=[5 5 5 5 6.5 8 8 8 8 8 8 8 8 6.5 5 4 3.2 1.6 2.6 2.6 2 1.6 1.3 1 0.65 0.65 0.65 1 1.6 2.6 4 1 1.3 1.6 2.6] ;

E=[8 8 5 5 4 2.6 2 2.6 2.6 3.2 4 4 4 5 5 6.5 8 6.5 8 10.5 10.5 10.5 10.5 10.5 8 8 6.5 6.5 6.5 6.5 1.3 1.6 2 2.6];


Tran = interp1(A,[E zeros(1,length(A)-length(E))]);

nanLocations = isnan(Tran);
Tran = Tran(~nanLocations);


Fs = 11025;

tau=1/(2*Fs);
c=345;

b=zeros(2,14);

allVowels = {A E Tran};


%%%%%%%%%%%%%%%%%%%%
%this loop will get us our denominators for filtering the glottal pusle
%%%%%%%%%%%%%%%%%%%%

for r=1:3
    currentVowel = cell2mat(allVowels(r));
    currentArea=fliplr(currentVowel);

    %create vector to reference position of each area segment
    currentLength=(0:length(currentArea)-1)*.005;

    % total number of tubes needed
    N=ceil(max(currentLength)/(tau*c))+1;

    % create position vector for only the needed number of tubes
    tubePosition = linspace(0,max(currentLength),N);

    areaFunction=spline(currentLength,currentArea,tubePosition);

    % "imaginary tube"
    areaFunction=[areaFunction max(areaFunction)+20 max(areaFunction)+20];
    tubePosition=[tubePosition 2*tubePosition(end)-tubePosition(end-1) tubePosition(end)+2*(tubePosition(end)-tubePosition(end-1))];


    % reflection ceofficients
    Af=areaFunction;
    for k=1:length(areaFunction)-3
        rCoeff(k)=(Af(k+1)-Af(k))/(Af(k+1)+Af(k));
    end

    % Then the lossy TF
    alpha=.75;
    r_lossy=[rCoeff alpha];
    Dz_lossy=1;
    for k=2:length(r_lossy)
        Dz_lossy=[Dz_lossy 0];
        Dz_lossy=Dz_lossy+r_lossy(k)*fliplr(Dz_lossy);
    end


    b(r,:)=[Dz_lossy zeros(1,length(b)-length(Dz_lossy))];

    a=zeros(1,(floor(N/2)));
    a=[a 1];

end


%%%%%%%%%%%%%%%%%%%%
%now we will define our glottal pulse
%%%%%%%%%%%%%%%%%%%%

freqs = [130 128 126];
freqs = freqs - 25;
pulse = {0,0,0,0,0};
for r = 1:length(freqs)

    %Total # of samples in one pulse
    N0=round(1/freqs(r)*Fs);

    %This gives the 4:1 ratio and 50% duty cycle
    N1=round(.4*N0);
    N2=round(1/4*N1);
    p = [.5*(1-cos(pi*(1:N1)/N1)),cos(pi*((N1:N1+N2)-N1)/(2*N2)),zeros(1,N0-(N1+N2)-1)];
    pulse(r) ={p + .05*randn(1,length(p))};

end

%define total times for each part of wow (U A U)
tNoise = .08;
tVowelA = .15;
tTran = .05;
tVowelE = .15;

%then get appropriate glottal pulse
gpA = repmat(cell2mat(pulse(1)),1,round(Fs*tVowelA/length(cell2mat(pulse(1)))));
gpE = repmat(cell2mat(pulse(3)),1,round(Fs*tVowelE/length(cell2mat(pulse(3)))));
gpTran = repmat(cell2mat(pulse(2)),1,round(Fs*tVowelE/length(cell2mat(pulse(2)))));
wgn = randn(1,round(Fs*tNoise));
wgn = wgn./max(abs(wgn));

source = [wgn gpA];


s = filter([a -alpha],b(1,:),source);

s = [s filter([a -alpha],b(3,:),gpTran) filter([a -alpha],b(2,:),gpE)];

s = s.*hamming(length(s))';

s = smooth(s);

soundsc(s,Fs)

Published with MATLAB® 7.7