课程设计基于双线性变换法的IIR数字低通滤波器设计.docx

课程设计基于双线性变换法的IIR数字低通滤波器设计.docx

《课程设计基于双线性变换法的IIR数字低通滤波器设计.docx》由会员分享，可在线阅读，更多相关《课程设计基于双线性变换法的IIR数字低通滤波器设计.docx（19页珍藏版）》请在冰豆网上搜索。

课程设计基于双线性变换法的IIR数字低通滤波器设计

课程设计---基于双线性变换法的IIR数字低通滤波器设计

课程设计

题目：

基于双线性变换法的IIR数字低通滤波器设计

姓名：

院系：

电气信息工程学院

专业班级：

电子信息工程11-02

学号：

541101030206

指导教师：

成绩：

时间：

2014年6月9日至2014年6月13日

课程设计任务书

题目基于双线性变换法的IIR数字低通滤波器设计

专业、班级电子信息工程11级2班学号541101030206姓名冯慧琦

主要内容、基本要求、主要参考资料等：

主要内容：

首先依据给定的性能指标，采用双线性变换法设计IIR数字低通滤波器；然后利用MATLAB软件的wavread函数读取.wav格式的语音信号，并利用所设计的滤波器对音频信号进行滤波处理，画出滤波前后信号的时域波形及频谱；最后回放语音信号，分析滤波前后的语音变化。

基本要求：

1、滤波器技术指标为：

fp=3000Hz;Ap=2dB;fs=4000Hz;As=45dB

2、采用双线性变换法设计IIR数字低通滤波器；

3、掌握利用wavread函数读取.wav格式语音信号的方法；

4、对语音信号进行滤波，并画出滤波前后信号的时域波形及频谱；

5、回放语音信号，分析滤波前后的语音变化。

主要参考资料：

1、从玉良.数字信号处理原理及其MATLAB实现[M].北京：

电子工业出版社.2009.7

2、胡广书.数字信号处理理论、算法与实现[M].北京：

清华大学出版社.2003,8

完成期限：

2010.6.24—2010.6.28

指导教师签名：

课程负责人签名：

2013年6月24日

摘要

根据IIR滤波器的特点，在MATLAB坏境下用双线性变换法设计IIR数字滤波器。

利用MATLAB设计滤波器，可以随时对比设计要求和滤波器特性调整参数，直观简便，极大的减轻了工作量，有利于滤波器设计的最优化。

关键词：

双线性变换法，数字滤波器，MATLAB，IIR

1.概述

数字滤波器对信号滤波的方法是：

用数字计算机对数字信号进行处理，处理就是按照预先编制的程序进行计算。

数字滤波器的原理如图所示，它的核心是数字信号处理器。

如果采用通用的计算机，随时编写程序就能进行信号处理的工作，但处理的速度较慢。

如果采用专用的计算机芯片，它是按运算方法制成的集成电路，连接信号就能进行处理工作，处理的速度飞快，但功能不易更改。

如果采用可编程的计算机芯片，那么，装入什么程序机器就能具有什么功能。

这种可编程芯片的优点很多，是现代电子产品的首选。

如果是对模拟信号进行处理，则需要添加模数转换器和数模转换器。

数字滤波器是按照程序计算信号，达到滤波的目的。

通过对数字滤波器的存储器编写程序，就可以实现各种滤波功能。

对数字滤波器来说，增加功能就是增加程序，不用增加元件，不受元件误差的影响，对低频信号的处理也不用增加芯片的体积。

用数字滤波方法可以摆脱模拟滤波器被元件限制的困扰。

数字滤波器一词出现在60年代中期。

由于电子计算机技术和大规模集成电路的发展，数字滤波器已可用计算机软件实现，也可用大规模集成数字硬件实时实现。

数字滤波器是一个离散时间系统（按预定的算法，将输入离散时间信号（对应数字频率）转换为所要求的输出离散时间信号的特定功能装置）。

应用数字滤波器处理模拟信号（对应模拟频率）时，首先须对输入模拟信号进行限带、抽样和模数转换。

