dpcheby1bilinear.m

来自「ｉｉｒ数字滤波器契比雪夫设计方法」· M 代码 · 共 44 行

%input signal
t=0:0.001:0.1;
x=sin(2*pi*100*t)+sin(2*pi*400*t);

xlabel('time/s')
fx=fft(x,512);
Pxx=fx.*conj(fx)/512;    %做功率谱分析
f=1000*(0:256)/512;      %坐标转换
xlabel('frequency/Hz')
grid;
save  x.dat -ASCII

%Impulse invariance method Design a cheby1 digital lowpass filter
wp=100*2*pi; ws=300*2*pi;Fs=1000;Rp=0.25;Rs=50;
[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');%估计滤波器最小阶数
[z,p,k]=cheb1ap(N,Rp);           %模拟滤波器函数引用
[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);
subplot(2,2,1);
plot(W*Fs/(2*pi),abs(H));
grid;
subplot(2,2,2);
plot(W*Fs/(2*pi),20*log10(abs(H)));
grid;
xlabel('frequency/Hz')
ylabel('Gain dB')

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