经典滤波器的MATLAB仿真源程序

1、％巴特沃斯低通模拟圆形滤波器

clear all;

n=0:0.01:2;

for i=1:4

switch i

case 1

N=2;

case 2

N=5;

case 3

N=10;

case 4

N=20;

end

[z,p,k]=buttap(N); %函数buttap--设计巴特沃斯低通滤波器

[b,a]=zp2tf(z,p,k); %函数zp2tf--零极点增益模型转换为传递函数模型[H,w]=freqs(b,a, n); %函数freqs--求解模拟滤波器频率响应

magH2=(abs(H)).A2; %函数abs--取模值函数

hold on %函数hold--控制是否保持当前图形

plot(w,magH2) %函数plot--画二维线性图

axis([0 2 0 1]);

%函数axis--控制坐标轴比例和外观

end

xlabel('w/wc');

ylabel('|H(jw)$2');

title('巴特沃斯低通模拟滤波器');

text(0.72,0.63,'N=2') %对不同曲线做标记

text(0.98,0.85,'N=20')

grid on;

2、％绘制切比雪夫I型低通模拟滤波器的平方幅频响应曲线，滤波器的阶数分别为2,4,6,8.

clear all;

n=0:0.01:2;

for i=1:4

switch i

case 1

N=2;

case 2

N=4;

case 3

N=6;

case 4

N=8;

end

Rs=10;

[z,p,k]=cheb1ap(N,Rs);

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

[H,w]=freqs(b,a, n);

magH2=(abs(H)).A2;

posplot=['22' nu m2str(i)];

subplot(posplot)

plot(w,magH2)

axis([0 2 0 1]);

xlabel('w/wc');

ylabel('H(jw)A2');

title(['N=' nu m2str(N)]);

grid on end

3、％切比雪夫II型低通模拟滤波器

clear all;

n=0:0.01:2;

for i=1:2

switch i

case 1

N=7;

case 2

N=8;

end

Rs=10; %阻带文波系数为10dB

[z,p,k]=cheb2ap(N,Rs); %函数cheb2---设计切比雪夫II型低通滤波器

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

[H,w]=freqs(b,a, n);

magH2=(abs(H)).A2;

%输出图形

posplot=['12' num2str(i)];

subplot(posplot)

plot(w,magH2)

axis([0 2 0 1.1]);

xlabel('w/wc');

ylabel('|H(jw)|A2');

title(['N=' nu m2str(N)]);

end

4、％运用冲击响应不变法设计一个低通Chebshevl型数字滤波器，其通带上限临界频率是3Hz,阻带临界频率是5H，采样频率是1000Hz,在通带内的最大衰减为0.3dB，阻带内的最小衰减为80dB。

MATLAB程序如下：

clc;

clear all;

%把数字滤波器的频率特征转换成模拟滤波器的频率特征

wp=300*2*pi;

ws=400*2*pi;

rp=0.3;

rs=80;

Fs=1000;

%选择滤波器的最小阶数。

[N,Wn]=cheb1ord(wp,ws,rp,rs,'s');

%创建Chebyshev1低通滤波器的原型

[Z,P,K]=cheb1 ap(N,rp);

[A,B,C,D]=zp2ss(Z,P,K);

%实现低通向低通的转换

[AT,BT,CT,DT]=lp2lp(A,B,C,D,Wn);

[nu m1,de n1]=ss2tf(AT,BT,CT,DT);

%运用冲击响应不变法把模拟滤波器转换成数字滤波器[nu m2,de n2]=impi nvar( num1,de n1,1000);

%绘出频率响应曲线

[H,W]=freqz( nu m2,de n2);

plot(W*Fs/(2*pi),abs(H));

grid;

xlabel('幅值');

ylabel('频率');

title('冲击响应不变法低通滤波器');

clc;

clear all;

%把数字滤波器的频率特征转换成模拟滤波器的频率特征wp=300*2*pi;

ws=400*2*pi;

rp=0.3;

rs=80;

Fs=1000;

%选择滤波器的最小阶数。

[N,Wn]=cheb1ord(wp,ws,rp,rs,'s');

%创建Chebyshevl低通滤波器的原型

[Z,P,K]=cheb1 ap(N,rp);

[A,B,C,D]=zp2ss(Z,P,K);

%实现低通向低通的转换

[AT,BT,CT,DT]=lp2lp(A,B,C,D,Wn);

