基于MATLAB的滤波器设计

最大衰减不超过0.3dB，抽样频率为20kHz。

其MATLAB程序见附录，频率响应如图6

图5 Butterworth带阻滤波器增益图

图6 ChebyshevⅡ高通滤波器的频率响应

可以看到经过离散采样、数字滤波后分离出了83.3hz的频率分量。之所以选取上面的叠加信号

作

为原始信号，是由于在实际工作中是要对已经经过差分滤波的信号进一步做带通滤波，信号的各分量基本

同一致，可以反映实际的情况。本例设计的滤波器已在实际工作中应用，取得了不错的效果。

结论

采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是Matlab工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。MATLAB 信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具, 在设计数字滤波器时, 善于应用MATLAB进行辅助设计, 能够大大提高设计效率。

参 考 文 献

[1] 陈德树. 计算机继电保护原理与技术.北京：水利电力出版社，1992. [2] 蒋志凯. 数字滤波与卡尔曼滤波.北京：中国科学技术出版社，1993

[3] 楼顺天、李博菡. 基于MATLAB的系统分析与设计－信号处理.西安：西安电子 [4] 董长虹等. MATLAB信号处理与应用.北京：国防工业出版社，2005 [5] [美] M.H.海因斯 著，张建华等 译.数字信号处理.北京：科学出版社，2002 [6] 张葛祥，李 娜. MATLAB仿真技术与应用.北京：清华大学出版社，2002

附录

1、ChebyshevⅡ高通滤波器的频率响应 ws=2*pi*1000; ws1=ws*2*pi; wp=2*pi*1400; wp1=wp*2*pi; rp=0.3; rs=15; fs=20000;

[N,Wn]=cheb2ord(wp1,ws1,rp,rs,’s’); [z,p,k]=cheb2ap(N,rs); [A,B,C,D]=zp2ss(z,p,k);

[At,Bt,Ct,Dt]=lp2hp(A,B,C,D,Wn); [At1,Bt1,Ct1,Dt1]=bilinear(At,Bt,Ct,Dt,fs); [num,den]=ss2tf(At1,Bt1,Ct1,Dt1); freqz(num,den); [H,W]=freqz(num,den); plot(W*fs/(2*pi),abs(H));grid; xlabel(‘频率/Hz’);

ylabel(‘幅值’);

2、带阻滤波器幅频特性实现 n=55-1; w=[0.4 0.6];

beta=0.1102*(60-8.7); kaiw=Kaiser(n+1,beta); b=fir1(n,w,’stop’,kaiw); [h,w]=freqz(b,1,512,2); plot(w,20*log10(abs(h)));grid; xlabel(‘频率（归一化）’); ylabel((‘幅度dB’));

3、Chebyshev 低通（I 型）滤波器 T=0.1; FS=1/T; fp=3;fs=4; Rp = 1; As = 15;

wp=fp/FS*2*pi; ws=fs/FS*2*pi; OmegaP = (2/T)*tan(wp/2); OmegaS = (2/T)*tan(ws/2); ep = sqrt(10^(Rp/10)-1); A = 10^(As/20); OmegaC = OmegaP; OmegaR = OmegaS/OmegaP; g = sqrt(A*A-1)/ep;

N = ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1))); [z,p,k] = cheb1ap(N,Rp); a = real(poly(p)); aNn = a(N+1); p = p*OmegaC; a = real(poly(p)); aNu = a(N+1); k = k*aNu/aNn; b0 = k;

B = real(poly(z)); b = k*B;

[bz,az] = bilinear(b,a,FS); %freqz(bz,az,200,FS,'whole'); H=freqz(bz,az,200,'whole'); plot(abs(H)); xlabel(‘频率/Hz’); ylabel(‘幅值’);

4、直接设计数字滤波器

FS=1;

[n,Wn]=buttord(0.3/(FS/2),0.2/( FS /2),1,25); [b,a]=butter(n,Wn,'high'); freqz(b,a,512, FS);

5、ButterWorth 低通滤波器 T=0.1; FS=1/T; fp=3;fs=4; wp=fp/FS*2*pi; ws=fs/FS*2*pi; Rp = 1; As = 20;

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

N = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS))); OmegaC = OmegaP/((10^(Rp/10)-1)^(1/(2*N))); [z,p,k] = buttap(N); p = p*OmegaC; k = k*OmegaC^N; B = real(poly(z)); b0 = k; cs = k*B; ds = real(poly(p)); [b,a] = bilinear(cs,ds,FS); freqz(b,a,512,FS);

6、Butterworth带阻滤波器增益图 function [g,w] = gain(num,den)

% Computes the gain function in dB of a % transfer function at 256 equally spaced points % on the unit circle

% Numerator coefficients are in vector num % Denominator coefficients are in vector den % Frequency values are returned in vector w % Gain values are returned in vector g w = 0:pi/255:pi; h = freqz(num,den,w); g = 20*log10(abs(h));

% Design of a Butterworth Bandstop Digital Filter Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50; % Estimate the Filter Order

[N1, Wn1] = buttord(Wp, Ws, Rp, Rs); % Design the Filter

[num,den] = butter(N1,Wn1,'stop'); % Display the transfer function

disp('Numerator Coefficients are ');disp(num); disp('Denominator Coefficients are ');disp(den); % Compute the gain response [g, w] = gain(num,den); % Plot the gain response plot(w/pi,g);grid axis([0 1 -60 5]);

xlabel('\\omega /\\pi'); ylabel('Gain in dB');

title('Gain Response of a Butterworth Bandstop Filter');