数字滤波器输入信号的数字频率（2π*f/fs,f为模拟信号的频率，fs为采样频率，注意区别于模拟频率）,按照奈奎斯特抽样定理，要使抽样信号的频谱不产生重叠，应小于折叠频率（ws/2=π），其频率响应具有以2π为间隔的周期重复特性，且以折叠频率即ω=π点对称。

为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换、平滑。

数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。

数字滤波器在语言信号处理、图像信号处理、医学生物信号处理以及其他应用领域都得到了广泛应用。

数字滤波器有低通、高通、带通、带阻和全通等类型。

它可以是时不变的或时变的、因果的或非因果的、线性的或非线性的。

应用最广的是线性、时不变数字滤波器，以及f.i.r滤波器。

2.系统总体设计方案

2.1设计原理与步骤

2.1.1设计原理

滤波器的种类很多，从功能上可分为低通、高通、带通和带阻滤波器，每一种又有模拟滤波器和数字滤波器两种形式。

如果滤波器的输人和输出都是离散时间信号，则该滤波器的冲击响应也必然是离散的，这种滤波器称之为数字滤波器。

数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。

数字滤波器也是具有一定传输选择特性的数字信号处理装置，其输入、输出均为数字信号，实质上是一个由有限精度算法实现的线性时不变离散系统。

IIR数字滤波器采用递归型结构，即结构上带有反馈环路。

IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。

数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应（IIR）数字滤波器和有限长冲激响应（FIR）数字滤波器。

IIR数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型

来实现，其差分方程为：

（2-1）

（2-2）

设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H（z），使其频率响应H（z）满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。

本次课程设计采用双线性变换法。

2.1.2设计步骤：

（1）将给出的数字滤波器的技术指标转换为模拟滤波器的技术指标；

（2）根据转换后的技术指标设计模拟低通滤波器H（s）；

（3）在按一定规则将H（s）转换为H（z）；

若所设计的数字滤波器是低通的，那么上述设计工作可以结束，若所设计的是高通、带通或者带阻滤波器，那么还有步骤：

（4）将高通、带通或者带阻数字滤波器的技术指标先转化为低通滤波器的技术指标，然后按上述步骤

（2）设计出模拟低通滤波器H（s），再将H（s）转换为所需的H（z）。

2.2设计方案

IIR数字滤波器是一种离散时间系统，其系统函数为

（2-3）

假设M≤N，当M＞N时,系统函数可以看作一个IIR的子系统和一个（M-N）的FIR子系统的级联。

IIR数字滤波器的设计实际上是求解滤波器的系数ak和bk，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。

如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。

实现IIR数字滤波器的设计有双线性变换法和脉冲响应不变法两种基本方案，现在就对两种基本方案的优劣进行具体论证，从而说明选择方案一的理由。

方案一：

双线性变换法设计IIR数字滤波器

双线性变换法主要是采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。

也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。

这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图2-1所示。

图2-1双线性变换的映射关系

为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现

（2-4）

式中,T仍是采样间隔。

当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。

将式（2-4）写成

将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得（2-5）

再将S1平面通过以下标准变换关系映射到Z平面

z=es1T

从而得到S平面和Z平面的单值映射关系为：

（2-6）

（2-7）

式（2-6）与式（2-7）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换

式（2-5）与式（2-6）的双线性变换符合映射变换应满足的两点要求。

首先,把z=ejω，可得

（2-8）

即S平面的虚轴映射到Z平面的单位圆。

其次，将s=σ+jΩ代入式（2-8），得

因此

由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。

也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。

因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。

方案二：

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

脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h（n）模仿模拟滤波器的冲激响应ha（t），即将ha（t）进行等间隔采样，使h（n）正好等于ha（t）的采样值，满足h（n）=ha（nT）式中,T是采样周期。

如图2-2所示。

如果令Ha（s）是ha（t）的拉普拉斯变换，H（z）为h（n）的Z变换，利用采样序列的Z变换与模拟信号的拉普拉斯变换的关系得

（2-9）

可看出，脉冲响应不变法将模拟滤波器的S平面变换成数字滤波器的Z平面，这个从s到z的变换z=esT是从S平面变换到Z平面的标准变换关系式。

图2-2脉冲响应不变法的映射关系

