数字滤波器研究课程论文

研 究 生 课 程 论 文
(2010-2011学年第一学期)
数字滤波器的研究
研究生：
中文题目：数字滤波器的研究
研究生姓名
摘要：无限冲激响应（IIR）数字滤波器和有限冲激响应（FIR）数字滤波器广泛应用于数字信号处理系统中。

FIR 滤波器具有很好的线性相位特性，因此线性相位特性在工程实际中如数据通信、图像处理等应用领域具有比较重要的意义。

IIR 数字滤波器具结构简单和实现容易的特点, 这对于相位要求不是特别严格的场合有很重要的意义。

利用MATLAB 实现IIR 低通数字滤波器的设计,并对被干扰的信号进行滤波,达到了预期结果。

关键词: IIR；FIR；线性相位；数字滤波器；MATLAB
1 引言
数字滤波器是现代数字信号处理的不可或缺的一部分，数字滤波器按其响应形式可分为无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器和两种。

当系统无严格相位要求时，IIR 数字滤波器与FIR 数字滤波器相比，可用较低的阶数获得较高的选择性[1]。

对于相同的设计指标，FIR 滤波器所要求的阶数比IIR滤波器高5～10 倍，而且信号的延迟也较大[2]。

MATLAB 作为一种矩阵运算为基础的交互式程序语言，是进行科学研究常用且必不可少的工具。

MATLAB着重针对科学计算、工程计算和绘图的需求。

它用解释方式工作，键入程序立即得出结果，人机交互性能好，为科学人员所乐于接受。

MATLAB 提供了丰富的函数，其中fir1 函数实现了加窗线性相位FIR 滤波器设计的经典方法，fir1 主要用于常用的标准通带滤波器设计，包括低通、带通、高通和带阻数字滤波器[3]。

一般来讲，数字滤波器的设计要经过三步：确定设计指标、模拟逼近和数字转换。

通常在设计滤波器之前，先根据具体的应用确定一些技术指标，
然后就可以根据数学知识和滤波器的基本原理提出一个滤波器的模型来逼近给定的指标，逼近的结果通常是得到以差分方程或脉冲响应描述的滤波器，最后可以根据这个描述用硬件或软件实现，至此完成一个滤波器设计的全过程[3]。

