%Demo done during Session 6 illustrating
%the utility of a simple second-order difference
%equation to notch out an undesired tone while
%preserving the integrity of the desired signal
%as best as possible.  This demo also illustrates
%the effects of zeroes and poles on the frequency
%response of an LTI System described by a difference 
%eqn. FIR vs IIR notch filtered is explored as well.
%Although this code is not commented, the input lines
%serve as comments that guide one through the code.
clear
clf
set(0,'defaultaxesfontsize',16);

% read speech signal
data=load('erf1s1t0')';
Fs=8192;
omega0=pi/4;
NL=100;
plot(data)
title('Original speech signal')
%end1=input('end1?');
%end2=input('end2?');
%datar=data(end1:end2);
datar=data;

% generate corrupted signal
x=datar+500*cos(omega0*(1:length(datar)));
input('Utterance played back at orignal sampling rate, Fs = 8192 Hz');
soundsc(datar,Fs)
plot(x)
title('Corrupted speech signal')
input('Utterance with 1 KHz tone interference played at Fs = 8192 Hz');
soundsc(x,Fs)

% compute and plot filter frequency response of IIR filter
r=0.95; c0=2*cos(omega0); rs=r^2; c0r=r*c0;
h2poles=r.^(1:NL).*sin(omega0*(1:NL));
hiir=conv(h2poles,[1 -c0 1]);
poles=roots([1 -c0r rs]);
zeroes=roots([1 -c0 1]);
polar(1,0,'.');
hold on
plot(real(poles),imag(poles),'x','Linewidth',3,'MarkerSize',18)
plot(real(zeroes),imag(zeroes),'o','Linewidth',3,'MarkerSize',18)
title('Polo-zero plot of filter transfer function')
hold off

input('Pole-zero plot -- NEXT: magnitude of freq. response for notch filter')
domega=2*pi/256;
omega=-pi:domega:pi-domega;
hf1=abs(fftshift(fft(hiir,256)));
plot(omega,hf1,'LineWidth',4)
axis([-pi pi 0 max(hf1)]);
title('Magnitude of DTFT of IIR notch filter');
xlabel('omega (radians/sample)');

% filter speech signal with IIR filter
r=0.95; c0=2*cos(omega0); rs=r^2; c0r=r*c0;
disp('Run tone corrupted speech through simple second');
input(' order difference equation...');
y=zeros(1,length(x));
for n=3:length(x),
 y(n)=c0r*y(n-1)-rs*y(n-2)+x(n)-c0*x(n-1)+x(n-2);
end
%
%y=conv(x,hiir);
input('IIR Notch-filtered utterance plus 1.0 KHz tone interference');
soundsc(y,Fs)

% apply FIR filter
z=conv(x,[1 -c0 1]);
hf2=abs(fftshift(fft([1 -c0 1],256)));
plot(omega,hf2,'LineWidth',4)
axis([-pi pi 0 max(hf2)]);
title('Magnitude of DTFT of FIR notch filter');
xlabel('omega (radians/sample)');
input('FIR Notch-filtered utterance plus 1.0 KHz tone interference');
soundsc(z,Fs)
input('Pole-zero plot -- NEXT: magnitude of freq. response for notch filter')
domega=2*pi/8192;
omega=-pi:domega:pi-domega;
daf1=abs(fftshift(fft(datar,8192)));
plot(omega,daf1)
axis([-pi pi 0 max(daf1)]);
title('Magnitude of DTFT of original utterance');
xlabel('omega (radians/sample)');