%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 橢圓濾波器設計：無論信號頻率怎么改變，橢圓濾波器系數都不隨之改變，因此可以跟蹤
%                上信號頻率，做到頻率跟蹤濾波。而且Goertzel算法里的
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
N=3048;N1=2048;
f1=50;f2=3*f1;f3=5*f1;fs=f1*256;Rp=0.5;Rs=60;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%左路信號
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pn=0.000001;
%pn=0.0001;
u=randn(1,N);u=u*sqrt(pn);
real_left=8.00*pi/180; % hudu pi/180=0.0174532925
fi2=2*pi*rand(1);
fi3=2*pi*rand(1);
n=1:N;
t=2*pi*n/fs;
xl=sin(f1*t+real_left*ones(1,N))+0.1*sin(f2*t+fi2*ones(1,N))+0.01*sin(f3*t+fi3*ones(1,N))+u;%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%右路信號
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=randn(1,N);u=u*sqrt(pn);
real_right=12.00*pi/180;
fi2=2*pi*rand(1) ;
fi3=2*pi*rand(1) ;
n=1:N;
t=2*pi*n/fs;
xr=sin(f1*t+real_right*ones(1,N))+0.1*sin(f2*t+fi2*ones(1,N))+0.01*sin(f3*t+fi3*ones(1,N))+u;%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% filter design
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %[n1,wn]=ellipord([1.6/256,2.4/256],[1.2/256,3/256],Rp,Rs);
 [b,a]=ellip(4,0.1,60,3/256);
 xl0=filter(b,a,xl);
 xr0=filter(b,a,xr);
 Y1=abs(fft(xl));
 Y2=abs(fft(xl0));
 subplot(2,1,1);
 n=1:N;
 plot(n,abs(Y1),n,Y2,'r');% 觀察濾波效果，結果顯示頻率改變，濾波器系數不用改變，都有很好的濾波器效果。
 for n=1:N1
    xl0(n)=xl0(n+1000);
    xr0(n)=xr0(n+1000);
end