% E04_4.m             09nov15   05apr22
% allpass filters

clear all;
load allpass_MIT.mat;
 sound(speech,11025); pause(4);
return

%% simple allpass filters
r = .9;  phi = .3*pi;
 for r = .1:.2:.9
% for phi = .02*pi:.4*pi:.82*pi
  po = [r*exp(j*phi), r*exp(-j*phi)];
  ze = [(1/r)*exp(j*phi), (1/r)*exp(-j*phi)];
% a different example with more singularities:
%  po = [.8*r*exp(j*phi), .8*r*exp(-j*phi), r*exp(j*phi*.1), r*exp(-j*phi*.1)];
%  ze = [(1/(.8*r))*exp(j*phi), (1/(.8*r))*exp(-j*phi), (1/r)*exp(j*phi*.1), (1/r)*exp(-j*phi*.1)];
  figure(1); zplane(ze',po');
  a = poly(po);
  b = fliplr(a);
  figure(2); freqz(b,a);
  figure(3); grpdelay(b,a); drawnow;

  yc = filter(b,a,speech);
  for i=1:100;  yc=filter(b,a,yc);  end; % 3, 300
  sound(yc,11025); pause();

  pause(.3);
 end;

%%  Oppenheim's homework
clear all;
load allpass_MIT.mat;

N = 50;  % 1, 50
x = [1, zeros(1,N*1000)]; % poor estimate if x too short!
yc = filter(b,a,x);
for i=1:N-1;  yc=filter(b,a,yc);  end;
figure(1); plot(yc);  % impulse response
Yc = fft(yc);  % frequency response
Yc = Yc(1:round(length(Yc)/2)); % from dc to pi
figure(2);
subplot(311); plot(20*log10(abs(Yc))); axis('tight'); % magnitude (dB)
aYc = unwrap(angle(Yc))*180/pi;
subplot(312); plot(aYc); axis('tight');  % phase (deg)
subplot(313); plot(-(aYc(2:end)-aYc(1:end-1))); % group delay (samples)
figure(3); zplane(b,a); % for a single filter pass

sound(speech,11025);
yc = filter(b,a,speech);
for i=1:N-1;  yc=filter(b,a,yc);  end;
pause(4); sound(yc,11025);
 
% note: for N>3 a direct realization of H(z)^N yields huge errors



