数字信号处理结课论文

基于MATLAB的数字滤波器设计

摘要

数字滤波器的实现是数字信号处理中的重要组成部分，设计过程较为复杂，牵涉到模型逼近、指标选择、计算机仿真、性能分析及可行性分析等一系列的工作，本文从设计原理以及数学软件matlab出发阐述数字滤波器的设计原理与方法。

关键词：MATLAB，数字信号处理，数字滤波器

1绪论

数字滤波器是数字信号处理的重要应用，21世纪数字滤波器及其相关技术广泛的应用于通信、电子、自动控制等领域，是一种有效的抑制噪音、提取有用信号的方法。它本身可以用硬件实现，也可以通过软件来实现，还可以通过专用的DSP处理器配合相应的软件，即软硬结合的方式来实现。数字滤波器可以分为有限冲击响应（FIR）和无限冲激响应(IIR)两种。通过MATLAB程序，实现输入相应技术指标及滤波器模型，输出相应数字滤波器的参数的功能。

2 无限长单位脉冲响应滤波器（IIR）的设计

根据模拟滤波器设计数字滤波器，就是通过已知的模拟滤波器系统的系统函数H(s)来设计数字滤波器的系统函数H(z)，主要是通过脉冲响应不变法，或双线性变换法完成S平面到Z平面的转换。通过典型的模拟滤波器(诸如：巴特沃斯滤波器、切比雪夫滤波器等)可以实现一定参数要求的数字滤波器。

根据已有的数字滤波器设计不同参数或者不同频带通断类型的数字滤波器。例如已知数字低通滤波器的模型，通过变量代换得到不同截止频率的数字低通滤波器，或通过已知低通滤波器的模型设计高通、高阻、带通、带阻滤波器，这里主要是通过来完成相应的变量代换来实现滤波器类型的变换和参数的变换。

3 有限长单位脉冲响应滤波器(FIR)的设计

IIR滤波器可用于较少的阶数达到所要求的幅度特性，且实现时所需的运算次数及存储单元都很少，十分适合于对于相位特性没有严格要求的场合，如果对相位特性有要求，这时选用FIR滤波器较好。

3.1 窗函数法

从时域出发，把理想的窗口函数h d(n)截取成有限长的，以此h(n)来逼近理想的窗口函数h d(n)，从而频率响应H(jw)也近似于理想的频率响应H d(jw)。主要窗函数有：矩形窗，汉宁(Hanning)窗，汉明(Hamming)窗，布莱克曼(Blackman)窗，凯塞(Kaiser)窗。3.2 频率采样法

从频率出发，对理想的频率响应H d(e jw)加以等间隔采样，以此H d(k)实际FIR滤波器的频率特性的离散样本H(k)。

4 利用MATLAB具体实现数字滤波器

4.1 MATLAB中IIR数字滤波器设计的设计函数

信号处理工具箱提供的IIR经典设计方法是基于经典的低通滤波器到具有相同性能指标的数字滤波器的变换。基本原理就是先根据滤波器的技术指标设计出相应的模拟滤波器，然后再将设计好的模拟滤波器变换为满足指标的数字滤波器。利用表1的设计函数，可以很容易地产生任何阶数的高通、低通、带通、带阻数字滤波器。

表1 滤波器设计函数

表中，n为滤波器阶次，Wn为滤波器的归一化截止频率；函数默认为低通或带通滤波器。b,a分别为滤波器传递函数的分子和分母的系数向量；z,p,k分别为滤波器的零点、极点和增益。options为滤波器类型参数，high为高通滤波器，截止频率为Wn；stop为带阻滤波器，截止频率为Wn = [W1,W2]。

4.2 脉冲响应不变法设计数字滤波器

调用格式：[bz,az] = impinvar(b,a,Fs)，再给定模拟滤波器参数b,a和取样频率Fs的前提下，计算数字滤波器的参数。两者的冲激响应不变，即模拟滤波器的冲激响应按Fs取样后等同于数字滤波器的冲激响应。

4.3 利用双线性变换法设计数字ButterWorth滤波器

调用格式：[bz,az] = bilinear[b,a,Fs]，根据给定的分子b,分母系数a和取样频率Fs，根据双线性变换将模拟滤波器变换成离散滤波器，具有分子系数向量bz和分母向量az。