由（2-9）式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为

（2-10）

这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。

正如采样定理所讨论的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率以内时，即

（2-11）

才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真，即

（2-12）

但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真，如图2-3所示。

这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。

当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。

这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。

图2-3脉冲响应不变法中的频响混叠现象

对某一模拟滤波器的单位冲激响应ha（t）进行采样，采样频率为fs，若使fs增加，即令采样时间间隔（T=1/fs）减小，则系统频率响应各周期延拓分量之间相距更远，因而可减小频率响应的混叠效应。

两种方案优缺点比较：

脉冲响应不变法的最大缺点是有频率响应的混叠效应。

所以，脉冲响应不变法只适用于限带的模拟滤波器（例如，衰减特性很好的低通或带通滤波器），而且高频衰减越快，混叠效应越小。

至于高通和带阻滤波器，由于它们在高频部分不衰减，因此将完全混淆在低频响应中。

如果要对高通和带阻滤波器采用脉冲响应不变法，就必须先对高通和带阻滤波器加一保护滤波器，滤掉高于折叠频率以上的频率，然后再使用脉冲响应不变法转换为数字滤波器。

当然这样会进一步增加设计复杂性和滤波器的阶数。

双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。

这是因为S平面与Z平面是单值的一一对应关系。

S平面整个jΩ轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。

这个关系如式（2-8）所示，重写如下：

上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图2-4所示。

由图2-4看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。

图2-4双线性变换法的频率变换关系

但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（2-8）及图2-2所示。

由于这种频率之间的非线性变换关系，就产生了新的问题。

首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变。

对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。

也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。

由于用脉冲响应不变法设计IIR数字滤波器，存在从S平面到Z平面的多值映射的缺点，而且还有频谱混叠现象，所以本次课程设计就采用双线性变换法来设计IIR数字滤波器。

2.3设计步骤

要求：

根据以上IIR数字滤波器设计方法，采用双线性变换法基于MATLAB设计一个巴特沃斯数字低通滤波器，其中要求Fp=3000Hz;通带最大衰减Rp=2dB;阻带最小衰减As=45dB;滤波器采样频率Fs=4000Hz.

2.3.1设计步骤

（1）根据任务,确定性能指标:

在设计低通滤波器之前,首先根据工程实际的需要确定滤波器的设计指标:

（2）根据公式Ω=2/T*tan（w/2）进行预修正，将数字滤波器指标转换成模拟滤波器设计指标，对双线性变换法一般T=2s

（3）将模拟高通滤波器指标转换为模拟低通滤波器指标。

（4）设计模拟低通原型滤波器。

用模拟低通滤波器设计方法得到模拟低通滤波器的传输函数Ha（s）并调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。

（5）利用双线性变换法将模拟低通滤波器Ha（s）转换成数字低通滤波器H（z）.

2.3.2程序流程框图

开始

↓

读入数字滤波器技术指标

↓

将指标转换成归一化模拟低通滤波器的指标

↓

设计归一化的模拟低通滤波器阶数N和1db截止频率

↓

模拟域频率变换，将G（P）变换成模拟低通滤波器H（s）

↓

用双线性变换法将H（s）转换成数字低通滤波器H（z）

↓

输入信号后显示相关结果

↓

结束

2.3.3运行结果及分析

低通滤波器的幅频响应

根据滤波器技术指标要求：

fp=3000HZ，Ap=2dB，fs=4000HZ，As=45dB.采用双线性变换法设计IIR数字低通滤波器，滤波器幅频响应和幅频响应（dB）曲线如2-5图所示：

其中，横坐标为归一化频率0对应0Hz，1对应Fs/2.

图2-5幅频响应及幅频响应（dB）

语音信号频谱图

对原始，加噪音及滤波后的语音信号调用Matlab库函数fft（）进行傅里叶变换即可得到对应的语音信号频谱图。

由于频率特性只能限于折叠频率Fs/2以内，所以对横坐标归一化处理时：

0对应0Hz，1对应Fs/2，如图2-6所示：

其中，从上到下依次为原始语音信号频谱图，加噪声后语音信号频谱图，滤波后语音信号频谱图

