%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%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');