% 2007-10-30-apfilt.diary
% Looking at allpass filters
% dpwe@ee.columbia.edu

% Define a simple pair of conjugate poles
z0 = .9*exp(j*.1*pi);
% poly() gives the polynomial coefficients for these roots
a = poly([z0 z0'])
%a =
%    1.0000   -1.7119    0.8100
% So, for an allpass filter, numerator coefficients are simply the reverse
b = fliplr(a);
% Look at the response
freqz(b,a)
% Magnitude response is flat to numerical resolution (10^-14)
% Where are poles/zeros?
zplane(b,a)
% zeros visibly at reciprocals (reflected in unit circle) of poles
% What is the group delay (matlab built-in does d theta/d omega for us)
grpdelay(b,a)
% Peak delay at point of peak rate of phase change


% Apply the filter repeatedly, to accentuate group delay
[x,sr]=wavread('mpgr1_sx419.wav');
y = x;
for i = 1:100; y = filter(b,a,y); end
soundsc(x,sr)
soundsc(y,sr)
% Weird sqeakiness from local group delay
% Check that long-term spectra unchanged
subplot(211)
plotspec(x,sr,16384);
% (plotspec.m is in the course matlab directory)
subplot(212)
plotspec(y,sr,16384);
% but spectrograms show clear effect
subplot(211)
specgram(x,512,sr)
subplot(212)
specgram(y,512,sr)