%RESONATOR FILTER APPLIED TO TRIANGLE+NOISE
%Triangle (red)  period 0.02->harmonics 50k Hertz, k odd
%Sample=10000/second. Duration=1 second.
clear;%GOAL: Keep triangle and eliminate white noise
T=[[-50:49]'*ones(1,50);[50:-1:-49]'*ones(1,50)]/50;T=T(:);
N=randn(10000,1);SNR=10*log10(sum(T.*T)/sum(N.*N))%-4.82
Y=T+N;FY=abs(fft(Y))/10000;
t=[0:9999]/10000; f=[1:500];%for plotting
t1=[1:1000];t10=(t1-1)/10000;%for plotting
figure %Triangle+noise and its spectrum
subplot(211),plot(t10,Y(t1)),title('TRIANGLE+NOISE')
subplot(212),plot(f-1,FY(f)),title('SPECTRUM OF TRIANGLE+NOISE')
%RESONATOR: Emphasize harmonics at 50,150,250 Hertz
%Compute frequency response H(f):
A=input('Type real part of resonator filter poles:')%-5
Z=j*2*pi*[-250:100:250]';P=Z+A;
NUM=poly(Z);DEN=poly(P);
H=polyval(NUM,j*2*pi*f)./polyval(DEN,j*2*pi*f);H=1-H;
%Compute impulse response h(t):
[R P K]=residue(NUM,DEN);
h=-real(R.'*exp(P*t));
figure %Plot resonator filter H(f) and h((t)
subplot(211),plot(f-1,FY(f),'r',f-1,abs(H(f)),'g')
subplot(212),plot(t,h),title('IMPULSE RESPONSE OF RESONATOR FILTER')
%Convolve triangle+noise Y with impulse response h.
%Cyclic convolution since periodic.
%conv(T1+T2,h)/100000 only works if longer (transients)
Z=real(ifft(fft(Y).*fft(h')/10000));
figure %Plot triangle+noise and filtered triangle+noise
subplot(211),plot(t10,Y(t1)),title('TRIANGLE+NOISE')
subplot(212),plot(t10,Z(t1)),title('FILTERED TRIANGLE+NOISE')
