【数字信号处理课程设计】基于MATLAB实现语音信号的采集与处理（偏重滤波）

设计数字滤波器和画出其频率响应给出各滤波器的性能指标；给定滤波器的性能指标如下：（1）低通滤波器的性能指标：fb=1000Hz,fc=1200Hz,As=100dB,Ap=1dB.（2）高通滤波器的性能指标：fc=4800Hz,fb=5000Hz,As=100dB,Ap=1dB.（3）带通滤波器的性能指标： fb1=1200Hz, fb2=3000Hz,fc1=1000Hz, fc2=3200H

舞果sight

20120人浏览 · 2022-07-20 20:04:20

舞果sight · 2022-07-20 20:04:20 发布

这个项目在我的B站上有专门的视频演示：

【数字信号处理课程设计】基于MATLAB实现语音信号的采集与处理（侧重滤波）+【通信原理课程设计】基于MATLAB实现的信源编码系统仿真 +【GUI设计】_哔哩哔哩_bilibili

一、目标与任务

1、语音信号的采集。利用 Windows 下的录音机，录制一段自己的话音，时间在 1s 内，然后在 Matlab 软件平台下，利用函数 wavread 对语音信号进行采样，记住采样频率和采样点数。

2、语音信号的频谱分析。在 Matlab 中，可以利用函数 fft 对信号进行快速傅立叶变换，得到信号的频谱特性，要求学生首先画出语音信号的时域波形，然后对语音信号进行频谱分析。

3、设计数字滤波器和画出其频率响应给出各滤波器的性能指标；给定滤波器的性能指标如下：

（1）低通滤波器的性能指标：fb=1000Hz,fc=1200Hz,As=100dB,Ap=1dB.

（2）高通滤波器的性能指标：fc=4800Hz,fb=5000Hz,As=100dB,Ap=1dB.

（3）带通滤波器的性能指标： fb1=1200Hz, fb2=3000Hz,fc1=1000Hz, fc2=3200Hz,As=100dB,Ap=1dB. 采用窗函数法和双线性变换法设计上面要求的 3 种滤波器，并画出滤波器的频率响应；

4、用滤波器对信号进行滤波，然后用自己设计的滤波器对采集到的信号进行滤波，画出滤波后信号的时域波形及频谱，并对滤波前后的信号进行对比，分析信号的变化；

5、回放语音信号，分析滤波前后的语音变化;设计系统界面为了使编制的程序操作方便，设计处理系统的用户界面，在所设计的系统界面上可以实现上述要求中的包括采集、分析、滤波等全部内容，并能够选择滤波器的类型，输入滤波器的参数、显示滤波器的频率响应等。

二、原理介绍

2.1 录音原理

Matlab 中有两个函数可以用于实现录音功能，一个是waverecord 函数调用 windows 音频设备进行录音，主要用于实时语音采集。第二个是使用Matlab 中的 wavread 函数也可以进行语音采集。不过 wavread 函数只能读取计算机中已经有的“.wav” 文件，也就是说必须事先使用 windows 的录音机等设备把有关的语音资料录制并存储为“.wav” 文件。显然 wavread 函数不能用于实时语音信号采集。另外语音采集的频率和精度不能改变。回放时也有两个函数可以选用。 Sound 函数和 wavplay，只不过后者的功能没有前者强大。不过在本设计中主要是对已有的信号进行处理，所以选用wavread 录放。同时，由于信号的时域描述只能反映信号的幅值随时间的变化情况，除只有一个频率分量的简谐波外般很难明确揭示信号的频率组成和各频率分量的大小，因此，我们要把它转换成频域来分析，它能够提供比时域信号波形更直观，丰富的信息。

2.2 滤波器的设计原理及设计方法

FIR 滤波器的设计问题在于寻求一系统函数H(z) ，使其频率响应H(ejω) 逼近滤波器要求的理想频率响应Hd(ejω) 其对应的单位脉冲响应hd(n) 。

1、用窗函数设计FIR 滤波器的基本方法

设计思想：从时域从发，设计h(n) 逼近理想hd(n) 。设理想滤波器Hd(ejω) 的单位脉冲响应为hd(n) 。以低通线性相位 FIR 数字滤波器为例。

hd(n) 一般是无限长的，且是非因果的，不能直接作为FIR 滤波器的单位脉冲响应。要想得到一个因果的有限长的滤波器 h(n) ，最直接的方法是截断h(n)=hd(n)w(n) ，即截取为有限长因果序列，并用合适的窗函数进行加权作为FIR 滤波器的单位脉冲响应。按照线性相位滤波器的要求，h(n) 必须是偶对称的。对称中心必须等于滤波器的延时常数，即

用矩形窗设计的FIR 低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的 9%，这个现象称为吉布斯（Gibbs）效应。为了消除吉布斯效应，一般采用其他类型的窗函数。

2、典型的窗函数

（1）矩形窗(Rectangle Window)

其频率响应和幅度响应分别为：

（2）三角形窗(Bartlett Window)

其频率响应为：

（3）汉宁(Hanning)窗，又称升余弦窗

其频率响应和幅度响应分别为：