4.4 IIR设计实例

利用脉冲响应不变法设计一个巴特沃斯数字低通滤波器，要求：ωp=0.25π，R p=1dB，

ωs=0.4π，A s=15dB，滤波器采样频率F s=2000Hz。

程序如下：

wp = 0.25*pi;

ws=0.4*pi;

Rp=1;As=15;

ripple=10^(-Rp/20); %滤波器的通带衰减对应的幅度值

Attn=10^(-As/20); %滤波器的阻带衰减对应的幅度值

%转换为模拟滤波器的技术指标

Fs=2000;T=1/Fs;Omgp=wp*Fs;Omgs=ws*Fs;

%模拟原型滤波器的计算

[n,Omgc] = buttord(Omgp,Omgs,Rp,As,’s’) %计算阶数n和截止频率[z0,p0,k0] = buttap(n); %设计归一化的巴特沃斯模拟滤波器原型

ba1 = k0*real(poly(z0)); %求原型滤波器的系数b

aa1 = real(poly(p0)); %求原型滤波器的系数a

[ba,aa] = lp2lp(ba1,aa1,Omgc); %变换为模拟低通滤波器

%用脉冲响应不变法计算数字滤波器的系数

[bd,ad] = impinvar(ba,aa,Fs)

[C,B,A] = dir2par(bd,ad)

[H,w] = freqz(bd,ad);