图2-6语音信号频谱图

结论与展望

通过本次课程设计，不仅使我更加熟悉Matlab软件的应用，而且对数字信号处理中相关的各类滤波器的设计又进行了一次系统的复习；也使我加深了对课堂抽象概念的理解，巩固了课堂上所学的理论知识，理解和掌握了数字信号处理中的基本概念、基本原理和基本分析方法。

在课程设计过程中，通过自己搜集材料，扩大了自己的知识面，系统地掌握了编程方法和解决实际问题的一些技巧，让我对MATLAB有了更深刻的了解，对数字滤波器的设计流程有了大致的了解，掌握了一些设计滤波器的基本方法，提高了理论用于实践的能力，掌握了更多专业相关的使用知识与技能。

同时暴露出了自己许多的不足之处，在以后的学习过程中还有待继续提高。

通过学习MATLAB相关知识，发现MATLAB具有相当强大的矩阵运算和操作功能，其程序编写也比其他语言要简单。

参考文献

[1]蒋志凯.数字滤波与卡尔曼滤波[M].北京：

中国科学技术出版社，1993

[2]施阳．等．MATLAB语言工具箱[M]..西安：

西北工业大学出版社,1999

[3]楼顺天、李博菡.基于MATLAB的系统分析与设计－信号处理[M].西安：

西安电子科技大学出版社，1998.

[4]蒙以正.MATLAB5.X应用与技巧[M].北京：

科学出版社，1999.

[5]胡广书.数字信号处理：

理论、算法与实现[M].北京：

清华大学出版社，西安：

1997.

[6]倪养华．数字信号处理与实现．上海：

上海交通大学出版社,1998  

[7]施阳.MATLAB语言工具箱．西安：

西北工业大学出版

[8]边肇祺.模式识别（第二版）.北京：

清华大学出版社,1988，25～35 

[9]胡广书.数字信号处理—理论、算法与实现[M].北京:

清华大学出版社,1997.

附录

MATLAB程序

clear;clc;closeall;

[S,Fs,Bit]=wavread（'SHE.wav'）;%读取音频信号

sound（S,Fs）;

fnoise1=10000;%频率为fnoise的正弦噪声信号

fnoise2=50;%频率为fnoise的正弦噪声信号

N=length（S）;

T=N/Fs;

t=T/N:

T/N:

T;

forj=1:

N

Noise1（j）=sin（2*pi*fnoise1*t（j））;

Noise2（j）=sin（2*pi*fnoise2*t（j））;

end

%对于低通滤波器：

S1=S+0.1*Noise1';

fp=3000;

fs=4000;

Rp=2;

Rs=45;

%求数字频率

wp=2*pi*fp/Fs;

ws=2*pi*fs/Fs;

%预畸变处理

Wp=2*Fs*tan（wp/2）;

Ws=2*Fs*tan（ws/2）;

%开始设计模拟低通滤波器

[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）;

[B,A]=butter（N,Wn,'s'）;

[b,a]=bilinear（B,A,Fs）;%impinvar冲激响应不变法

[h,w]=freqz（b,a,1000,'whole'）;

h=h（1:

501）;

w=w（1:

501）;

mag=abs（h）;

db=10*log10（（mag+eps）/max（mag））;

subplot（1,2,1）;plot（w/pi,mag）;

title（'幅频响应'）;

axis（[0,1,0,1.1]）;

subplot（1,2,2）;plot（w/pi,db）;

title（'幅频响应（dB）'）;

axis（[0,1,-40,1]）;

set（gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]）;

set（gca,'YTickMode','manual','YTick',[-40,-20,-3,0]）;grid;

X1=abs（fft（S））;

Z1=（0:

150000）/150000;

subplot（3,1,1）;plot（Z1,X1（1:

150001））;

X2=abs（fft（S1））;

Z2=（0:

150000）/150000;

subplot（3,1,2）;plot（Z2,X2（1:

150001））;

S2=filter（b,a,S1）;

sound（S2,Fs）;

X3=abs（fft（S2））;

Z3=（0:

150000）/150000;

subplot（3,1,3）;plot（Z3,X3（1:

150001））;