（4）汉明(Hamming)窗，又称改进的升余弦窗

其幅度响应为：

（5）布莱克曼(Blankman)窗，又称二阶升余弦窗

其幅度响应为：

（6）凯泽(Kaiser)窗

其中：β是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，β 越大，过渡带越宽，阻带越小衰减也越大。是第一类修正零阶贝塞尔函数。

若阻带最小衰减表示为As=-20log10⁡δ ，β 的确定可采用下述经验公式：

若滤波器通带和阻带波纹相等即δp=δs 时，滤波器节数可通过

下式确定：

式中：

3．利用窗函数设计FIR 滤波器的具体步骤如下：

（1）按允许的过渡带宽度△ω 及阻带衰减 AS 选择合适的窗函数，并估计节数N：其中A 由窗函数的类型决定。

（2）由给定的滤波器的幅频响应参数求出理想的单位脉冲响应hd(n) 。

（3）确定延时值。

（4）计算滤波器的单位取样响应h(n) ，

（5）验算技术指标是否满足要求。

2.3 IIR 数字滤波器设计原理

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

假设 M≤ N，当M＞ N时, 系统函数可以看作一个 IIR 的子系统和一个 (M-N)的FIR子系统的级联。 IIR 数字滤波器的设计实际上是求解滤波器的系数 ak 和 bk ，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在 S平面上去逼近，就得到模拟滤波器；如果在 z平面上去逼近，就得到数字滤波器。

设计通数字滤波器通常可以归纳为如图所示的两种常用方法。

图 2-1 数字带通滤波器设计的两种方法

方法 1: 首先设计一个模拟原型低通滤波器，然后通过频率变换成所需要的模拟高通滤波器，最后再使用冲激不变法或双线性变换成相应的数字高通滤波器。

方法 2: 先设计一个模拟原型低通滤波器，然后采用冲激响应不变法或双线性变换法将它转换成数字原型低通滤波器，最后通过频率变换把数字原型低通滤波器变换成所需要的数字高通滤波器。

本课程设计采用第一种设计方法，先构造一个巴特沃斯模拟低通滤波器，然后将模拟低通滤波器转换成模拟高通滤波器，最后利用双线性变换将模拟高通滤波器转换成数字高通滤波器。

2.4 双线性变换法

为了克服冲激响应法可能产生的频率响应的混叠失真，这是因为从 S 平面到Ｚ平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到 - π / T～π / T 之间，再用 z=esT 转换到 Z 平面上。也就是说，第一步先将整个 S平面压缩映射到 S1 平面的 - π / T～π/ T 一条横带里；第二步再通过标准变换关系 z=esT 将此横带变换到整个 Z 平面上去。这样就使 S平面与 Z 平面建立了一一对应的单值关系，消除了多值变换性，

也就消除了频谱混叠现象，映射关系如图 2-2

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

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

式中 , T 仍是采样间隔。

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

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

再将 S1平面通过以下标准变换关系映射到 Z 平面z=esT 从而得到 S 平面和 Z 平面的单值映射关系为：

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

首先 , 把z=eiw 可得

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

其次，将 s=σ+jw 代入，得

因此

由此看出，当 <0 时， | z|<1 ；当 >0 时， | z|>1 。也就是说， S 平面的左半平面映射到 Z 平面的单位圆内， S平面的右半平面映射到 Z 平面的单位圆外，S平面的虚轴映射到 Z 平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。

双线性变换法优缺点：双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为 S平面与 Z平面是单值的一一对应关系。S 平面整个 j Ω 轴单值地对应于 Z 平面单位圆一周，即频率轴是单值变换关系。

这个关系重写如下：

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

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

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

三、GUI界面设计与实现

页面设计如图3-1所示：

图 3-1 GUI界面

四、基于MATLAB仿真

4.1实验过程

（1）加入噪声。

向原始信号中加入正弦干扰噪声或高斯噪声，加入正弦干扰噪声后，如图4-1所示。

图 4-1 加入正弦干扰噪声前后

加入高斯干扰噪声后，如图4-2所示。

图 4-2 加入高斯噪声前后

（2）滤波器的频谱特性

IIR线性滤波器

巴特沃斯滤波器低通幅频相频特性如图4-3所示

图 4-3 巴特沃斯低通滤波器幅频相频特性

车比雪夫滤波器高通幅频相频特性如图4-4所示

图 4-4 车比雪夫高通滤波器幅频相频特性

椭圆滤波器带通幅频相频特性如图4-5所示

图 4-5 椭圆型带通滤波器幅频相频特性

FIR窗函数滤波器

Blackman窗低通抽样响应和频谱特性如图4-6所示

图 4-6 Blackman窗低通抽样响应和频谱

Kaiser窗高通抽样响应和频谱特性如图4-7所示

图 4-7 kaiser窗低通抽样响应和频谱

（3）滤波前后波形

低通滤波后，如图4-8所示。

图 4-8 低通滤波

高通滤波后，如图4-9所示。

图 4-9 高通滤波

带通滤波后，如图4-10所示。

图 4-10 带通滤波

4.2 结果分析