dbH = 20*log10((abs(H)+eps)/max(abs(H));

subplot(2,2,1);plot(w/pi,abs(H));

ylable(‘|H|’);title(‘幅度响应’);axis([0,1,0,1.1]);

set(gca,’XTickMode’,’manual’,’XTick’,[0,0.25,0.4,1]);

set(gca,’YTickMode’,’manual’,’YTick’,[0,Attn,ripple,1]);grid

subplot(2,2,2);plot(w/pi,dbH);title(‘幅度响应(dB)’)；

ylabel(‘dB’);xlable(‘频率’);axis([0,1,-40,5]);

set(gca,’XTickMode’,’manual’,’XTick’,[0,0.25,0.4,1]);

set(gca,’YTickMode’,’manual’,’YTick’,[-50,-15,-1,0]);grid

subplot(2,2,4);zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title(‘零极点图’);

自定义dir2par函数：

function I=cplxcomp(p1,p2)

I=[];

for j=1:length(p2)

for i=1:length(p1)

if(abs(p1(i)-p2(j))<0.0001)

I=[I,i];

end

end

end

I=I';

function [C,B,A]=dir2par(num,den)

M=length(num);

N=length(den);

[r1,p1,C]=residuez(num,den);

p=cplxpair(p1,10000000*eps);

I=cplxcomp(p1,p);

r=r1(I);

K=floor(N/2);

B=zeros(K,2);

A=zeros(K,3);

if K*2==N;

for i=1:2:N-2;

Brow=r(i:1:i+1,:);

Arow=p(i:1:i+1,:);

[Brow,Arow]=residuez(Brow,Arow,[]);

B(fix(i+1)/2,:)=real(Brow);

A(fix(i+1)/2,:)=real(Arow);

end

[Brow,Arow]=residuez(r(N-1),p(N-1),[]);

B(K,:)=[real(Brow),0];

A(K,:)=[real(Arow),0];

else

for i=1:2:N-1;

Brow=r(i:1:i+1,:);

Arow=p(i:1:i+1,:);

[Brow,Arow]=residuez(Brow,Arow,[]);

B(fix(i+1)/2,:)=real(Brow);

A(fix(i+1)/2,:)=real(Arow);

end

end

程序运行结果如下：

n = 5

Omgc = 103.2016

bd = 0.0072 0.0362 0.0725 0.0725 0.0362 0.0072 ad = 1.0000 -1.9434 1.9680 -1.0702 0.3166 -0.0392 sos =

1.0000 0.9956 0 1.0000 -0.3193 0

1.0000

2.0072 1.0072 1.0000 -0.6984 0.2053

1.0000 1.9972 0.9973 1.0000 -0.9257 0.5976

g = 0.0072

频率特性如图1所示：

0.250.4

1

00.1778

0.8913

1|H |

幅度响应

0.250.4

1

-10

1

φ

相位响应

-15

-1

0幅度响应(dB)

d B

频率(π)

-1-0.500.51

Real Part

I m a g i n a r y P a r t

零极点图

图1 频率特性曲线

4.5 MATLAB 中FIR 数字滤波器的设计的函数

在MATLAB 中产生窗函数十分简单： 根据长度n 产生一个相应类型窗函数的类型：

（1）矩形窗(Rectangle Window)：调用格式：w = boxcar(n)。 （2）三角窗(Triangular Window)：调用格式：w = triang(n)。 （3）汉宁窗(Hanning Window)：调用格式：w = hanning(n)。 （4）哈明窗(Hamming Window)：调用格式：w = hamming(n)。 （5）布莱克曼窗(Blackman Window)：调用格式：w = blackman(n)。 （6）凯塞窗(Kaiser Window)：调用格式：w = kaiser(n,beta)。

基于窗函数的FIR 滤波器设计利用MATLAB 提供的函数firl 来实现。 调用格式：firl(n,Wn,’ftype’,Window)，n 为阶数，Wn 是截止频率（如果输入是形如[W1 W2]的矢量时，本函数将设计带通滤波器，其通带为W1<ω

是滤波器的类型(低通省略该参数，高通high，带阻stop)，Window是窗函数。

格式中：n为FIR滤波器的阶数，对于高通、带阻滤波器，n取偶数。

Wn为滤波器的截止频率，0~1；对于带通、带阻滤波器，Wn=[W1,W2]，且W1

‘ftype’为滤波器类型，默认为低通滤波器或带通滤波器。其选择如下：

‘high’为高通滤波器，’stop’为带阻滤波器，’DC-1’使多带的第一频带为通带，’D C-0’使多带的第一频带为阻带。

Window为窗函数，列向量，其长度为n+1。

4.6 FIR设计实例

用矩形窗设计一个FIR数字低通滤波器，要求：N=64，截止频率wc=0.4，描绘理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。

程序如下：

wc=0.4*pi;

N=64;n=0:N-1;

hd=ideal_lp(wc,N);

windows=(boxcar(N))';

b=hd.*windows;

[H,w]=freqz(b,1);

dbH=20*log10(abs(H)+eps)/max(abs(H));

subplot(2,2,1),stem(n,hd);

axis([0,N,1.1*min(hd),1.1*max(hd)]);title('理想脉冲响应');

xlabel('n');ylabel('hd(n)');

subplot(2,2,2);stem(n,windows);

axis([0,N,0,1.1]);title('窗函数特性');

xlabel('n');ylabel('wd(n)');

subplot(2,2,3);stem(n,b);

axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(2,2,4);plot(w/pi,dbH);

axis([0,1,-80,10]);title('幅度频率响应');

xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})');

set(gca,'XTickMode','manual','XTick',[0,wc/pi,1]);

set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid

程序运行结果如图2所示：

20

40

60

理想脉冲响应

n

h d (n )

020

4060

窗函数特性

n

w d (n )

20

40

60

实际脉冲响应

n

h (n )

0.4

1

-50

-20

-3

0幅度频率响应

频率（单位：π）

H (e j ω)

图2 FIR 数字低通滤波器

总结

本文结合数字信号处理技术与计算机技术，从理论上分析了各种数字滤波器的设计思路，以及介绍了利用MATLAB进行各种数字滤波器的设计的基本方法。利用MATLAB 提供的各种库函数，及仿真工具大大的简化了设计步骤，使设计更方便、灵活。体现了电脑辅助设计软件的优越性，为以后的数字滤波器设计学习与工作提供了新思路。

参考文献

[1] chi-Tsong chen.数字信号处理频谱计算与滤波器设计[M].北京:电子工业出版社，2002.

[2] 方勇.数字信号处理——原理与实践[M].北京:清华大学出版社，2005.

[3]吴镇杨.数字信号处理[M].北京:高等教育出版社,2005.

[4]刘波，文忠，曾涯.MATLAB信号处理[M].北京:电子工业出版社，2002.