%Project 2, Part A

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

I=[4 4 3.2 1.6 1.3 1 0.65 0.65 0.65 0.65 0.65 0.65 0.65 1.3 2.6 4 6.5 8 8 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 8 8 2 2 2.6 3.2];

O=[3.2 3.2 3.2 3.2 6.5 13 13 16 13 10.5 10.5 8 8 6.5 6.5 5 5 4 3.2 2 1.6 2.6 1.3 0.65 0.65 1 1 1.3 1.6 2 3.2 4 5 5 1.3 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];

Fs = 11025;

allVowels = {A,E,I,O,U};
vowelChars = ['a' 'e' 'i' 'o' 'u'];


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


for r=1:5
    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
    A=areaFunction;
    for k=1:length(areaFunction)-3
        rCoeff(k)=(A(k+1)-A(k))/(A(k+1)+A(k));
    end

    % plot vocal tract x-sectional area
    figure(3*(r-1)+1)
    subplot(211)
    stairs(tubePosition,areaFunction);
    title(['Cross-Sectional Area Function for /',vowelChars(r),'/ , N = ',num2str(N)]);
    xlabel('Position (m)'),ylabel('Area (cm^2)'),grid on,axis tight;


    % plot reflection coefficients
    subplot(212);
    stem(rCoeff)
    title('Relfection Coefficients');
    xlabel('k'),ylabel('Amplitude'),grid on;

    % First we'll find the lossless TF from the book
    lossless=[rCoeff 1];
    Dz=1;
    for k=2:length(lossless)
        Dz=[Dz 0];
        Dz=Dz+lossless(k)*fliplr(Dz);
    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


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

    [Vlossless,w]=freqz(num,Dz);
    [Vlossy,w]=freqz(num,Dz_lossy);

    % plot V(z)
    figure(3*(r-1)+2)
    clf;
    plot(w/pi,db(Vlossless),'r'); hold on;
    plot(w/pi,db(Vlossy),'b--'); axis tight, grid on;
    legend('Lossless',['Lossy \alpha = ' num2str(alpha)]);
    xlabel('Normalized Frequency (\omega/\pi)'), ylabel('20 log|V(\omega)|');
    title(['V(z) for /',vowelChars(r),'/']);

    % with R(z)
    [R,w]=freqz([1 -alpha]);
    H=Vlossy.*R;

    figure(3*(r-1)+3)
    plot(w/pi,db(H));
    axis tight; grid on;
    title(['TF with Radiation model for /',vowelChars(r),'/']);
    xlabel('Normalized Frequency (\omega/\pi)'), ylabel('20 log|V(\omega)|');


end

Published with MATLAB® 7.7