[nu m1,de n1]=ss2tf(AT,BT,CT,DT);

%运用冲击响应不变法把模拟滤波器转换成数字滤波器[nu m2,de n2]=impi nvar( num1,de n1,1000);

%绘出频率响应曲线

[H,W]=freqz( nu m2,de n2);

plot(W*Fs/(2*pi),abs(H));

grid;

xlabel('幅值');

ylabel('频率');

title('冲击响应不变法低通滤波器');

5、％使用双线性Z变换设计一低通数字滤波器，使数字滤波器满足以下参数：

采样频率Fs=1000HZ,通带临界频率fl =200Hz，通带内衰减小于1dB( a p=1);

阻带临界频率fh=300Hz，阻带内衰减大于25dB ( a s=25)。

FS=1000;

FI=200;Fh=300; %通带、阻带截止频率

Rp=1;Rs=25;

wp=FI*2*pi; %临界频率采用角频率表示

ws=Fh*2*pi; %临界频率采用角频率表示

wp1=wp/FS; %求数字频率

WS仁ws/FS; %求数字频率

OmegaP=2*FS*ta n(wp1/2);% 频率预畸

OmegaS=2*FS*ta n(ws1/2);%频率预畸

%选择滤波器的最小阶数

[n, Wn]=buttord(OmegaP,OmegaS,Rp,Rs,'s');%此处是代入经预畸变后获得的归一化模拟频率参数

[bt,at]=butter(n,Wn,'s'); %设计一个n阶的巴特沃思模拟滤波器

[bz, az]=bi lin ear(bt,at,FS); %双线性变换为数字滤波器

[H,W] = freqz(b z,az); %求解数字滤波器的频率响应pIot(W*FS/(2*pi),abs(H));grid; xlabel('频率/Hz');ylabel('幅值');

title('双线性Z变换设计低通数字滤波器')

6基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字高通滤波器，要求数字截止频率为-= 0.2二ra d ,通带内衰减不大于3dB,阻带起始频率为?‘S =0.05二rad,阻带内衰减不小于48dB。

wp=0.2*pi;ws=0.05*pi;Ap=3;As=48;

Wp=ta n(wp/2);Ws=ta n(ws/2);

[N,w n]=buttord(Wp,Ws,Ap,As,'s'); fprintf('滤波器阶数N=%.0f\n',N)

[b,a]=butter(N,1,'s');

[nu ma,de na]=lp2hp(b,a,Wp);

[nu md,de nd]=bili near( nu ma,de na,0.5); disp('分子系数b');

fprin tf('%.4e ',nu md);

fprin tf('\n');

disp('分母系数a');

fprin tf('%.4e ',de nd);

fprin tf('\n');

w=li nspace(0,pi,1024);

h=freqz( nu md,de nd,w);

plot(w/pi,20*log10(abs(h)));

axis([0 1 -50 0]);grid;

xlabel('归一化频率');

