dpellipbilinear.m

ｉｉｒ数字滤波器椭圆设计方法

%input signal
t=0:0.001:0.1;
%x=sin(2*pi*100*t)+sin(2*pi*400*t);
x=sin(2*pi*100*t)+sin(2*pi*200*t)+sin(2*pi*300*t);
subplot(3,2,1);
plot(t,x);
xlabel('time/s')
fx=fft(x,512);
Pxx=fx.*conj(fx)/512;    %做功率谱分析
f=1000*(0:256)/512;      %坐标转换
subplot(3,2,2);
plot(f,Pxx(1:257));
xlabel('frequency/Hz')
grid;
save  x.dat -ASCII

%using bilinear method Design a ellip digital lowpass filter
wp=100*2*pi; ws=300*2*pi;Fs=1000;Rp=0.25;Rs=50;
[N,Wn]=ellipord(wp,ws,Rp,Rs,'s');%估计滤波器最小阶数
[z,p,k]=ellipap(N,Rp,Rs);              %模拟滤波器函数引用
[A,B,C,D]=zp2ss(z,p,k);
[AT,BT,CT,DT]=lp2lp(A,B,C,D,Wn); %频率转换
[b,a]=ss2tf(AT,BT,CT,DT);
[bz,az]=bilinear(b,a,Fs);
[H,W]=freqz(bz,az);
omega1=linspace(0,ws,500);
omega2=linspace(wp,ws,200);
omega3=linspace(ws,2*pi*500,500);
h1=-20*log10(abs(freqs(bz,az,omega1)));
h2=-20*log10(abs(freqs(bz,az,omega2)));
h3=-20*log10(abs(freqs(bz,az,omega3)));
fprintf('Rp=%.4f\n',max(h1));
fprintf('Rs=%.4f\n',min(h3));
subplot(3,2,3);
axis([0 600 -150 1]);
plot(W*Fs/(2*pi),20*log10(abs(H)));
%plot([omega1 omega2 omega3]/(2*pi),[h1 h2 h3]);
grid;
xlabel('frequency/Hz')
ylabel('Gain in dB')

% using filter to do with input signal 
y=filter(bz,az,x);
subplot(3,2,5);
plot(t,y);
xlabel('time/s')
fy=fft(y,512);
Pyy=fy.*conj(fy)/512;
f=1000*(0:256)/512;
subplot(3,2,6);
plot(f,Pyy(1:257));
xlabel('frequency/Hz')
grid;