matlab的数字滤波器的仿真怎么实现

fp = 30;

fs = 35;

Fs = 100;

wp = 2*pi*fp/Fs;

ws = 2*pi*fs/Fs;

wp = tan(wp/2);

ws = tan(ws/2); % 通带最大衰减为0.5dB，阻带最小衰减为40dB

[N, wn] = buttord(wp, ws, 0.5, 40, 's'); % 模拟低通滤波器极零点

[z, p, k] = buttap(N); % 由极零点获得转移函数参数

[b, a] = zp2tf(z, p, k); % 由原型滤波器获得实际低通滤波器

[B, A] = lp2lp(b, a, wp);

[bz, az] = bilinear(B, A, .5);

[h, w] = freqz(bz, az, 256, Fs);

figure

plot(w, abs(h))

grid on

图1 巴特沃斯数字低通滤波器

1-2基于Butterworth模拟滤波器原型，使用双线性状换设计数字滤波器：各参数值为：通带截止频率Omega=0.2*pi,阻带截止频率Omega=0.3*pi,通带波动值Rp=1dB,阻带波动值Rs=15dB,设Fs=4000Hz。

代码：

wp=0.2*pi;ws=0.3*pi;

Fs=4000;T=1/Fs;

OmegaP=(2/T)*tan(wp/2);

OmegaS=(2/T)*tan(ws/2);

rp=1;rs=15;as=15;

ripple=10^(-rp/20);attn=10^(-rs/20);

[n,wn]=buttord(OmegaP,OmegaS,rp,rs,'s');

[z,p,k]=Buttap(n);

[b,a]=zp2tf(z,p,k);

[bt,at]=lp2lp(b,a,wn);

[b,a]=bilinear(bt,at,Fs);

[db,mag,pha,grd,w]=freqz_m(b,a);

%

%下面绘出各条曲线

subplot(2,2,1);plot(w/pi,mag);title('Magnitude Frequency幅频特性');

xlabel('w(/pi)');ylabel('|H(jw)|');

axis([0,1,0,1.1]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);

set(gca,'YTickMode','manual','YTick',[0 attn ripple 1]);grid

subplot(2,2,2);plot(w/pi,db);title('Magnitude Frequency幅频特性(db)');

xlabel('w(/pi)');ylabel('dB');

axis([0,1,-30,5]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);

set(gca,'YTickMode','manual','YTick',[-60 -as -rp 0]);grid

subplot(2,2,3);plot(w/pi,pha/pi);title('Phase Frequency相频特性');

xlabel('w(/pi)');ylabel('pha(/pi)');

axis([0,1,-1,1]);

subplot(2,2,4);plot(w/pi,grd);title('Group Delay群延时');

xlabel('w(/pi)');ylabel('Sample');

axis([0,1,0,15]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);grid

运行结果：

图2巴特沃思数字低通滤波器幅频-相频特性

1-3设计一巴特沃斯高通数字滤波器，要求通带截止频率0.6*pi，通带衰减不大于1dB，阻带衰减15DB，抽样T=1。

Wp=0.6*pi;

Ws=0.4*pi;

Ap=1;

As=15;

[N,wn]=buttord(Wp/pi,Ws/pi,Ap,As); %计算巴特沃斯滤波器阶次和截止频率

%频率变换法设计巴特沃斯高通滤波器

[db,mag,pha,grd,w]=freqz_m(b,a); %数字滤波器响应

plot(w,mag);

title('数字滤波器幅频响应|H(ej\Omega)|')

图3巴特沃斯数字高通滤波器

2-1用窗函数法设计一个线性相位FIR低通滤波器，并满足性能指标:通带边界频率

Wp=0.5*pi,阻带边界频率Ws=0.66*pi,阻带衰减不小于40dB,通带波纹不大于3dB。选择汉宁窗。

代码：

wp =0.5*pi;

ws=0.66*pi;

wdelta =ws-wp;

N= ceil(8*pi/wdelta)

if rem(N,2)==0

N=N+1;

end

);