2.数字滤波器的设计
2.1模拟椭圆滤波器对加噪语音信号的滤波：
设计模拟滤波器的一般方法是：首先根据技术指标确定滤波器的传输函数H(s）,然后综合电路网络实现该传递函数。

而设计滤波器传输函数的关键是找到逼近函数，在各种滤波器响应中，椭圆函数响应最优越[4-6]。

首先，打开的是一个摩托车引擎信号的波形文件，它是一个频率范围大致分布在5000-6000Hz的一个低频语音信号。

对于该语音信号，施加高斯白噪声，设计一个1dB截止频率6000Hz 最小阻带衰减为100dB的椭圆低通滤波器。

其原始信号波形、加噪声后信号波形与滤波后信号波形分别图1和图2：
图1.加噪声后语音信号波形与频谱
图2.原始语音信号与滤波后语音信号波形与频谱对比
通过matlab仿真可以发现，原始语音信号基本可以从高斯白噪声中恢复，但是滤波后的语音信号比原始语音信号有明显的波形延迟，这是模拟椭圆低通滤波器不足的地方。

2.2IIR滤波器的结构与设计
这里采用的是直接II型结构实现IIR 滤波器，如式（1）所示：
（1）
此结构便于准确地实现滤波器的零点、极点，也便于调整滤波器的频率响应性能[7]；另一个优点是所需的存储单元较少，在硬件实现时甚至还可以用一个二阶节进行时分复用，进一步降低了对现场可编程逻辑阵列硬件资源的要求[8]。

IIR滤波器的系数在做强噪音下的语音增强，需要通过共振峰对语音进行
端点识别。

前两个共振峰对区别不同语音有非常重要的作用[9]；1 kHz以下基本为第一共振峰的范围,对语音感知、语意识别作用比较重要的第二共振峰基本在1 kHz之外[9]。

因此对带噪语音信号加入预处理环节,即进行数字滤波(通带下限1 kHz,Butterworth滤波器高通,阻带衰减3 dB，语音信号采样频率为8000 Hz)。

利用MATLAB提供的FDATool便可直接得到直接II型各个子系由于硬件当中对小数进行运算耗费资源比较大，需要在计算精度与速度之间优化选择[10]。

下面用双线性变换设计一个0.4dB截止频率为10KHz且在30KHz处有最小阻带衰减为50dB的数字巴特沃兹低通滤波器，其抽样率为100KHz；
设计的具体步骤：
1：先设计一个符合以上指标的模拟巴特沃兹低通滤波器，作为与要设计的数字滤波器的参照滤波器。

先用手算计算出滤波器所需的阶数是Ns=7, 与[Ns,wn]=buttord(wp,ws,Rp,Rs,'s');得出的结果一样。

2：双线性变换公式：
)
2
tan(
2ω
T
=
Ω（2）
其中，Ω是模拟的频率，ω是数字的归一化频率，T是积分步长。

用公式（2）预畸所求数字滤波器的的数字频率指标（数字滤波器通带截止频率0.2π，阻带截止频率0.6π），得到一个等价的模拟低通滤波器的频率指标。

3：得到原型模拟低通滤波器的指标后作频率响应，从而得模拟的传输函数H(s)，然后对传输函数进行双线性变换，从而得到数字滤波器的传输函数H(z).
结果分析：
图3.模拟原型与数字滤波器幅度相应比较
图4. 模拟原型与数字滤波器相位相应比较
由图1可以看出模拟原型低通滤波器的通带截止频率为10KHz对应着数字低通滤波器的截止频率0.2π rad与双线性变换结果一致。

由图2可以看出，相位响应就有很大差别了，这是由双线性变换的频率畸变所致。

2.3 FIR滤波器的结构与设计
对于IIR滤波器必须明确要求推导出来的传输函数是稳定的，而另一方面，对于FIR滤波器的设计，由于滤波器的传输函数是以1 /z的多项式表示的，所以FIR数字滤波器都是稳定的[11]。

利用傅里叶加窗级数法设计线性相位低通FIR滤波器（没有使用到fir1函数）设计指标：通带截止频率4rad/s，阻带截止频率6rad/s，最大通带衰减为0.2dB,最小阻带衰减为42dB，抽样率为18rad/s.
图5.由截短得到的理想FIR低通滤波器的冲击响应
图6.各种窗函数设计的低通滤波器结果分析：
Hamming窗函数的长度是30
通带波纹是Rp = 0.0394dB
最小阻带衰减时As = 52dB
Hanning窗函数的长度是29
通带波纹是Rp= 0.0711dB
最小阻带衰减时As = 44dB
Blackman窗函数的长度是52
通带波纹是Rp= 0.0027dB
最小阻带衰减时As = 75dB
均满足滤波器的性能指标。

结论
在电子系统中，具有严格的线性相位特性的FIR 数字滤波器被广泛使用，是信号处理的基本构成之一。

利用matlab的强大运算功能，基于matlab信号处理工具箱（signal processing toolbox）的数字滤波器设计法可以快速有效的设计由软件组成的常规数字滤波器，设计方便、快捷，极大的减轻了工作量。

在设计过程中可以对比滤波器特性，随时更改程序参数，以达到滤波器设计的最优化。

利用matlab设计数字滤波器在信号处理软件和微机保护中，有着广泛的应用前景。

本文通过调用simulink中的功能模块构成数字滤波器的仿真框图，编写matlab语言进行仿真，从结果可以看出各种滤波器的性能并反映实际的情况。

参考文献
[1] 李惠琼，孟令军，王宏涛.基于Nios Ⅱ的IIR数字滤波器的设计[J].通信技术,2009,42(08)
[2] 王秀敏,汪毓铎,张洋,杨世华. 通信系统中FIR数字滤波器的设计研究[J]. 通信技术, 2009,(09) .
[3梁辰基于MATLAB的FIR数字滤波器的设计[J].机械设计与制造，2010，（12）
[4] 史燕,杨小雪. 基于Matlab和Multisim的综合性实验——椭圆滤波器设计与仿真[J]. 北华航天工业学院学报, 2008,(S1) .
[5] 王靖,李永全. 数字椭圆滤波器的Matlab设计与实现[J]. 现代电子技术, 2007,(06) .
[6] 刘抒珍,童子权,任丽军,刘小红. DDS波形合成技术中低通椭圆滤波器的设计[J]. 哈尔滨理工大学学报, 2004,(05) .
[7] 杨晓琳. 循环码理论及其译码算法研究——设计距离为n 的二元BCH 码构造及其
B-M 迭代译码算法实现[D].成都：成都理工大学,2008.
[8] 王宇. BCH 编码在GPS 探空仪中的应用[J].信息安全与通信保密,2010，(07)
[9] 赵亚梅, 杜红棉, 张志杰. 基于MATLAB一种IIR数字带通滤波器的设计与仿真[J]. 微计算机信息, 2007,(13)
[10]杨世华,王秀敏,陈豪威;基于DSP Builder和FPGA的IIR滤波器设计[J],2010,(12)
[11]LAWRENCE.R.R FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation PROCEEDINGS OF THE IEEE, VOL. 63, NO. 4, APRIL 1975.
附仿真源程序：
程序1：
[y1,Fs,bits]=wavread('engine.wav',2048);
%Fs是音频采样率，bits是每个采样点包含的bit数目sound(y1,Fs,bits);
hold on
Y1=fft(y1,4096);
f=Fs*(0:4095)/4096;
y2=y1+0.2*randn(2048,1);
figure(1)
subplot(211)
plot(y2);
xlabel ('时间序号n'); ylabel ('振幅');
legend('叠加噪声后信号波形',1)
axis ([0 ,2100, -1 , 1]);
Y2=fft(y2,4096);
subplot(212)
plot(f,fftshift(abs(Y2)));
xlabel ('信号频率Hz'); ylabel ('振幅');
legend('叠加噪声后信号频谱',1)
axis ([4000,7000, 0 , 150]);
Ap=1;%通带波纹
As=100;%最小阻带衰减
fp=200;
fs=300;
Ft=20000;
wp=2*pi*fp/Ft;%归一化通带截止角频率
ws=2*pi*fs/Ft;%归一化阻带截止角频率
[n,wn]=ellipord(wp,ws,Ap,As);
[num,den]=ellip(n,Ap,As,wn);
y3=filter(num,den,y2);
Y3=fft(y3,4096);
figure(2)
subplot(211)
plot(y1,'r-');
hold on
plot(y3,'g-');
xlabel ('时间序号n'); ylabel ('振幅');
legend('原始信号波形','滤波后信号波形',1)
axis ([0 ,2100, -1 , 1]);
subplot(212)
plot(f,fftshift(abs(Y1)),'r-');
hold on
plot(f,fftshift(abs(Y3)),'g-');
xlabel ('信号频率Hz'); ylabel ('振幅');
legend('原始信号频谱','滤波后信号频谱',1)
axis ([4000,7000, 0 , 150]);
程序2：
%用双线性变换设计一个0.4dB截止频率为10KHz且在30KHz处有最小阻带衰减为50dB 的数字巴特沃兹低通滤波器，其抽样率为100KHz；
clear all;
fp=10000;
fs=30000;
Fs=100000;
Rp=0.4;
Rs=50;
wp=2*pi*fp;
ws=2*pi*fs;
[Ns,wn]=buttord(wp,ws,Rp,Rs,'s');
[bs,as]=butter(Ns,wn,'s');
%计算模拟原型滤波器的频率响应
w = [0: 200: 80000*pi];
h = freqs(bs,as,w);
%预畸所求数字滤波器的数字频率指标，得到一个等价的的模拟滤波器的指标
Wp=tan(wp/(2*Fs))
Ws=tan(ws/(2*Fs))
[Nz,Wn]=buttord(Wp,Ws,Rp,Rs,'s');
[b,a]=butter(Nz,Wn,'s');
%进行双线性变换，得到所求的数字IIR传输函数
[Bz,Az]=bilinear(b,a,0.5);
[H,omega] = freqz(Bz,Az,1024);
figure(1)
subplot (211);
plot (w/(2*pi*1000),20*log10(abs(h)));
xlabel('Frequency, KHz'); ylabel('Gain, dB');
title('模拟原型滤波器幅度响应')
subplot (212);
plot (omega/pi,20*log10(abs(H)));
xlabel('\omega/\pi'); ylabel('Gain, dB');
title('数字滤波器幅度响应')
figure(2)
subplot (211);
phase=180*angle(h)/pi;
plot (w/(2*pi*1000),phase);
xlabel('Frequency, KHz'); ylabel('phase,°');
title('模拟原型滤波器相位响应')
subplot (212);
Phase=180*angle(H)/pi;
plot(omega/pi,Phase);
xlabel('\omega/\pi'); ylabel('Phase,°');
title('数字滤波器相位响应')
程序3：
%利用窗函数发设计线性相位低通FIR滤波器
wp=2*pi*4/18; %归一化通带截止频率
ws = 2*pi*6/18; %归一化阻带截止频率
tr_width =ws-wp; %确定过度带宽
N1=ceil (6.64*pi/tr_width)+1;% 确定滤波器阶数，ceil是向上取整函数
N2=ceil (6.22*pi/tr_width)+1;% 确定滤波器阶数，ceil是向上取整函数
N3=ceil (11.12*pi/tr_width)+1;% 确定滤波器阶数，ceil是向上取整函数
n1 =[0:1:N1-1];
n2 =[0:1:N2-1];
n3 =[0:1:N3-1];
wc=(ws+wp)/2; %理想低通的截止频率
hd1=ideal_lp(wc, N1);%调用理想低通滤波器的单位冲击响应函数
hd2=ideal_lp(wc, N2);%调用理想低通滤波器的单位冲击响应函数
hd3=ideal_lp(wc, N3);%调用理想低通滤波器的单位冲击响应函数
W1=(hamming(N1))';
W2=hanning(N2)';
W3=blackman(N3)';
h1= hd1.*W1;
h2= hd2.*W2;
h3= hd3.*W3;
[db1,mag1,pha1,w1]=freqz_m(h1,[1]);%计算hamming窗低通滤波器的幅度响应[db2,mag2,pha2,w2]=freqz_m(h2,[1]);%计算hanning窗低通滤波器的幅度响应[db3,mag3,pha3,w3]=freqz_m(h3,[1]);%计算blackman窗低通滤波器的幅度响应delta_w=2*pi/18;
Rp1= - (min(db1(1:1: wp/delta_w +1))) %hamming窗实际通带波动
As1= - round(max(db1(ws/delta_w +1: 1: 501 )))% hamming窗最小阻带衰减
Rp2= - (min(db2(1:1: wp/delta_w +1))) %hanning窗实际通带波动
As2= - round(max(db2(ws/delta_w +1: 1: 501 )))% hanning窗最小阻带衰减
Rp3= - (min(db3(1:1: wp/delta_w +1))) %blackman窗实际通带波动
As3= - round(max(db3(ws/delta_w +1: 1: 501 )))% blackman窗最小阻带衰减figure(1)
stem(n1, h1,'r-s');
hold on
stem(n2, h2,'g-d' );
hold on
stem(n3, h3,'b-*' );
legend('hamming窗','hanning窗','blackman窗',1)
xlabel ('时间序号n'); ylabel ('振幅');
axis ([0,52,-0.2,0.6]);
title ('由截短得到的理想FIR低通滤波器的冲激响应');
figure(2)
plot (w1/pi, db1,'r-s' );
hold on
plot (w2/pi, db2,'g-d' );
hold on
plot (w3/pi, db3,'b-*' );
axis ([0 ,1, -150 , 10]);
xlabel ('归一化频率（单位：π）'); ylabel ('增益（单位：dB）');
legend('hamming窗','hanning窗','blackman窗',1)
title ('各种窗幅度设计的滤波器');
%自己编写的理想低通滤波器的单位冲击响应函数
function[hd]=ideal_lp(wc,N)
alpha=(N-1)/2;
n=[0:1:(N-1)];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);。