%NOTCH FILTER APPLIED TO SUM OF TRIANGLE AND SINE WAVEFORMS
%Sinusoid (red) frequency 50 Hertz.
%Triangle (blue) period 0.04->harmonics 25k Hertz, k odd
%Sample=10000/second. Duration=1 second.
clear;%GOAL:Reject sinusoid (red) keep triangle (blue)
t=[0:9999]/10000;S=cos(2*pi*50*t);S=S';%50 Hertz sinusoid
T=[[-100:99]'*ones(1,25);[100:-1:-99]'*ones(1,25)]/100;T=T(:);
FS=abs(fft(S))/10000;FT=abs(fft(T))/10000;
t1=[1:1000];t10=(t1-1)/10000;f=[1:200];%for plotting
figure %Sinusoid and triangle and their spectra
subplot(211),plot(t10,S(t1),'r',t10,T(t1),'b'),axis([0 0.1 -2 2])
subplot(212),plot(f-1,FS(f),'r',f-1,FT(f),'b'),title('SPECTRA')
%NOTCH: Eliminate 50 Hertz. Compute frequency response H(f):
A=input('Type real part of comb filter poles:')%-5
Z=j*2*pi*[-50 50]';P=Z+A;
NUM=poly(Z);DEN=poly(P);
H=polyval(NUM,j*2*pi*f)./polyval(DEN,j*2*pi*f);
%Compute impulse response h(t):
[R P K]=residue(NUM,DEN);
h=real(R.'*exp(P*t));
figure %Plot notch filter H(f) and h((t)
subplot(211),plot(f-1,FS(f),'r',f-1,FT(f),'b',f-1,abs(H(f)),'g')
subplot(212),plot(t,h),title('IMPULSE RESPONSE OF NOTCH FILTER')
%Convolve sum of sine and triangle S+T with impulse response h.
%Add in impulse separately (S+T+)
%Cyclic convolution since periodic.
%conv(S+T,h)/100000 only works if longer (transients)
Z=(S+T)+real(ifft(fft(S+T).*fft(h')/10000));
figure %Plot sum of two triangles and filtered sum
subplot(211),plot(t10,S(t1)+T(t1)),title('SINUSOID+TRIANGLE')
subplot(212),plot(t10,Z(t1)),title('FILTERED SINUSOID+TRIANGLE')
