%Project 2, Part C - WOW

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] ;

U=[0.65 0.65 0.32 0.32 2 5 10.5 13 13 13 13 10.5 8 6.5 5 3.2 2.6 2 2 2 1.6 1.3 2 1.6 1 1 1 1.3 1.6 3.2 5 8 8 10.5 10.5 10.5 2 2 2.6 2.6];

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

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


Fs = 11025;

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

b=zeros(2,14);

allVowels = {A U 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 = [127 130 133 132 132];
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)
tU1 = .1;
tA = .35;
tU2 = .2;
tTransition = .05;

%then get appropriate glottal pulses for each
gpU1 = repmat(cell2mat(pulse(1)),1,round(Fs*tU1/length(cell2mat(pulse(1)))));
gpA = repmat(cell2mat(pulse(3)),1,round(Fs*tA/length(cell2mat(pulse(3)))));
gpU2 = repmat(cell2mat(pulse(5)),1,round(Fs*tU2/length(cell2mat(pulse(5)))));
gpTransition1 = repmat(cell2mat(pulse(2)),1,round(Fs*tTransition/length(cell2mat(pulse(2)))));
gpTransition2 = repmat(cell2mat(pulse(4)),1,round(Fs*tTransition/length(cell2mat(pulse(4)))));


U1 = filter([a -alpha],b(2,:),gpU1);
Aa = filter([a -alpha],b(1,:),gpA);
U2 = filter([a -alpha],b(2,:),gpU2);
Transition1 = filter([a -alpha],b(3,:),gpTransition1);
Transition2 = filter([a -alpha],b(3,:),gpTransition2);

U1 = U1./max(abs(U1));
Aa = Aa./max(abs(Aa));
U2 = U2./max(abs(U2));
Transition1 = Transition1./max(abs(Transition1));
Transition2 = Transition2./max(abs(Transition2));

s = [U1 Transition1 Aa Transition2 U2];

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

s = smooth(s);

soundsc(s,Fs)

Published with MATLAB® 7.7