%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%DIGITAL WAVEGUIDE

%

%Rebecca Kleinberger

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

 

c = 34300; % in cm/s

nbSection = 44;

x = [1:44];

tube = zeros(nbSection);

areaFunction = 3;

lenghtTube = 17.5; % in cm

        %when the tube is straight it gives resonance at frequencies

        % 500 1500 and 2500 Hz -> formants of the voyel e

l = lenghtTube/nbSection;

 

Aschwa = [0.431328241979789;0.316379065522625;0.243148068695599;0.250409764613464;0.231989112094445;0.308139234156677;0.698765089367914;1.16984954923104;1.61474757297871;1.45632861413407;1.29433697071196;1.49963769625094;1.59469289792708;1.75127143809991;1.81865183800231;1.81924917180826;1.71986538750220;1.81619478701684;1.89678726284165;1.96000115184812;2.03579078843900;1.91932253113801;1.90175906840462;1.92038186975954;1.84707813678174;1.80887380983335;1.64765544920289;1.56985518037548;1.64985354985459;1.90641322036454;2.09957308321170;2.27285000521050;2.55968692599228;2.91515943155468;3.29500790040448;3.49166562945021;3.72830002798708;3.72336741105462;3.20774857637792;2.48502579498848;2.14496653113607;1.79301544179273;1.47889814281895;1.23556103985363];

 

Ai = [0.33 0.30 0.36 0.33 0.64 0.46 1.70 3.14 2.89 2.45 2.87 3.71 3.77 3.92 4.50 4.44 4.71 4.44 4.47 4.15 4.17 3.51 2.98 2.10 1.69 1.44 1.13 0.72 0.39 0.33 0.21 0.10 0.08 0.27 0.21 0.26 0.45 0.21 0.43 0.77 1.69 2.06 2.01 1.58];

       

Aa = [0.45 0.20 0.26 0.21 0.32 0.30 0.33 1.05 1.12 0.85 0.63 0.39 0.26 0.28 0.23 0.32 0.29 0.28 0.40 0.66 1.20 1.05 1.62 2.09 2.56 2.78 2.86 3.02 3.75 4.60 5.09 6.02 6.55 6.29 6.27 5.94 5.28 4.70 3.87 4.13 4.25 4.27 4.69 5.03];

 

Au = [0.40 0.36 0.29 0.44 0.69 2.15 3.00 2.72 2.15 2.48 4.95 5.91 5.49 5.05 4.60 4.41 3.77 3.39 3.18 3.29 3.24 3.33 2.08 2.04 1.42 0.62 0.18 0.17 0.22 0.25 0.46 0.71 0.75 1.33 2.23 2.45 3.16 5.16 4.92 2.73 1.21 0.79 0.42 0.86];

 

Ae = [0.23 0.13 0.14 0.19 0.04 0.26 1.08 1.26 1.21 0.96 0.72 0.74 0.91 1.64 1.91 2.20 2.62 2.77 2.98 3.00 2.83 2.84 2.86 2.44 2.15 2.10 1.84 1.77 1.84 1.72 1.45 1.37 1.36 1.43 1.72 2.08 2.36 2.66 2.38 1.95 2.68 2.61 2.19 1.60 ];

A= Ae;

 

f1 = figure('name','Section of the acoustic tubes');

set(f1, 'OuterPosition', [0 0 1100 650]);

bar(A);

xlabel('distance from the glottis');

ylabel('midsagittal distance to cross-sectional area');

 

%sampling frequency

Fs = 44100;

rho = 0.00114;

 

%number of semple we will extract to be able to reconstruct

 

nSamples = 4096*10;

%conditions at the limites??

 

% No of speech samples as output

nSamples = 4096*10;

 

 

%%

% The unit impulse/glottal volume velocity input (then use of rho in

% computation of glottis entry

u = zeros(nSamples, 1);

u(1)= 200;

 

%Creation of a Rosenberg pulse

F0 = 200 %Hz

v = [1:nSamples];

v = 200*sin(2*3.1415*218*v/nSamples);

 

glottal_tense = 0.85;

pitch = 120;

fs = Fs;

signal = [];

while (length(signal)<nSamples)

    signal = [signal rosenberg(glottal_tense, .65, pitch, fs)];

end

 

u =signal;

 

f2 = figure('name','Input signal from the glottis (either impulst signal of time varying, self sustained');

set(f2, 'OuterPosition', [0 0 1100 650]);

plot(u);

 

xlabel('sample number');

ylabel('glottis input');

 

 

 

%% delay arrays

%

p_plus = zeros(nbSection, 1);

p_minus = zeros(nbSection, 1);

Plips = zeros(nSamples, 1);

 

 

%% Coefficients of scattering

%

r = (A(1:(nbSection-1))-A(2:nbSection))./(A(1:nbSection-1)+A(2:nbSection));

r_glottis = 0.99; r_lips = -0.99;

 

 

 

%% COmputation of samples (scattering equations)

%

%

for n=1:nSamples

 

   

    %for the section 1: reflexion + glottis source

    p_plus(1) = r_glottis*p_minus(1) + u(n)*(rho*c)/A(1);

 

    %for all the other junctions

    % different between the even and the odd (for not infinite

    % recursion-ray tracing)

   

   

    % Scattering equations for even dD

    for i=2:2:nbSection-2

        theta  = r(i)*( p_plus(i)-p_minus(i+1) );

        p_plus(i+1) = p_plus(i) + theta;

        p_minus(i)   = p_minus(i+1) + theta;

    end

 

 

    % Scattering equations for odd dD

    for i=1:2:nbSection-1

        theta  = r(i)*( p_plus(i)-p_minus(i+1) );

        p_plus(i+1) = p_plus(i) + theta;

        p_minus(i)   = p_minus(i+1) + theta;

    end

 

 

    % last one = reflexion +lips

    p_minus(44) = r_lips*p_plus(44);

 

 

    % summ to get P at the lips

    Plips(n) = p_plus(44)+p_minus(44);

end

 

% Plips now contains the impulse response of the vocal tract filter

 

 

%% Frequency response and formant plot

 

f3 = figure('name','Frequency response and formant plot');

set(f3, 'OuterPosition', [0 0 1100 650]);

 

Y = fft(Plips, nSamples);

f = 0:Fs/nSamples:(nSamples-1)*Fs/nSamples;

Z = abs(Y);

yaxis = 20*log10(Z); % dB scale

ha = axes;

 

plot(ha, f,yaxis);

set(ha, 'Xlim',[0 5000],'Ylim', [-inf inf]);

grid on; xlabel('Frequency (Hz)'); ylabel('Frequency response of filter (dB)');

 

%%

 

wavwrite(Plips, fs, 16,'voyelle.wav');