DSP课程设计语音信号的数字滤波处理.docx

DSP课程设计语音信号的数字滤波处理.docx

《DSP课程设计语音信号的数字滤波处理.docx》由会员分享，可在线阅读，更多相关《DSP课程设计语音信号的数字滤波处理.docx（34页珍藏版）》请在冰豆网上搜索。

DSP课程设计语音信号的数字滤波处理

摘要

数字信号处理（DSP）是20世纪60年代以来，随着信息学科和计算机学科的高速发展而迅速发展起来的一门新兴学科。

语音信号处理是研究用数字信号处理技术对语音信号处理的一门学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。

Matlab是一套用于科学工程计算的可视化高性能软件，是一种交互式的以矩阵为基本数据结构的系统，具有强大的矩阵运算能力，是本次课程设计的计算机辅助分析与设计部分的基础。

本设计通过运用MATLAB软件采集一段语音信号、回放语音信号并画出语音信号的时域波形和频谱图。

再在Matlab中分别设计不同形式的IIR和FIR数字滤波器。

之后对采集的语音信号经过不同的滤波器（低通、高通、带通）后，再对其进行频谱分析，通过分析滤波前后频谱图的区别，比较不同滤波器的滤波效果及其对语音信号的影响。

关键词：

DSP；语音信号处理；MATLAB；滤波器

1绪论

随着信息化的推进，数字信号处理的地位和作用变得越来越重要。

因为信息化的基础是数字化，而数字化的核心技术就是数字信号处理。

半个世纪以来，在如此强有力的需求牵引下，伴随着计算机技术、微电子技术日新月异的突破，数字信号处理的方法和应用越来越广泛和深入，发展十分迅速。

“数字信号处理（DSP）”是研究数字序列信号的表示方法，并对信号进行运算，以提取包含在其中的特殊信息的一门学科。

DSP是一门理论和实践密切结合的理论性和工程性都很强的学科，其理论性体现在，它综合应用数学、电路理论、信号与系统等领域的基础理论和方法，发展并形成了自己的理论体系，成为通信、雷达、声纳、电声、电视、测控、生物医学工程等众多学科和领域的重要理论基础与技术基础。

其工程性体现在，它的应用极为广泛，从科学技术的各个领域到国民经济建设的各个行业；从国防建设的各种武器装备到林林总总的消费类电子产品的设计与生产，都是数字信号处理技术的应用领域。

数字滤波器是数字信号处理的基础，用来对信号进行过滤、检测与参数估计等处理，在通信、图像、语音、雷达等许多领域都有着十分广泛的应用。

尤其在图像处理、数据压缩等方面取得了令人瞩目的进展和成就。

由此可见，数字滤波器的设计尤为重要。

目前数字滤波器的设计有许多现成的高级语言设计程序，但他们都存在设计效率较低，不具有可视图形，不便于修改参数等缺点，而Matlab为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。

他以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式的工作环境中。

尤其是Matlab工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。

其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。

本次课程设计主要介绍基于Matlab环境的IIR数字滤波器和FIR数字滤波器的设计，然后将设计结果应用于语音信号的处理，针对处理结果进一步分析了滤波器性能的。

2课程设计目的与内容

2.1设计目的

1.掌握数字信号处理的基本概念，基本理论和基本方法。

2.熟悉离散信号和系统的时域特性。

3.掌握序列快速傅里叶变换方法。

4.学会MATLAB的使用，掌握MATLAB的程序设计方法。

5.掌握利用MATLAB对语音信号进行频谱分析。

6.掌握滤波器的网络结构。

2.2设计内容

本课程设计的内容分为两个部分：

预习题：

设计卷积运算的演示程序

是编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应。

设计题：

用切比雪夫模拟滤波器（低通、带通、高通）对有噪声信号的语音进行滤波

用hanning窗（低通、带通、高通）对有噪声信号的语音进行滤波。

3设计原理

3.1用窗函数法设计FIR滤波器

根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。

根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n），如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：

在一般情况下，hd（n）是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。

用窗函数wd（n）将hd（n）截断，并进行加权处理，得到

如果要求线性相位特性，则h（n）还必须满足：

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。

要根据所设计的滤波特性正确选择其中一类。

