% M08-phasedelay.diary
% Design a filter with nearly flat frequency response, but lumpy phase
% response, to illustrate phase delay effects
% dpwe@ee.columbia.edu
% 2007-10-16

% But first: Visualization of shifting phase and group delay
% press a key once to see group delay increased without changing 
% phase delay, then press a key again to see phase delay changed 
% without changing group delay.
phasegroupdelay(100)

% We'll look at allpass filters later in the semester, but they 
% have flat frequency response with nonflat phase response.  They are 
% characterized by reciprocal poles and zeros
%
% Two poles, quite close together, quite close to the unit circle, for 
% large and abrupt phase change -> group delay
z0 = .9*exp(j*.1*pi);
z1 = .95*exp(j*.11*pi);
% Add another root on the real axis, just to make the frequency 
% response not completely flat
p = 0.99;
% Generate denominator polynomial from these roots + cpx conjugates
% (single-quote operator is conjugate *and* transpose)
a = poly([z0 z0' z1 z1' p*p]);
% Numerator of the allpass comes from reciprocal roots
b = poly(1./[z0 z0' z1 z1' p]);  % Real zero p just outside real pole p^2
% Check the pole/zero layout
zplane(b,a)
% zeros at reciprocal-pole locations
% Check the frequency response
freqz(b,a)
% very flat gain, far from flat phase response
% Compare phase response, phase delay, and group delay
[H,W] = freqz(b,a);
subplot(311)
plot(W/pi, unwrap(angle(H)));
title('Phase response (\theta(\omega))');
subplot(312)
plot(W/pi, unwrap(angle(H))./W);
title('Phase delay');
subplot(313)
grpdelay(b,a)
title('Group delay');
% peak group around w = 0.1pi is about 55 samples

% Apply to a soundfile.  First read the sound
[d,sr]=wavread('mpgr1_sx419.wav');
% Plot the spectrogram 
% (short-time fourier transform, shows energy in time and frequency)
subplot(211)
specgram(d,256,sr)
% Pass the sound through the allpass filter 100 times
e = d;
for i = 1:100; e = filter(b,a,e); end
% Spectogram post filtering
subplot(212)
specgram(e,256,sr)
% Clear delay around 0.1 pi
% What does it sound like?
soundsc(e,sr)
% weird