%COMB FILTER APPLIED TO SUM OF TWO TRIANGLE WAVEFORMS
%Triangle#1 (red)  period 0.02->harmonics 50k Hertz, k odd
%Triangle#2 (blue) period 0.04->harmonics 25k Hertz, k odd
%Sample=10000/second. Duration=1 second.
clear;%GOAL:Reject triangle#1 (red) keep triangle#2 (blue)
T=[[-50:49]'*ones(1,50);[50:-1:-49]'*ones(1,50)]/50;T1=T(:);
T=[[-100:99]'*ones(1,25);[100:-1:-99]'*ones(1,25)]/100;T2=T(:);
FT1=abs(fft(T1))/10000;FT2=abs(fft(T2))/10000;
t=[0:9999]/10000; f=[1:500];%for plotting
t1=[1:1000];t10=(t1-1)/10000;%for plotting
figure %Two triangles and their spectra
subplot(211),plot(t10,T1(t1),'r',t10,T2(t1),'b'),axis([0 0.1 -2 2])
subplot(212),plot(f-1,FT1(f),'r',f-1,FT2(f),'b'),title('SPECTRA')
%COMB: Eliminate harmonics at 50,150,250,350 Hertz
%Compute frequency response H(f):
A=input('Type real part of comb filter poles:')%-5
Z=j*2*pi*[-350:100:350]';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 comb filter H(f) and h((t)
subplot(211),plot(f-1,FT1(f),'r',f-1,FT2(f),'b',f-1,abs(H(f)),'g')
subplot(212),plot(t,h),title('IMPULSE RESPONSE OF COMB FILTER')
%Convolve sum of two triangles T1+T2 with impulse response h.
%Add in impulse separately (T1+T2+)
%Cyclic convolution since periodic.
%conv(T1+T2,h)/100000 only works if longer (transients)
Z=(T1+T2)+real(ifft(fft(T1+T2).*fft(h')/10000));
figure %Plot sum of two triangles and filtered sum
subplot(211),plot(t10,T1(t1)+T2(t1)),title('SUM OF TWO TRIANGLES')
subplot(212),plot(t10,Z(t1)),title('FILTERED SUM OF TWO TRIANGLES')