例如，要设计线性相位低通特性可选择h（n）=h（N-1-n）一类，而不能选h（n）=-h（N-1-n）一类。

验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。

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

脉冲响应不变法的主要缺点是产生频率响应的混叠失真。

这是因为从S平面到Ｚ平面是多值的映射关系所造成的。

为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。

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

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

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

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

（3-1）

式中,T仍是采样间隔。

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

将式（3-1）写成

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

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

z=es1T

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

（3-2）

（3-3）

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

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

首先,把z=ejω，可得

（3-4）

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

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

因此

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

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

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

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

双线性变换法优缺点如下：

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

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

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

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

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

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

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

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

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

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

的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3-3所示。

图3-3双线性变换法幅度和相位特性的非线性映射

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

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

4设计程序的调试和运行结果

4.1预习题部分

分别动态演示两个序列进行线性卷积x1（n）﹡x2（n）和圆周卷积x1（n）⊙x2 （n）的过程；要求分别动态演示翻转、移位、乘积、求和的过程，程序及波形图如下所示；

翻转程序：

x1=[200957050202];

y=fliplr（x1）

subplot（2,1,1）;stem（x1）;

xlabel（'n'）;ylabel（'x1[n]'）;

subplot（2,1,2）;stem（y）;

xlabel（'n'）;ylabel（'y[n]'）;

翻转程序波形图如图4-1所示：

图4-1翻转波形图

卷积程序如下：

x1=[200957050202];

x2=[1,2.43,6.17,12.93,22.17,32.25,40.88,45.87,45.87,40.88,32.25,22.17,12.93,6.17,2.43,1.0000]

y=conv（x1,x2）

stem（y）;xlabel（'n'）;ylabel（'y[n]'）;

卷积图形如图4-2所示：

图4-2卷积波形图

圆周卷积程序如下：

m1=[200957050202];

m2=[12.436.1712.9322.1732.2540.8845.8745.8740.8832.2522.1712.936.172.431];

subplot（711）;stem（m1）;ylabel（'x1（n）'）;

subplot（712）;stem（m2）;ylabel（'x2（n）'）;

m2=fliplr（m2）;subplot（713）;stem（m2）;ylabel（'shift'）;

k=length（m1）;j=length（m2）;N=max（j,k）;M=N;

x2=[zeros（1,M-j-1）m20];x1=[m1zeros（1,M-k）];y=ones（1,N）;

whileN>0

x2=[x2（M）x2（1:

M-1）];

y（N）=sum（x1.*x2）;

ifN>M-3

subplot（7,1,4+M-N）;stem（x2）;ylabel（'shift'）;

end

N=N-1;

end

subplot（717）;stem（y）;ylabel（'y'）;

图形如图4-3所示：

图4-3圆周卷积图

采样程序如下所示：

t=0:

0.0005:

0.3;

xa=8*exp（-8*sqrt

（2）*pi*t）.*sin（8*sqrt

（2）*pi*t）;

subplot（611）;plot（t,xa）;ylabel（'xa'）;

t=0.015;

n=0:

t:

0.3;

pt=ones（0.3/t+1）;

subplot（612）;stem（n,pt）;ylabel（'samplesignal'）;

xaa=8*exp（-8*sqrt

（2）*pi*n）.*sin（8*sqrt

（2）*pi*n）;

subplot（613）;stem（n,xaa）;ylabel（'sampledsignal'）;

y1=fft（xa,4096）;

subplot（614）;plot（abs（y1））;ylabel（'samplesignal'）;

y2=fft（xaa,4096）;

subplot（615）;plot（abs（y2））;ylabel（'samplesignal'）;

t=0.005;

n=0:

t:

0.3;

xaa=8*exp（-8*sqrt

（2）*pi*n）.*sin（8*sqrt

（2）*pi*n）;

y2=fft（xaa,4096）;

subplot（616）;plot（abs（y2））;ylabel（'samplesignal'）;

波形图如图4-4所示：

图4-4采样图

卷积翻转采样求和程序如下所示：

m1=wavread（'c:

\windows\media\notify.wav',10）;m1=m1';m1=m1（1,:

）

m2=[12.436.1712.9322.1732.2540.8845.8745.8740.8832.2522.1712.936.172.431];

