数字信号处理滤波器的设计.doc

数字信号处理滤波器的设计

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

个人收集整理 勿做商业用途
对噪声下的男女声信号进行滤波处理
摘 要：本次设计当中，我们以声音信号的滤波工作作为核心.主要的内容就是滤波器的选取过程。根据我们所学的知识，分别设计不同的FIR，IIR滤波器。设计过程中，我们不仅对滤波器的本身特性(包括四个主要的参数)进行了一些研究，还从滤波前后的时域、频域图形进行了对比。
本次课题的目的就是在针对具体实例下区分和选取最优的滤波器，并对设计的滤波器进行评价。我们在设计的过程中正是通过设计滤波器的阶数(对应占用的资源),设计后的实际(通阻带纹波、通阻带截止频率）进行考量，最终可以获得最优滤波效果的滤波器。
关键词： FIR IIR 最优的滤波器 通阻带纹波 通阻带截止频率

目 录
1 采集噪声下的男女声信.。....。.。......。。..。。.。。....。..。.。...。.....3
2 滤波器的选取问题。。...。..。。.。。..。.。。。。..。..。。.。。。。。.。.....。。。。。。3
3 滤波器的参数的设计....。。。....。。。。..。。。..。。.。。...。。。.。。...。..。。。4
4 滤波器的设计内容
4.1 FIR类型的滤波
4.1。1窗函数法...。。。..。.。。。。。。..。...。.。。。..。。。..。。。。。。。。。。.4
4.1.2 最优等纹波线性相位法..。.。.。。...。。。......。..。......。.。5
4。2 IIR类型的滤波
4.2。1 导数逼近法
4。2.2 冲激不变法
4。2.3 双线性变换法..。..。..。.。...。。...。。。。。。。。..。。。..。。..。.。8
5 滤波前后的特性对比.。..。..。。.。。。。..。...。.。...。。.。.。。...。。..。。.10
6 对工作的概括.。。。.。.。.。.....。..。..。。。.。。。.。。.。.....。.。。。。..。.。.11
7 滤波器设计的结论..。。.。。..。。..。。.。。。。。..。。..。..。.。。。.。.。。.。。.。.11
7。1滤波器的设计过程
8 参考的文献..。。.。。...。。。。..。。。.。。。。。.。.。。。..。.。..。.。....。.。。。.。12
9 附录。。。.。.。。......。.。..。...。。。。..。。。.。..。。.。...。。.。。..。...。。.。12

1. 采集噪声下的男女声信号
在半期测验的报告中我们已经简述过声音信号的采集过程了。这一次我们采集自然噪声下的语音信号。这一次，我们仍采用wavrecord函数采集声音信号,在这之后将会生成一个wav文件。
同样的，我们可以通过半期测验的工作成果得到初始状态下时域和频域上的结果。得到的图形如下：

2. 滤波器的选取问题
就我们现阶段学习的内容而言，我们选取数字滤波器考虑的第一个问题就是该采用IIR滤波器还是FIR滤波器。
根据程佩青《数字信号处理教程》P369页描述的对FIR和IIR滤波器的比较

可以看出，IIR和FIR在各自的设计过程中都有一定的针对性，具体的情况下两者分别会展现出各自的优势.
就上述的语音信号而言，我们并不能通过IIR和FIR的一些异同就直接判定哪一种更好，该选取使用哪一种。还是要通过滤波之后展现出来的滤波效果进行评判。
3. 滤波器的参数的设计
从上述的时域和频域图形来看，我们可以知道人声(我们认为的有效声音)分布在频率较低的部分，引入的噪声则分布在频率较高的部分。我们设计的滤波器应该是一个低通滤波器.
设计低通滤波器，我们还给出了如下的参数指标：

4。滤波器的设计内容
简单来说,我们学习过的设计方法可以用下表来概括：

类型
滤波器

I

II

III
FIR
窗函数设计法
频率采样法
最优等纹波线性相位法
IIR
导数逼近法
冲激不变法
双线性变换法
我们知道:
1) 频率采样法能够设计具有线性相位的FIR滤波器,但是它的主要优点在于大多频率采样为零时得到的有效频率采样结构。
2) 双线性变换法一般适用于设计满足高通和带通条件的IIR滤波器.
基于以上的原因，我们不会考虑上述两种设计方法。
下面，我们对此外的几种方法对滤波器进行设计。对应的，我们会给出仿真的结果。
4.1 FIR类型的滤波
4。1。1 窗函数法
窗函数滤波是FIR滤波器设计中最常用也最简单的一种方式，它的设计原理也很简单。通过与的傅里叶变换关系:
(其中）
可以知道，这里的单位冲激响应在时间上是无限的，一定要在某一点上截断，才可以得到长度为M的FIR滤波器。这里，我们就用到了窗函数，截取得到一段长度的。
在常用的窗函数中，凯撒窗可以提供60db的衰减，我们先选取凯撒窗，并利用matlab内部的工具箱设计得到下列有关幅频响应、相频响应、冲激响应的图形:

并且设计的阶数N=30。
4.1.2 最优等纹波线性相位法
这一类的滤波器设计方法可以表述为切比雪夫逼近问题。它将理想频率响应和实际频率响应之间的加权逼近误差均匀的分散到滤波器的整个阻带和通带,并且最小化误差。
通过matlab内部的工具箱得到的结果是：

设计的阶数N=15。
4。2 IIR类型的滤波
4。2。1 导数逼近法
导数逼近法的核心思想是用等效的差分方程逼近模拟滤波器表示的微分方程.通过关系式:

将s域的内容映射到z域当中去。

4.2。2 冲激不变法
我们的目的是得到这一函数。由变换关系可以知道，我们只需要得到相应的就好了。我们这里的思想就是用单位冲激响应的不变性，将模拟域的过渡到来。关系式为：。
切比雪夫滤波器I:

N=6;
得到的指标：

切比雪夫滤波器II：

N=6；
得到的指标：

4。2.3 双线性变换法
双线性变换是一个保形映射,它仅仅将轴变换到平面的单位圆一次，从而避免了频率成分的混叠。
在这里，从s平面到z平面的映射为：。
这就是双线性变换。
巴特沃斯滤波器：

N=9；
得到的指标：

椭圆滤波器：

N=5；
得到的指标：

5。滤波前后的特性对比
我们为了可以看到所有涉及的滤波器的效果,所以将全部的滤波器时频特性能够同时间进行对比， 具体的图形因为篇幅所限,所以我们不在本论文上一一列出，主要在课堂上演示.下面是凯撒窗滤波效果的图形(仅男声)：

6。对工作的概括
本次设计当中，我们以声音信号的滤波工作作为核心。主要的内容就是滤波器的选取过程.根据我们所学的知识，分别设计不同的FIR，IIR滤波器.设计过程中,我们不仅对滤波器的本身特性(包括四个主要的参数）进行了一些研究，还从滤波前后的时域、频域图形进行了对比。
本次课题的目的就是在针对具体实例下区分和选取最优的滤波器，并对设计的滤波器进行评价。我们在设计的过程中正是通过设计滤波器的阶数(对应占用的资源），设计后的实际(通阻带纹波、通阻带截止频率）进行考量，最终可以获得最优滤波效果的滤波器.
7。滤波器设计的结论
7.1 滤波器的设计过程
1) 确定指标
在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际运用中，数字滤波器常常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以两种方式给出。第一种是绝对指标。它提供对幅度响应函数的要求,一般应用于FIR滤波器的设计.第二种指标是相对指标。它以分贝值的形式给出要求。在实际工程中，这种指标最受欢迎。对于相位响应指标形式,通常希望系统在通频带中有线性相位。运用线性相位响应指标进行滤波设计具有如下优点：只包含实数算法，不涉及复数运算；不存在延迟失真（线性相位的含义就是不同频率的成分通过系统的时间是相同的)，只有固定数量的延迟；长度为N的滤波器(阶数为N—1)，计算量为N/2数量级。
2) 逼近
确定了技术指标之后，就可以建立一个目标的数字滤波器模型.通常采用理想的滤波器模型。之后，利用数字滤波器的设计方法，设计出一个世纪滤波器模型来逼近给定的目标。
3) 性能分析和计算机仿真
上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性,以验证设计结果是否满足指标要求;或者利用计算机仿真实现设计的滤波器，再分析滤波器结果来判断。
8。参考文献
1。 程佩青；数字信号处理； 2003
2. Sanjit K.Mitra；数字信号处理；2005
3. 百度文库
9。附录
程序代码：
fs=11000；
y=wavrecord(5＊fs,fs，’int16’);%设置采样频率
wavwrite(y，fs，’liming.wav')；%写入语音信号(要求在这一时刻开始录音）
［y1，fs］=wavread（'liming。wav’)；%读取该声音信号

fs=11000；
y=wavrecord(5＊fs，fs,’int16'）；％设置采样频率
wavwrite(y,fs,'shengqiuyi。wav'）；％写入语音信号（要求在这一时刻开始录音）
［y2，fs］=wavread（'shengqiuyi。wav’)；％读取该声音信号

%运行时的m文件
[y1,fs］=wavread(’liming。wav')；％读取该声音信号
［y2，fs]=wavread(’shengqiuyi.wav’）;％读取该声音信号
%-—-——---———-—-—-——-—-—---—-———--—---—---—-—---—-—-———-—-——--—---———-—-—--—
%—--—-----——--——————-—-—-----做原始的男声信号的时域和频域图形
%——---—-—--—-—--———--—-—-—--———---—-——-—---——--—-—-—--—--—-----————-——--——-
t=0：1/fs:（length(y1）-1)/fs； %确定时间轴的横轴坐标
plot（t，y1） %作原始语音信号的时域图形
%图像格式设置
title（'原始语音信号（男声)'）;
xlabel('时间（s)')；
ylabel（'幅度’）;

x=fft（y1。＊hamming(length（y1)));％作加窗傅里叶变换（汉明窗）
fm=5000＊length（x）/fs; %限定频率范围
f=(0:fm）＊fs/length(x）； ％确定频率刻度
plot（f,abs(x（1：length（f）)+eps)）;％画频谱图
%eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2.2204＊10^（-16)
％若某个量的绝对值小于eps，就认为这个量为0
title（’频谱图(男声)')；
xlabel（’频率(Hz）'）；
ylabel(’频谱幅度(db)'）；
%—-—-——---—————-——-—-—-———-—----—-——-———--—-————----—---————-——--———————--—
%——-----——-------—--—----——-—做原始的女声信号的时域和频域图形
%-——-——--——-————--——-———-—-—--———————--—--——-——--———-—---——---—-—-—-——-—---
t=0:1/fs:(length（y2)—1)/fs； %确定时间轴的横轴坐标
plot（t,y2) ％作原始语音信号的时域图形
％图像格式设置
title（’原始语音信号(女声）')；
xlabel（'时间(s)’);
ylabel(’幅度’);

x=fft（y2.＊hamming(length(y2))）；%作加窗傅里叶变换(汉明窗)
fm=5000＊length（x）/fs; ％限定频率范围
f=（0：fm）＊fs/length（x); ％确定频率刻度
plot（f,abs(x（1：length（f）)+eps））;%画频谱图
%eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2.2204*10^(—16）
%若某个量的绝对值小于eps，就认为这个量为0
title（'频谱图（女声）’)；
xlabel(’频率（Hz)’）;
ylabel（’频谱幅度(db)')；
%——--——-—--—-—-—---—-————--——-—————-—----———---—-—————----——————-—-——----—-
%—————---------——-——--——-——-—窗函数滤波器
％---——-----—-——-—--————-----—-—-—--——----—---------———----—---——-—-—-———---
％—-—-—----—-———————-——-—--—-——--—---——-—-—---————---—-——-
%—-——-——-——-—-———-—--—--——------—--————---—--——-—-—-————男
wavplay(y1)；
Hd = source1；
z1=filter（Hd，y1）;
wavplay(z1);
plot(t，z1） ％作原始语音信号的时域图形
％图像格式设置
title('滤波后的语音信号(男声)')；
xlabel（'时间（s)')；
ylabel（’幅度’）；
x=fft（z1。＊hamming（length（z1)））;%作加窗傅里叶变换(汉明窗）
fm=5000＊length(x）/fs; ％限定频率范围
f=（0:fm）＊fs/length（x）; %确定频率刻度
plot（f,abs(x(1：length（f)）+eps));%画频谱图
％eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2。2204＊10^（—16)
％若某个量的绝对值小于eps，就认为这个量为0
title('滤波后的频谱图（男声）’）;
xlabel（'频率（Hz)’)；
ylabel('频谱幅度（db）'）;
％-—---—-—--—-—-——--———-—---—--——-————-----———----————-——-
％-—--———-—-—-—-------—--——--—————--——---—---———-—---—-——女
wavplay（y2);
Hd = source1；
z11=filter（Hd，y2);
wavplay(z11);
plot(t,z11） ％作原始语音信号的时域图形
%图像格式设置
title('滤波后的语音信号（女声)’);
xlabel(’时间(s）');
ylabel('幅度’）;
x=fft(z11.＊hamming(length(z11）)）;％作加窗傅里叶变换（汉明窗）
fm=5000＊length(x）/fs； ％限定频率范围
f=(0：fm）*fs/length（x); ％确定频率刻度
plot（f,abs（x(1:length（f)）+eps））;%画频谱图
%eps 叫做机器的浮点运算误差限!
%PC机上eps的默认值为2.2204*10^（—16）
%若某个量的绝对值小于eps，就认为这个量为0
title('滤波后的频谱图（女声）')；
xlabel（’频率(Hz)’）;
ylabel(’频谱幅度（db）’);

％-—---—-—-—-———--——-—-——-—--——--———--———--—--———-——--—-——--——--———----—-——-
％——--——-—-——---—-———-----—--—等纹波滤波器
％----——----—-—--—--—-——--—-—-------—-————-——-----———----—-—--————---———----
wavplay(y1)；
Hd = source2；
z2=filter(Hd,y1）；
wavplay（z2);
plot(t，z2) ％作原始语音信号的时域图形
%图像格式设置
title('滤波后的语音信号(男声)'）;
xlabel(’时间(s）’);
ylabel('幅度’)；
x=fft(z2.＊hamming(length(z2）))；%作加窗傅里叶变换(汉明窗）
fm=5000*length(x）/fs； %限定频率范围
f=(0:fm）＊fs/length（x); ％确定频率刻度
plot(f，abs（x(1：length（f）)+eps）);%画频谱图
%eps 叫做机器的浮点运算误差限！
%PC机上eps的默认值为2.2204＊10^(—16)
％若某个量的绝对值小于eps，就认为这个量为0
title（’滤波后的频谱图(男声)’）;
xlabel('频率(Hz)’);
ylabel（’频谱幅度（db）’)；
%—--——-----—-——----—---——-———--———-—---—-—-——-————-—-——--
%---———-——-—--———---—-----———————-—-—-—-———-—-——-————--—女
wavplay（y2）；
Hd = source2；
z22=filter(Hd,y2);
wavplay（z22)；
plot（t,z22) %作原始语音信号的时域图形
%图像格式设置
title('滤波后的语音信号(女声)’）；
xlabel（'时间(s）'）;
ylabel（’幅度’)；
x=fft（z22。*hamming（length（z22）h）);%作加窗傅里叶变换(汉明窗）
fm=5000＊length（x)/fs; %限定频率范围
f=(0：fm）*fs/length(x); %确定频率刻度
plot(f,abs(x（1:length(f))+eps)）;％画频谱图
%eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2。2204*10^（—16)
%若某个量的绝对值小于eps，就认为这个量为0
title('滤波后的频谱图(女声）’）；
xlabel（'频率(Hz)’)；
ylabel（’频谱幅度(db）');

%—-—-——————-—---————---——-———--——----————-——--—---—-----—-—----————————-———
％—----——------—-——-—----—----巴特沃斯滤波器
%--————-——-——---——-——--—-——---————————---—-——-———-—-——--—--—-----——————--—-
wavplay（y1）；
Hd = source3;
z3=filter（Hd,y1)；
wavplay(z3）;
plot（t,z3） ％作原始语音信号的时域图形
%图像格式设置
title（'滤波后的语音信号(男声)'）；
xlabel（’时间（s）')；
ylabel('幅度’);
x=fft（z3.*hamming(length（z3))）；%作加窗傅里叶变换(汉明窗）
fm=5000＊length（x)/fs; ％限定频率范围
f=(0：fm）＊fs/length（x)； %确定频率刻度
plot(f，abs（x（1：length(f））+eps））;%画频谱图
％eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2.2204*10^(-16)
%若某个量的绝对值小于eps，就认为这个量为0
title（’滤波后的频谱图(男声）’);
xlabel（’频率(Hz)’)；
ylabel（’频谱幅度（db)’)；
％—-—-—-——--—-——---—--——-—-—-———-———--——--—---—---—-—---——
%-—-—————-—--——-——--———-—-—-——---—-------—————---——————-女
wavplay(y2）;
Hd = source3；
z33=filter(Hd,y2);
wavplay(z33）;
plot(t,z33） %作原始语音信号的时域图形
％图像格式设置
title（’滤波后的语音信号（女声)');
xlabel(’时间(s)');
ylabel（’幅度’);
x=fft（z33.＊hamming(length（z33)）);％作加窗傅里叶变换（汉明窗）
fm=5000＊length(x）/fs； %限定频率范围
f=(0:fm)＊fs/length（x); %确定频率刻度
plot（f，abs(x（1：length（f））+eps）);％画频谱图
％eps 叫做机器的浮点运算误差限！
%PC机上eps的默认值为2.2204*10^(-16)
％若某个量的绝对值小于eps，就认为这个量为0
title(’滤波后的频谱图(女声)’）；
xlabel（’频率（Hz）’）;
ylabel（’频谱幅度(db）');

%—--——-——--—--———-—-———--------——--—-——------—-—-------—-—-—----——---——----
％--————————-—---———-—-—------切比雪夫滤波器II
%--—-—--—-———----———----—-————-——-—--—-——--————-——-------——-------————--——-
wavplay(y1)；
Hd = source4；
z4=filter（Hd,y1);
wavplay(z4)；
plot（t，z4） ％作原始语音信号的时域图形
%图像格式设置
title（’滤波后的语音信号（男声)'）；
xlabel（'时间(s）’)；
ylabel(’幅度’）;
x=fft(z4。*hamming（length（z4）)）；％作加窗傅里叶变换(汉明窗）
fm=5000*length(x）/fs； ％限定频率范围
f=(0：fm)＊fs/length(x）; %确定频率刻度
plot（f，abs(x(1：length(f）)+eps）)；％画频谱图
％eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2.2204＊10^(—16)
%若某个量的绝对值小于eps，就认为这个量为0
title（'滤波后的频谱图（男声）’)；
xlabel(’频率(Hz)’)；
ylabel（'频谱幅度（db）’)；
％-—--——-—————-------——-————-——--————-———-———-----—-——-———
％-——--—--———-————-——--——--——--—-——-—-—----—-—-—-——-—--—-女
wavplay（y2）；
Hd = source4;
z44=filter（Hd,y2）；
wavplay(z44）;
plot(t，z44） %作原始语音信号的时域图形
%图像格式设置
title('滤波后的语音信号（女声）’);
xlabel('时间(s）')；
ylabel（'幅度’）;
x=fft（z44。＊hamming（length（z44))）;％作加窗傅里叶变换(汉明窗）
fm=5000*length(x）/fs； ％限定频率范围
f=(0:fm）*fs/length(x)； ％确定频率刻度
plot(f,abs(x(1:length(f))+eps)）；％画频谱图
%eps 叫做机器的浮点运算误差限！
%PC机上eps的默认值为2.2204＊10^(—16）
％若某个量的绝对值小于eps，就认为这个量为0
title('滤波后的频谱图(女声)'）；
xlabel（’频率（Hz）');
ylabel('频谱幅度(db)'）;

%—-——------———-—-————-——--——--——--——-—-——————---——-—-——-——--—---——-————--——
％-—--—-———--—--———--————-----椭圆滤波器
%--——-—-—-———————————————----—---—-——---——--—--—-—--—-————-—--—--—--———----
wavplay(y1)；
Hd = source5;
z5=filter(Hd,y1）；
wavplay(z5）；
plot（t,z5） ％作原始语音信号的时域图形
％图像格式设置
title（’滤波后的语音信号（男声）'）；
xlabel（’时间（s）')；
ylabel（'幅度’)；
x=fft（z5。＊hamming(length（z5）））;%作加窗傅里叶变换(汉明窗）
fm=5000＊length（x）/fs; %限定频率范围
f=(0：fm)＊fs/length（x)； ％确定频率刻度
plot(f，abs（x（1:length（f)）+eps））;%画频谱图
％eps 叫做机器的浮点运算误差限！
%PC机上eps的默认值为2。2204*10^（—16）
%若某个量的绝对值小于eps，就认为这个量为0
title（’滤波后的频谱图（男声)’);
xlabel（’频率（Hz）’)；
ylabel('频谱幅度（db)');
%---———----—-—-———-—-—-—-——--—-—-----————----——-—————--——
％-—---—----——-----------—--——-—-——----———-——-———--—---——女
wavplay(y2)；
Hd = source5;
z55=filter(Hd,y2)；
wavplay(z55);
plot(t，z55) %作原始语音信号的时域图形
％图像格式设置
title（’滤波后的语音信号（女声)’);
xlabel（’时间（s）')；
ylabel（'幅度’）;
x=fft（z55.＊hamming(length（z55)))；％作加窗傅里叶变换（汉明窗)
fm=5000*length(x)/fs； %限定频率范围
f=(0：fm)*fs/length(x）; ％确定频率刻度
plot（f,abs（x(1：length(f)）+eps)）;%画频谱图
％eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2.2204＊10^(—16）
％若某个量的绝对值小于eps，就认为这个量为0
title（'滤波后的频谱图（女声）')；
xlabel(’频率（Hz）')；
ylabel('频谱幅度（db）’)；

%—-—-——-—-——-----—---—--—--——--—-—-——---—————-—----————---—-——---—--—---—--
％----—---—-——-———--—-————--—-切比雪夫滤波器I
％————-------—-----—---————-—-——--———--—-—---—-—---—-—-——--—-—-—————-—-—-—-—
wavplay(y1)；
Hd = source6;
z6=filter（Hd，y1）;
wavplay(z6);
plot(t，z6） ％作原始语音信号的时域图形
%图像格式设置
title（’滤波后的语音信号（男声）');
xlabel（’时间（s)')；
ylabel（’幅度');
x=fft(z6.＊hamming(length(z6))）；%作加窗傅里叶变换（汉明窗)
fm=5000*length（x)/fs; ％限定频率范围
f=(0：fm）*fs/length(x）; ％确定频率刻度
plot(f,abs(x(1：length(f）)+eps)）；%画频谱图
%eps 叫做机器的浮点运算误差限！
％PC机上eps的默认值为2。2204＊10^（-16）
%若某个量的绝对值小于eps，就认为这个量为0
title(’滤波后的频谱图（男声）’)；
xlabel('频率（Hz)’);
ylabel（'频谱幅度(db）'）;
％————--———-——--—---——---—--—--—————-----—---—-———---—--——
％—-—-————--———--———-----—--———————--————--————-—-—-——---女
wavplay(y2);
Hd = source6；
z66=filter（Hd，y2)；
wavplay(z66)；
plot(t,z66) %作原始语音信号的时域图形
％图像格式设置
title('滤波后的语音信号（女声)')；
xlabel('时间（s)’)；
ylabel('幅度’）;
x=fft(z66.*hamming（length（z66）））;%作加窗傅里叶变换（汉明窗）
fm=5000＊length（x)/fs； %限定频率范围
f=（0:fm)＊fs/length（x）; %确定频率刻度
plot(f,abs(x（1：length(f)）+eps）);％画频谱图
％eps 叫做机器的浮点运算误差限!
%PC机上eps的默认值为2。2204＊10^（—16）
%若某个量的绝对值小于eps,就认为这个量为0
title（'滤波后的频谱图(女声)');
xlabel(’频率（Hz）'）；
ylabel（'频谱幅度（db)'）;
23