dpbutterwordbilinear.m

IIR数字滤波器巴特沃思设计方法

%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);
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

%using bilinear method Design a butterworth digital lowpass filter
wp=100*2*pi; ws=300*2*pi;Fs=1000;Rp=0.25;Rs=50;
[N,Wn]=buttord(wp,ws,Rp,Rs,'s');%估计滤波器最小阶数
[z,p,k]=buttap(N);           %模拟滤波器函数引用
[Bap,Aap]=zp2tf(z,p,k);
[b,a]=lp2lp(Bap,Aap,Wn);
[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('magnitude')


% 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,512);
Pyy=fy.*conj(fy)/512;
f=1000*(0:256)/512;
subplot(2,2,4);
plot(f,Pyy(1:257));
xlabel('frequency/Hz')
grid;