1.通过本次设计发现，当抽样率Fs大于15000Hz时，滤波器的设计除低通外，高通和带通都会出现异常，主要异常为矩阵奇异值，或者NAN、INF；而Fs小于10000Hz时，低通和带通设计正常，但高通设计会报错异常，而Fs取11025Hz时，可以实现设计需求。

2.使用巴特沃斯滤波器，低通、高通、带通滤波效果都不能实现，而切比雪夫1型滤波器和椭圆滤波器除高通外，低通和带通效果都很明显。

3.布莱克曼窗和凯塞窗在阻带衰落为100dB时设计出的滤波器性能几乎一致，且除了高通外，低通和带通对加噪后的音频滤波效果都很明显。

4.加正弦噪声时，由于本设计所加单频主要位于100Hz处，所以通过修改适当的滤波器参数指标，信噪比较大或较小时，高通和带通都能很好的滤除噪声。

5.加高斯噪声时，无论采用哪种滤波器都不能完美的滤除噪声，但由于高斯噪声的随机特性，会发现加高斯噪声后频谱会出现很多尖峰，所以采用低通滤波器，在信噪比较小，即噪声干扰很大的时候，滤波效果会很明显。

五、总结

5.1 函数用法总结

1.[n, Wnj= buttord(Wp, Ws,Rp,Rs)

其中Wp和ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

2.chebyl函数的用法为: [b,a]-chebyl1(n,Rp, Wn,/ftype/)

在使用chebyl函数设计IR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn， cheblord函数可在给定滤波器性能的情况下，选择契比雪夫1型滤波器的最小阶和截止频率Wn.

3.cheblord函数的用法为：[n, Wnj-cheblord(Wp, Ws,Rp,Rs)

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

4.butter函数的用法为：[b,a]=butter(n,Wn,/ftype/)

其中n代表滤波器阶数， Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn.

5，Butterworth模拟低通滤波器函数语句： [z，p.kj—buttap（n）

其中n为滤波器的阶数， z. p、k分别为零点、极点和增益

6. [bz, az]=impinvar (b, atC, FS t, Fp)

式中， b， a为模拟滤波器分子和分母多项式系数向量： Fs为采样频率（所滤波数据），单位Hz，缺省时为1Hz.Fp为预畸变频率（Prewarped frequency），是一个“匹配”频率，在该频率上，频率响应在变换前后和模拟频率可精确匹配

在MATLAB中，函数bilinear采用双线性变换法实现模拟s域至数字z域的映射，直接用于模拟滤波器变换为数字滤波器。其调用方式为：

[2zd, pd, kd]-bil inear (a, p, k,Fs)

[numd, dend]-bilinear (num, den, Fs)

式中， z，p分别为模拟滤波器零点、极点列向量： k为模拟滤波器的增益：Fs为采样频率，单位Hz. zd， pd， kd为数字滤波器的零极点和增益.num， den分别为模拟滤波器传递函数分子和分母多项式系数向量，模拟滤波器传递函数具有下面的形式：

H(s)-mum(5) mum(1)s+ +mum (mn)s+ num (nn+1)

den (s) den()s"..+den(nd )s + den (nd +1)

nund和dend分别为数字滤波器传递函数分子和分母多项式系数向量

7，函数trez用于求数字逃波器的频率响应，其调用格式为：on[h, w]-]fregz (b, a, nL,'whole']) ;或[h, f]=freqz (b, a, nt,' whole'],Fs); 式中， b，a为数字滤波器分子和分母多项式的系数， n为复数频率的响应点数，为整数，最好为2的幂，缺省时为512： Fs为采样频率，单位Hz，如果给定该值，则f位置输出为频率Hz，若没有给定，则按角频率（Angular frequency）给定f的频率矢量：， whole'表示返回的频率或值包含2平面整个单位圆频率矢量，即0—27：缺省时，频率f或值包含z平面上半单位圆（0—x）之间等间距n个点频率矢量。h为复频率响应： w为n点频率向量（单位rad）； f为n点频率向量（Hz），函数返回值缺省时，绘制幅频响应和相频响应图。

5.2 心得体会

通过这次课设，对设计数字滤波器的整个设计过程有了很好的掌握，懂得了设计滤波器的基本方法。对IIR滤波器和FIR滤波器的应用场景和区别有了一定了解，熟练了对巴特沃斯、车比雪夫、椭圆形滤波器和几种常见窗函数的使用。同时也熟悉了 MATLAB的环境，巩固了相关知识。初步掌握了 MATLAB语言在数字信号处理中一些基本库函数的调用和编写基本程序等应用；熟悉了滤波器设计的一般原理，对滤波器有了一个感性的认识；学会了数字高通滤波器设计的一般步骤；加深了对滤波器设计中产生误差的原因以及双线性变换法优缺点的理解和认识。总之，使理论联系了实际，巩固并深化了对课本基本知识的认识和理解，使理论得以升华。

在做本次课程设计的过程中，我和队友曾经遇到过问题，在解决问题的过程中锻炼了自己处理问题的能力，彼此相互沟通，通过查阅书籍、搜索资料等方法使问题得以解决。培养了我们的耐心，也深深感受到自身知识的有限。通过使所学应用于实践，提高了动手能力。