ylabel (幅度/dB');-

7、基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字带通滤波器，要求抽样频率Fs=2000HZ,通带范围为300~400Hz,在带边频率处衰减不大于3dB,

在200Hz以下和500Hz以上衰减不小于18dB。

Matlab程序如下：

clear all;

fp=[300 400];fs=[200 500];

rp=3; rs=18;

Fs=2000;

wp=fp*2*pi/Fs;

ws=fs*2*pi/Fs;

wap=2*Fs*ta n(wp./2)

was=2*Fs*ta n(ws./2);

[n,wn]=buttord(wap,was,rp,rs,'s');

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

[bp,ap]=zp2tf(z,p,k)

bw=wap(2)-wap(1) w0=sqrt(wap(1)*w ap(2));

[bs,as]=lp2bp(bp,ap,w0,bw) [h1,w1]=freqs(bp,ap);

figure(1)

plot(w1,abs(h1));grid;

ylabel('Ba ndpass AF and DF') xlabel('Hz')

8按下列要求用海明窗设计一数字低通滤波器：电力系统低频振荡频率在0.2?

2.5 Hz之间，故滤波器通带截至频率为3 Hz,阻带截止频率为5 Hz,通带内最大衰减不高于0.5 dB,阻带最小衰减不小于50 dB。采样频率为100Hz

fs=100; %采样频率

wp = 3*pi/50; ws = 5*pi/50; deltaw= ws - wp; % 过渡带宽△ w的计算N0 = ceil(6.6*pi/ deltaw) + 1; %按海明窗计算所需的滤波器阶数N0

N=N0+mod(N0+1,2); %为了实现第一类偶对称滤波器，应使其长度N为奇数w_ham = (hammi ng(N))'; %求窗函数

wc = (ws+wp)/2 ; %截止频率取为两边界频率的平均值

tao = (N-1)/2; %理想脉冲响应的对称中心位置

n = [0: (N-1)]; %设定脉冲响应长度

m = n - tao + eps; %加一个小数以避免零作除数

hd = si n( wc*m) ./ (pi*m); %理想脉冲响应

h = hd .* w_ham; %设计的脉冲响应(即系数)为理想脉冲响应与窗函数乘积[H,w]=freqz(h,1);

db=20*log10(abs(H));

dw = 2*pi/1000;

Rp = -(mi n(db(1:wp/dw+1))) % 检验通带波动

As = -round(max(db(ws/dw+1:501))) % 检验最小阻带衰减

%绘图

n=0:N-1;

plot(w*fs/(2*pi),db);title('幅度响应(单位：dB)');

grid;

axis([0 50 -100 10]);

xlabel('频率(单位：Hz)');

ylabel('分贝')

set(gca,'XTickMode','ma nual','XTick',[0,3,5,50])

set(gca,'YTickMode','ma nual','YTick',[-50,0])

set(gca,'YTickLabelMode','ma nual','YTickLabels',['50';' 0'])

9、要求：FIR高通数字滤波器指标为:

p =0.5二数字通带截止频率(弧度)

?，s=0.3二数字阻带截止频率(弧度)

Rp =1dB 通带衰减(dB)

A s =40d

B 阻带衰减(dB)

因为衰减为40dB,所以选择汉宁窗。

过渡带宽为Wp —Ws=0.2 n，由公式N > 6.2 n - 0.2 n =31,所以N=32。程序如下：

N=32;

Wc=pi/2;

wc=Wc/pi;%频率归一化

h=fir1(N,wc, 'high', Hanning(N+1));

[H,m]=freqz(h,[1],1024,'whole'); % 频率响应

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=a ngle(H);

subplot(2,2,1)

n=0:N;

stem( n,h,'.')

axis([0 N -0.1 0.3])

hold on

n=0:N-1;

x=zeros(N);

plot (n ,x,'-')

hold off

xlabel(' n')

ylabel('h( n)')

title('实际低通滤波器的h(n)')

subplot(2,2,2)

plot(m/pi,db)

axis([0 1 -100 0])

xlabel('w/pi') ylabel('dB') title('幅频衰减特性') grid on

subplot(2,2,3) plot(m,pha) hold on n=0:7;

x=zeros(8);

plot (n ,x,'-')

hold off axis([0 3.15 -4 4]) xlabel('频率(rad)') ylabel('相位(rad)') title('相频特性')

subplot(2,2,4) plot(m,mag) axis([0 6.15 0 1.5]) xlabel('频率W(rad)') ylabel('幅值') title('幅频特性')

10、要求如下：低端阻带截止频率wls = 02 pi ；

低端通带截止频率wlp 二0.35 pi ;

咼端通带截止频率whp = 0.65 pi ；

高端阻带截止频率whs =0.8“ pi ；用Matlab实现的程序为：

wls = 0.2*pi;

wlp = 0.35*pi;

whp = 0.65*pi;

wc = [wlp/pi,whp/pi];

B = wlp-wls;

N = ceil(8/0.15);

n=0:N-1;

win dow= hannin g(N);

[h1,w]=freqz(wi ndow,1);

figure(1);

stem(wi ndow);

axis([0 60 0 1.2]);

grid;

xlabel(' n');

figure(2);

plot(w/pi,20*log(abs(h1)/abs(h1(1))));

axis([0 1 -350 0]);

grid;

xlabel('w/pi');

ylabel('幅度(dB)');

hn = fir1(N-1,wc, ha nning (N));

[h2,w]=freqz(h n,1,512);

figure(3);

stem( n,h n);

axis([O 60 -0.25 0.25]);

grid;

xlabel(' n');

ylabel('h( n)');

figure(4);

plot(w/pi,20*log(abs(h2)/abs(h2(1)))); grid;

xlabel('w/pi');

ylabel('幅度(dB)');