% 2008-01-31 E6820 acoustics lecture matlab demos
% dpwe@ee.columbia.edu

% First, simple traveling wave simulation
help travel1.m
%  y = travel1(l,d,r,x,v)  Visualize reflected travelling wave in waveguide
%     l is the length of each waveguide in samples.  
%     d is the number of output samples to produce (default 10000).
%     r is the refelection FIR at the far end (default [-1])
%     x is optional initial input waveform ([1] by default)
%     v if present and zero suppresses plotting.
%  2001-02-01 dpwe@ee.columbia.edu.  Based on pluck1a.m

y = travel1(20,200,[-1],hanning(5)');  % 5 pt pulse, fixed end, 20 pt len

% Maybe have some 'dispersion' in reflection
y = travel1(20,200,[-.25 -.5 -.25],[1]);

% Waves travelling in both directions...
y = travel1(20,200,[-1],[1 0 0 0 0 0 0 0 0 0 1]);


% Animations of sinusoidal standing waves
help sintwavemov
%  m = sintwavemov(k,w,z,l,n)  Build a movie of sin travelling wave interference
%     For a pair of opposite-moving sinusoidal travelling waves of 
%     temporal frequency w and spatial frequency k where the ratio 
%     of forward to backward waves at *l* is z (complex), make n 
%     plots of l points of the wave, both forward, backward, sum and 
%     difference.  Record each plot in a movie, returned in m.
%  2001-02-02 dpwe@ee.columbia.edu e6820 acous

% 20 steps for one cycle physically, 40 time steps per cycle, open end
subplot(211)
m = sintwavemov(2*pi/20,2*pi/40,-1,33,40);
movie(m,10)
% You can see waves traveling each way 
%  - the sum (blue) is the pressure (zero at end)
%  - the diff (red) is the volume velocity u(x,t)


% Vowel simulation by simple pole pairs
% e.g. ah = 500 1500 2500
wc = pi/4000;  % so f*wc is normalized angular frq of pole
r = 0.97;
r2 = r*r;
aah = conv(conv([1 -2*r*cos(500*wc) r2],[1 -2*r*cos(1500*wc) r2]), [1 -2*r*cos(2500*wc) r2]);
freqz(1,aah,512,8000)
pt = (rem(1:8000,80)==0);
vah = filter(1,aah,pt);
soundsc(vah)

% 'ee' has formants at 250 Hz, 2300 Hz (say)
aee = conv([1 -2*r*cos(250*wc) r2],[1 -2*r*cos(2300*wc) r2]);
freqz(1,aee,512,8000)
vee = filter(1,aee,pt);
soundsc(vee)
soundsc(vah)

% oo around 330, 900
aoo = conv([1 -2*r*cos(wc*330) r2],[1 -2*r*cos(wc*900) r2]);
freqz(1,aoo,512,8000)
voo = filter(1,aoo,pt);
soundsc(voo)

soundsc([vah,voo/8,vee])
% Sounds OK... oo is louder because formants both low-frequency


% Plucked string simulation
help pluck1a
%  y = pluck1a(l,d,r,x,v) Plucked string synthesis via pair of waveguide 'rails'
%     l is the length of each waveguide in samples.  
%     d is the number of output samples to produce (default 10000).
%     r is the refelection FIR at the far end (default [-.25 -.5 -.25] if empty)
%     x is optional initial shape to plucked string (triangular by default)
%       if x is a scalar value, it is taken as the pluck point as a
%       proportion of the string length (default 0.5).
%     v if present and nonzero plots waves & string at each sample time.
%  2001-02-01 dpwe@ee.columbia.edu.  Based on Julius Smith's 1992 CMJ paper.

x = pluck1a(40,10000);
soundsc(x)
subplot(211)
specgram(x,512,8000)
% Odd harmonics only - plucked in the middle

% Frequency-dependent decay reflects gain of reflection...
freqz([.25 .5 .25],1,512,8000)

% Pure loss still decays, but no coloration
x0 = pluck1a(40,10000,[0 -.97 0]);
soundsc(x0)
subplot(211)
specgram(x0,512,8000)

plot(x0(1:300))
subplot(212)
plot(x(1:300))
% output is actually 'loss' at bridge - high-pass for disperse filt

% Pluck point affects initial excitation, hence harmonics
x2 = pluck1a(40,10000, [], 0.2);
subplot(211)
specgram(x,512,8000)
subplot(212)
specgram(x2,512,8000)

soundsc(x2)
soundsc(x)

% Visualization of string behavior
x1 = pluck1a(40,10000, [], 0.5,1);

% Notice how resultant wave (red) is fixed at ends, but 'unfolds' each way