subplot（711）;stem（m1）;ylabel（'x1（n）'）;

subplot（712）;stem（m2）;ylabel（'x2（n）'）;

y=conv（m1,m2）;

subplot（713）;stem（y）;ylabel（'y'）;

N=max（length（m1）,length（m2））;

X1=fft（m1,25）;X2=fft（m2,25）;

subplot（714）;stem（X1）;ylabel（'X1'）;

subplot（715）;stem（X2）;ylabel（'X2'）;

Y1=X1.*X2;

subplot（716）;stem（Y1）;ylabel（'Y1'）;

Y2=fft（y）;

subplot（717）;stem（Y2）;ylabel（'Y2'）;

波形图如图4-5所示：

图4-5卷积翻转采样求和程序图

4.2设计题部分

4.2.1语音信号的频谱分析

[y,fs,bits]=wavread（'e:

\over.wav',[116384]）;

sound（y,fs,bits）;

Y=fft（y,1024）;

f=fs*（0:

511）/1024;

figure

（1）;

subplot（211）;plot（abs（Y（1:

512）））;

plot（y）;title（'原始语音信号'）;

subplot（212）;plot（f,abs（Y（1:

512）））;title（'原始语音信号频谱'）;

figure

（2）;

freqz（y）;title（'频率响应图'）;

图4-6图4-7

4.2.2污染信号的频谱分析

[y,fs,bits]=wavread（'e:

\over.wav',[1 16384]）;

f=fs*（0:

511）/1024;

t=0:

1/16384:

（size（y）-1）/16384; %将所加噪声信号的点数调整到与原始信号相同

Au=0.03;

d=[Au*cos（2*pi*10000*t）]';   %噪声为5kHz的余弦信号 要用随机噪声

y2=y+d;

sound（y2,16384）;    %播放加噪声后的语音信号

Y2=fft（y2,1024）;

figure

（1）

subplot（2,1,1）;plot（t,y）;

subplot（2,1,2）;plot（t,y2）;

title（'加噪后的信号'）;

xlabel（'时间 n'）;

ylabel（'幅值'）;

figure

（2）

subplot（2,1,1）;

plot（f,abs（Y（1:

512）））;

title（'原始语音信号频谱'）;

xlabel（'Hz'）;

ylabel（'幅值'）;

subplot（2,1,2）;

plot（f,abs（Y2（1:

512）））;

title（'加噪后的信号频谱'）;

xlabel（'Hz'）;

ylabel（'幅值'）;

figure（3）; 

plot（d）;title（'所添加噪声'）;axis（[0 500 -0.04 0.04]）;

图4-8图4-9

图4-10

4.2.3IIR滤波器——切比雪夫型滤波器

1）切比雪夫Ⅱ型—低通滤波器

图4—11图4—12

图4—13图4—14

图4—15

2）切比雪夫Ⅱ型—高通滤波器

图4—16图4—17

图4—18图4—19

图4—20

3）切比雪夫Ⅱ型—带通滤波器

图4—21图4—22

图4—23图4—24

图4—25

4.2.4FIR滤波器—hanning窗滤波器

1）低通

图4—26图4—27

图4—28图4—29

图4—30图4—31

2）高通

图4—32图4—33

图4—34图4—35

图4—36图4—37

3）带通

图4—38图4—39

图4—40图4—41

图4—42图4—43

5课程设计的思考与体会

5.1设计思考

1．双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？

从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？

答：

在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。

如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。

2．能否利用公式完成脉冲响应不变法的数字滤波器设计？

为什么？

答：

IIR数字滤波器的设计实际上是求解滤波器的系数

和

，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。

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

但是它的缺点是，存在频率混迭效应，故只适用于阻带的模拟滤波器。

5.2设计体会

为期两周的课程设计结束了，在老师和同学的帮助下，我顺利完成了我的语音信号的数字滤波处理的课程设计。

MATLAB和DSP都是我们上学习所学课程，好几个月的时间没有接触，开学的第一个任务就是课程设计。

当我再一次看到这些曾经熟悉的知识时，只觉得措手不及。

通过小组同学的商量，我们决定先翻阅以前的课本以及实验指导书，再到网上以及图书馆查阅资料。

复习MATLAB以及DSP的知识，了解语音信号处理这门技术，从而使自己的设计过程更加的得心应手。

在这个学习的过程当中，我们分工明确，相互鼓励，以尽可能高的效率完成了本次的课程设计。

在设计的过程当中，我们也遇到了不少的问题，比如说语音信号的采集，我们从来没有接触过，通过查阅资料，我们学会了不少的新的MATLAB函数，也让我意识到MATLAB的功能相当强大，是一个非常好的学习工具。

另外，滤波器的种类繁多，很参数指标也容易搞混淆，在老师的帮助下，我们也攻克了滤波器这一难关。

在课程设计之前，我自认为自己的DSP和MATLAB学得不错，可是到了真正要应用的时候才知道自己学得很不扎实。

通过这次课程设计，让我懂得了在以后的学习当中，不能仅仅注重于学习书本上的知识，还要懂得将知识运用于实践。

参考文献

[1]丁玉美.数字信号处理（第二版）[M].西安：

西安电子科技大学出版社，2001．

[2]王创新,文卉.数字信号处理试验指导书[Z].长沙：

长沙理工大学印刷（内部使用），2006.

[3]陈怀琛.数字信号处理及其MATLAB实现[M].北京：

电子工业出版社，2001.

[4]陈怀琛.MATLAB及在电子信息课程中的应用[M].北京：

电子工业出版社，2002.

[5]A.V.奥本海姆，R.W.谢弗.数字信号处理[M].北京：

科学出版社,2001.

[6]胡广书.数字信号处理——理论、算法与实现（第二版）[M].北京：

电子工业出版社.

致谢

在本次课程设计中，我常常会遇到自己弄不明白的问题，比如说程序看不懂，调试时显示结果与预期不一致等等。

而每当我遇到困难的时候，老师总会耐心地给我讲解，帮助我发现问题，纠正错误；同学也会很热情的帮助我修改和完善设计方案。

在这样一个良好的环境下，我顺利地完成了本次数字信号处理课程设计，让我发现了自己所学知识的盲点和漏洞，提升了我的动手能力。

特在此对给予我帮助的老师和同学们表示衷心的感谢！

你们辛苦了！

附录MATLAB程序代码

1.加入噪声的语音信号的调用及分析程序

[y,fs,bits]=wavread（'e:

\over.wav',[116384]）;

f=fs*（0:

511）/1024;

t=0:

1/16384:

（size（y）-1）/16384;%将所加噪声信号的点数调整到与原始信号相同

Au=0.03;d=[Au*cos（2*pi*10000*t）]';%噪声为5kHz的余弦信号要用随机噪声

y2=y+d;sound（y2,16384）;%播放加噪声后的语音信号

Y2=fft（y2,1024）;

figure

（1）;

subplot（2,1,1）;plot（t,y）;

subplot（2,1,2）;plot（t,y2）;title（'加噪后的信号'）;

xlabel（'时间n'）;ylabel（'幅值'）;

figure

（2）;

subplot（2,1,1）;plot（f,abs（Y（1:

512）））;

title（'原始语音信号频谱'）;xlabel（'Hz'）;ylabel（'幅值'）;

subplot（2,1,2）;plot（f,abs（Y2（1:

512）））;

title（'加噪后的信号频谱'）;xlabel（'Hz'）;ylabel（'幅值'）;

figure（3）;

plot（d）;title（'所添加噪声'）;axis（[0500-0.040.04]）;

2.hanning窗带通滤波器程序

[y,fs,bits]=wavread（'e:

\over.wav',[116384]）;

f=fs*（0:

7500）/15000;

%d=0.1*rand（1,length（y））';%随机信号

t=0:

1/16384:

（size（y）-1）/16384;%将所加噪声信号的点数调整到与原始信号相同

Au=0.03;d=[Au*cos（2*pi*5000*t）]';

y2=y+d;sound（y2,16384）;%播放加噪声后的语音信号

Y=fft（y,15000）;

f=fs*（0:

7500）/15000;

Y2=fft（y2,15000）;

figure

（1）

subplot（2,1,1）;plot（y）;title（'原始语音信号'）;

subplot（2,1,2）;plot（y2）;

title（'加噪后的信号'）;xlabel（'时间n'）;ylabel（'幅值'）;

figure

（2）

subplot（2,1