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

1、 语音信号的数字滤波处理摘要数字信号处理（DSP）是 20 世纪 60 年代以来，随着信息学科和计算机学科的高速发展而迅速发展起来的一门新兴学科。语音信号处理是研究用数字信号处理技术对语音信号处理的一门学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。Matlab 是一套用于科学工程计算的可视化高性能软件，是一种交互式的以矩阵为基本数据结构的系统，具有强大的矩阵运算能力，是本次课程设计的计算机辅助分析与设计部分的基础。本设计通过运用 MATLAB 软件采集一段语音信号、回放语音信号并画出语音信号的时域波形和频谱图。再在 Matlab 中分别设计不同形式的 IIR 和 FIR 数字滤波器

2、。之后对采集的语音信号经过不同的滤波器（低通、高通、带通）后，再对其进行频谱分析，通过分析滤波前后频谱图的区别，比较不同滤波器的滤波效果及其对语音信号的影响。关键词：关键词： DSP；语音信号处理；MATLAB；滤波器 语音信号的数字滤波处理目录1 绪论.12课程设计目的与内容 .22.1 设计目的 .22.2 设计内容 .23设计原理 .23.1用窗函数法设计 FIR 滤波器.23.2用双线性变换法设计 IIR 数字滤波器.34 设计程序的调试和运行结果.74.1 预习题部分 .74.2 设计题部分 .114.2.1 语音信号的频谱分析 .114.2.2 污染信号的频谱分析 .124.2.3

3、 IIR 滤波器切比雪夫型滤波器.144.2.4 FIR 滤波器hanning 窗滤波器 .175 课程设计的思考与体会.205.1 设计思考 .205.2 设计体会 .20参考文献.22致谢.23附录 MATLAB 程序代码 .24 语音信号的数字滤波处理第 1 页 共 36 页1 绪论随着信息化的推进，数字信号处理的地位和作用变得越来越重要。因为信息化的基础是数字化，而数字化的核心技术就是数字信号处理。半个世纪以来，在如此强有力的需求牵引下，伴随着计算机技术、微电子技术日新月异的突破，数字信号处理的方法和应用越来越广泛和深入，发展十分迅速。 “数字信号处理(DSP)”是研究数字序列信号的表

4、示方法，并对信号进行运算，以提取包含在其中的特殊信息的一门学科。DSP 是一门理论和实践密切结合的理论性和工程性都很强的学科，其理论性体现在，它综合应用数学、电路理论、信号与系统等领域的基础理论和方法，发展并形成了自己的理论体系，成为通信、雷达、声纳、电声、电视、测控、生物医学工程等众多学科和领域的重要理论基础与技术基础。其工程性体现在，它的应用极为广泛，从科学技术的各个领域到国民经济建设的各个行业；从国防建设的各种武器装备到林林总总的消费类电子产品的设计与生产，都是数字信号处理技术的应用领域。数字滤波器是数字信号处理的基础，用来对信号进行过滤、检测与参数估计等处理，在通信、图像、语音、雷达等

5、许多领域都有着十分广泛的应用。尤其在图像处理、数据压缩等方面取得了令人瞩目的进展和成就。由此可见，数字滤波器的设计尤为重要。目前数字滤波器的设计有许多现成的高级语言设计程序，但他们都存在设计效率较低，不具有可视图形，不便于修改参数等缺点，而 Matlab 为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。他以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式的工作环境中。尤其是 Matlab 工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。本次课程设计主要介绍基于 Matla

6、b 环境的 IIR 数字滤波器和 FIR 数字滤波器的设计，然后将设计结果应用于语音信号的处理，针对处理结果进一步分析了滤波器性能的。 语音信号的数字滤波处理第 2 页 共 36 页2课程设计目的与内容2.1 设计目的1. 掌握数字信号处理的基本概念，基本理论和基本方法。2. 熟悉离散信号和系统的时域特性。3. 掌握序列快速傅里叶变换方法。4. 学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法。5. 掌握利用 MATLAB 对语音信号进行频谱分析。6. 掌握滤波器的网络结构。2.2 设计内容本课程设计的内容分为两个部分：预习题：设计卷积运算的演示程序是编写程序演示采样定理（时域1

7、2采样、频谱周期延拓） ，同时演示采样频率小于 2fc 时，产生的混叠效应。设计题：用切比雪夫模拟滤波器（低通、带通、高通）对有噪声信号的1语音进行滤波用 hanning 窗（低通、带通、高通）对有噪声信号的语音进行2滤波。3设计原理3.1 用窗函数法设计 FIR 滤波器根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度 N（或阶数 M=N-1） ，窗函数类型可根据最小阻带衰减 As 独立选择，因为窗口长度 N对最小阻带衰减 As 没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度 N，设待求滤波器的过渡带宽为 w，它与窗口长度 N 近似成反比，窗函

8、数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的 N，在 N 和窗函数类型确定后，即可调用 MATLAB 中的窗函数求出窗函数 wd（n） 。 语音信号的数字滤波处理第 3 页 共 36 页根据待求滤波器的理想频率响应求出理想单位脉冲响应 hd(n)，如果给出待求滤波器频率应为 Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从 w=0 到 w=2 采样 N 点，采用离散傅里叶反变换（IDFT）即可求出。用窗函数 wd(

9、n)将 hd(n)截断，并进行加权处理，得到如果要求线性相位特性， 则 h(n)还必须满足：根据上式中的正、 负号和长度 N 的奇偶性又将线性相位 FIR 滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择 h(n)=h(N-1-n)一类，而不能选 h(n)=-h(N-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用 freqz 子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。3.2 用双线性变换法设计 IIR 数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从

10、S 平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T 之间，再用z=esT 转换到 Z 平面上。也就是说，第一步先将整个 S 平面压缩映射到 S1 平面的-/T/T 一条横带里；第二步再通过标准变换关系 z=es1T 将此横带变换到整个 Z 平面上去。这样就使 S 平面与 Z 平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图 3-1 所示。1( )()2jj nddh nHeed( )( ) ( )dh nhnn( )(1)h nh Nn 语音信号的数字滤波处理第 4 页 共 36 页

11、图 3-1 双线性变换的映射关系为了将 S 平面的整个虚轴 j 压缩到 S1 平面 j1 轴上的-/T 到 /T 段上，可以通过以下的正切变换实现 （3-1）式中,T 仍是采样间隔。当 1 由-/T 经过 0 变化到 /T 时， 由-经过 0 变化到+，也即映射了整个 j 轴。将式（3-1）写成将此关系解析延拓到整个 S 平面和 S1 平面，令 j=s，j1=s1，则得再将 S1 平面通过以下标准变换关系映射到 Z 平面z=es1T从而得到 S 平面和 Z 平面的单值映射关系为： (3-2) (3-3) 语音信号的数字滤波处理第 5 页 共 36 页式（3-2）与式（3-3）是 S 平面与 Z

12、 平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（3-1）与式（3-2）的双线性变换符合映射变换应满足的两点要求。首先,把 z=ej，可得 (3-4)即 S 平面的虚轴映射到 Z 平面的单位圆。其次，将 s=+j 代入式（3-4） ，得因此由此看出，当 0 时，|z|0 时，|z|1。也就是说，S 平面的左半平面映射到 Z 平面的单位圆内，S 平面的右半平面映射到 Z 平面的单位圆外，S 平面的虚轴映射到 Z 平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点如下：双线性变换法与脉冲响应不变法相比，其主要的优点是避

13、免了频率响应的混叠现象。这是因为 S 平面与 Z 平面是单值的一一对应关系。S 平面整个 j轴单值地对应于 Z 平面单位圆一周，即频率轴是单值变换关系。这个关系如式（3-4）所示，重写如下：上式表明，S 平面上 与 Z 平面的 成非线性的正切关系，如图 2-2 所示。 语音信号的数字滤波处理第 6 页 共 36 页由图 2-2 看出，在零频率附近，模拟角频率 与数字频率 之间的变换关系接近于线性关系；但当 进一步增加时， 增长得越来越慢，最后当时， 终止在折叠频率 = 处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。图 3-2 双线性变换法

14、的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（3-4）及图 3-2 所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性） ，不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图 3-3 所示。图 3-3 双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，

15、仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通 语音信号的数字滤波处理第 7 页 共 36 页过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。4 设计程序的调试和运行结果4.1 预习题部分分别动态演示两个序列进行线性卷积 x1(n)x2(n)和圆周卷积 x1(n)x2 (n)的过程；要求分别动态演示翻转、移位、乘积、求和的过程，程序及波形图如下所示；翻转程序：x1=2 0 0 9 5 7 0 5 0 2 0 2;y=fliplr(x1)subplot(2,1,1);stem(x1);xlabe

16、l(n);ylabel(x1n);subplot(2,1,2);stem(y);xlabel(n);ylabel(yn);翻转程序波形图如图 4-1 所示： 语音信号的数字滤波处理第 8 页 共 36 页0246810120510nx1n0246810120510nyn图 4-1 翻转波形图卷积程序如下：x1=2 0 0 9 5 7 0 5 0 2 0 2;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.0000y=conv(x1,x2)stem(y

17、);xlabel(n);ylabel(yn);卷积图形如图 4-2 所示： 语音信号的数字滤波处理第 9 页 共 36 页051015202530020040060080010001200nyn图 4-2 卷积波形图圆周卷积程序如下：m1=2 0 0 9 5 7 0 5 0 2 0 2;m2=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;subplot(711);stem(m1);ylabel(x1(n);subplot(712);stem(m2);ylabel(x2(n)

18、;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) m2 0;x1=m1 zeros(1,M-k);y=ones(1,N);while N0 x2=x2(M) x2(1:M-1); y(N)=sum(x1.*x2); if NM-3 subplot(7,1,4+M-N);stem(x2);ylabel(shift); end N=N-1; 语音信号的数字滤波处理第 10 页 共 36 页endsubplot(717);stem(

19、y);ylabel(y);图形如图 4-3 所示：0246810120510 x1(n)0246810121416050 x2(n)0246810121416050shift0246810121416050shift0246810121416050shift0246810121416050shift0246810121416010002000y图 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=

20、0:t:0.3;pt=ones(0.3/t+1);subplot(612);stem(n,pt);ylabel(sample signal);xaa=8*exp(-8*sqrt(2)*pi*n).*sin(8*sqrt(2)*pi*n);subplot(613);stem(n,xaa);ylabel(sampled signal);y1=fft(xa,4096);subplot(614);plot(abs(y1);ylabel(sample signal);y2=fft(xaa,4096);subplot(615);plot(abs(y2);ylabel(sample signal); 语音信

21、号的数字滤波处理第 11 页 共 36 页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(sample signal);波形图如图 4-4 所示：00.050.10.150.20.250.30.35-505xa00.050.10.150.20.250.30.3500.51sample signal00.050.10.150.20.250.30.35-505sampled signal050010001500200

22、0250030003500400045000200400sample signal0500100015002000250030003500400045000510sample signal05001000150020002500300035004000450002040sample signal图 4-4 采样图卷积翻转采样求和程序如下所示：m1=wavread(c:windowsmedianotify.wav,10);m1=m1;m1=m1(1,:)m2=1 2.43 6.17 12.93 22.17 32.25 40.88 45.87 45.87 40.88 32.25 22.17 12.

23、93 6.17 2.43 1;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); 语音信号的数字滤波处理第 12 页 共 36 页subplot(715);stem(X2);ylabel(X2);Y1=X1.*X2;subplot(716);

24、stem(Y1);ylabel(Y1);Y2=fft(y);subplot(717);stem(Y2);ylabel(Y2);波形图如图 4-5 所示：12345678910-101x1(n)0246810121416050 x2(n)0510152025-101y0510152025-101X10510152025-5000500X20510152025-101Y10510152025-101Y2图 4-5 卷积翻转采样求和程序图4.2 设计题部分4.2.1 语音信号的频谱分析y,fs,bits=wavread(e:over.wav,1 16384);sound(y,fs,bits);Y=f

25、ft(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(原始语音信号频谱); 语音信号的数字滤波处理第 13 页 共 36 页figure(2);freqz(y);title(频率响应图); 图 4-6 图 4-74.2.2 污染信号的频谱分析y,fs,bits=wavread(e:over.wav,1 16384);f=fs*(0:511)/1024;t=0:1/16384:(size(y

26、)-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(原始语音信号频谱);xla

27、bel(Hz); 语音信号的数字滤波处理第 14 页 共 36 页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-104.2.3 IIR 滤波器切比雪夫型滤波器1）切比雪夫型低通滤波器 语音信号的数字滤波处理第 15 页 共 36 页图 411 图 412图 413 图 414图 4152）切比雪夫型高通滤波器 语音信号的数字滤波处理

28、第 16 页 共 36 页图 416 图 417图 418 图 419图 4203）切比雪夫型带通滤波器 语音信号的数字滤波处理第 17 页 共 36 页图 421 图 422图 423 图 424图 4254.2.4 FIR 滤波器hanning 窗滤波器1)低通 语音信号的数字滤波处理第 18 页 共 36 页图 426 图 427图 428 图 429图 430 图 4312)高通 语音信号的数字滤波处理第 19 页 共 36 页图 432 图 433图 434 图 435图 436 图 4373)带通 语音信号的数字滤波处理第 20 页 共 36 页图 438 图 439图 440 图

29、 441图 442 图 4435 课程设计的思考与体会 语音信号的数字滤波处理第 21 页 共 36 页5.1 设计思考1双线性变换法中 和 之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？ 答：在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。如以上实验过程中，采用双线性变化法设计的 butter 和 cheby1 数字滤波器，从图中可以看到这种非线性关系。2能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？答：IIR 数字滤波器的设计

30、实际上是求解滤波器的系数 和 ，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在 S 平面上去逼近，就得到模拟滤波器；如果在 z 平面上去逼近，就得到数字滤波器。但是它的缺点是，存在频率混迭效应，故只适用于阻带的模拟滤波器。5.2 设计体会为期两周的课程设计结束了，在老师和同学的帮助下，我顺利完成了我的语音信号的数字滤波处理的课程设计。MATLAB 和 DSP 都是我们上学习所学课程，好几个月的时间没有接触，开学的第一个任务就是课程设计。当我再一次看到这些曾经熟悉的知识时，只觉得措手不及。通过小组同学的商量，我们决定先翻阅以前的课本以及实验指导书，再

31、到网上以及图书馆查阅资料。复习 MATLAB 以及 DSP 的知识，了解语音信号处理这门技术，从而使自己的设计过程更加的得心应手。在这个学习的过程当中，我们分工明确，相互鼓励，以尽可能高的效率完成了本次的课程设计。在设计的过程当中，我们也遇到了不少的问题，比如说语音信号的采集，我们从来没有接触过，通过查阅资料，我们学会了不少的新的 MATLAB 函数，也让我意识到 MATLAB 的功能相当强大，是一个非常好的学习工具。另外，滤波器的种类繁多，很参数指标也容易搞混淆，在老师的帮助下，我们也攻克了滤波器这一难关。在课程设计之前，我自认为自己的 DSP 和 MATLAB 学得不错，可是到了真正 语音

32、信号的数字滤波处理第 22 页 共 36 页要应用的时候才知道自己学得很不扎实。通过这次课程设计，让我懂得了在以后的学习当中，不能仅仅注重于学习书本上的知识，还要懂得将知识运用于实践。参考文献1 丁玉美.数字信号处理（第二版）M.西安：西安电子科技大学出版社，20012 王创新,文卉.数字信号处理试验指导书Z.长沙：长沙理工大学印刷（内部使用） ，2006.3 陈怀琛.数字信号处理及其 MATLAB 实现M.北京：电子工业出版社，2001.4 陈怀琛.MATLAB 及在电子信息课程中的应用M.北京：电子工业出版社，2002. 语音信号的数字滤波处理第 23 页 共 36 页5 A.V.奥本海姆

33、，R.W.谢弗.数字信号处理M.北京：科学出版社,2001.6 胡广书.数字信号处理理论、算法与实现（第二版）M.北京：电子工业出版社.致谢在本次课程设计中，我常常会遇到自己弄不明白的问题，比如说程序看不懂，调试时显示结果与预期不一致等等。而每当我遇到困难的时候，老师总会耐心地给我讲解，帮助我发现问题，纠正错误；同学也会很热情的帮助我修改和完善设计方案。在这样一个良好的环境下，我顺利地完成了本次数字信号处理课程设计，让我发现了自己所学知识的盲点和漏洞，提升了我的动手能力。特在此对给予我帮助的老师和同学们表示衷心的感谢！你们辛苦了！ 语音信号的数字滤波处理第 24 页 共 36 页附录 MATL

34、AB 程序代码1.加入噪声的语音信号的调用及分析程序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); 语音信号的数字滤波处理第 25 页 共 36 页figure(1);subplot(2,1,1);plot(t,y);

35、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);2. hanning 窗带通滤波器程序y,fs,

36、bits=wavread(e:over.wav,1 16384);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(原始语音信号);s

37、ubplot(2,1,2);plot(y2);title(加噪后的信号);xlabel(时间 n);ylabel(幅值);figure(2)subplot(2,1,1);plot(abs(Y(1:7500);title(原始语音信号频谱);xlabel(Hz);ylabel(幅值); 语音信号的数字滤波处理第 26 页 共 36 页subplot(2,1,2);plot(abs(Y2(1:7500);title(加噪后的信号频谱);axis(0 8000 0 100);xlabel(Hz);ylabel(幅值); wlp=2000*2*pi/fs;wup=4000*2*pi/fs;Rp=1;R

38、s=100;wls=1500*2*pi/fs;wus=7000*2*pi/fs;Bt=wlp-wls; %计算过渡带宽度NO=ceil(6.2*pi/Bt); %汉宁窗计算所需 h(n)长度 NO %ceil(x)取大于等于 x 的最小整数N=NO+mod(NO+1,2); %确保 h(n) 长度 N 为奇数wc=(wlp+wls)/2/pi,(wup+wus)/2/pi; %计算理想高通滤波器通带截止频率（归一化）hn=fir1(N-1,wc,hanning(N);h,w=freqz(hn,1,15000);figure(3)freqz(hn,1,15000);xlabel(频率(HZ);y

39、label(幅度(db);title(带通汉宁窗数字滤波器幅频特性);x=fftfilt(hn,y2);X=fft(x,15000);figure(4);subplot(2,1,1);plot(abs(Y2(1:7500);%axis(0,12000,0,30);title(滤波前信号频谱);subplot(2,1,2);plot(abs(X(1:7500);%axis(0,12000,0,30);title(滤波后信号频谱); figure(5);subplot(2,1,1);plot(y2);title(滤波前信号波形);subplot(2,1,2);plot(x);%axis(0,1,-

40、1,1);title(滤波后信号波形);axis(0,18000,-0.2,0.2);sound(x,fs,bits);3. hanninghanning 窗低通滤波器程序 语音信号的数字滤波处理第 27 页 共 36 页y,fs,bits=wavread(e:over.wav,1 16384);f=fs*(0:7500)/15000;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,1500

41、0);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,1);plot(abs(Y(1:7500);title(原始语音信号频谱);xlabel(Hz);ylabel(幅值);subplot(2,1,2);plot(abs(Y2(1:7500);title(加噪后的信号频谱);axis(0 8000 0 100)

42、;xlabel(Hz);ylabel(幅值); wp=5500*2*pi/fs;Rp=1;Rs=100;ws=6000*2*pi/fs;Bt=ws-wp; %计算过渡带宽度NO=ceil(6.2*pi/Bt); %汉宁窗计算所需 h(n)长度 NO %ceil(x)取大于等于 x 的最小整数N=NO+mod(NO+1,2); %确保 h(n) 长度 N 为奇数wc=(wp+ws)/2/pi; %计算理想高通滤波器通带截止频率（归一化）hn=fir1(N-1,wc,hanning(N);h,w=freqz(hn,1,15000);figure(3)freqz(hn,1,15000);xlabel

43、(频率(HZ);ylabel(幅度(db); 语音信号的数字滤波处理第 28 页 共 36 页title(低通汉宁窗数字滤波器幅频特性); x=fftfilt(hn,y2);X=fft(x,15000);figure(4);subplot(2,1,1);plot(abs(Y2(1:7500);%axis(0,12000,0,30);title(滤波前信号频谱);subplot(2,1,2);plot(abs(X(1:7500);%axis(0,12000,0,30);title(滤波后信号频谱); figure(5);subplot(2,1,1);plot(y2);title(滤波前信号波形)

44、;subplot(2,1,2);plot(x);%axis(0,1,-1,1);title(滤波后信号波形);axis(0,18000,-0.2,0.2);sound(x,fs,bits);4. hanning 窗高通滤波器程序y,fs,bits=wavread(e:over.wav,1 16384);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,1

45、6384); %播放加噪声后的语音信号 Y=fft(y,15000);f=fs*(0:7500)/15000;Y2=fft(y2,15000); figure(1)subplot(2,1,1);plot(y);title(原始语音信号); 语音信号的数字滤波处理第 29 页 共 36 页subplot(2,1,2);plot(y2);title(加噪后的信号);xlabel(时间 n);ylabel(幅值); figure(2)subplot(2,1,1);plot(abs(Y(1:7500);title(原始语音信号频谱);xlabel(Hz);ylabel(幅值);subplot(2,1,

46、2);plot(abs(Y2(1:7500);title(加噪后的信号频谱);axis(0 8000 0 100);xlabel(Hz);ylabel(幅值); wp=6000*2*pi/fs;Rp=1;Rs=100;ws=5500*2*pi/fs;Bt=wp-ws; %计算过渡带宽度NO=ceil(6.2*pi/Bt); %汉宁窗计算所需 h(n)长度 NO %ceil(x)取大于等于 x 的最小整数N=NO+mod(NO+1,2);wc=(wp+ws)/2/pi; %计算理想高通滤波器通带截止频率（归一化）hn=fir1(N-1,wc,high,hanning(N);h,w=freqz(h

47、n,1,15000); figure(3)freqz(hn,1,15000);xlabel(频率(HZ);ylabel(幅度(db);title(高通汉宁窗数字滤波器幅频特性); x=fftfilt(hn,y2);X=fft(x,15000);figure(4);subplot(2,1,1);plot(abs(Y2(1:7500); title(滤波前信号频谱);subplot(2,1,2);plot(abs(X(1:7500);title(滤波后信号频谱); figure(5);subplot(2,1,1);plot(y2); 语音信号的数字滤波处理第 30 页 共 36 页title(滤波

48、前信号波形);subplot(2,1,2);plot(x);title(滤波后信号波形);axis(0,18000,-0.2,0.2);sound(x,fs,bits);5. 切比雪夫型带通滤波器y,fs,bits=wavread(e:over.wav,1 16384);f=fs*(0:7500)/15000;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:750

49、0)/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,1);plot(abs(Y(1:7500);title(原始语音信号频谱);xlabel(Hz);ylabel(幅值);subplot(2,1,2);plot(abs(Y2(1:7500);title(加噪后的信号频谱);axis(0 8000 0 100);xlabel(Hz);yl

50、abel(幅值); %滤波器设计fss=15000wp=600*2*pi/fss,4000*2*pi/fss;Rp=1;Rs=100;ws=300*2*pi/fss,7000*2*pi/fss;wo=sqrt(600*2*pi/fss)*(4000*2*pi/fss);Bw=4000*2*pi/fss-600*2*pi/fss; 语音信号的数字滤波处理第 31 页 共 36 页n,wn=cheb2ord(wp,ws,1,100,s);z,p,k=cheb2ap(n,99);b0,a0=zp2tf(z,p,k);bt,at=lp2bp(b0,a0,wo,Bw);bz,az=bilinear(bt

51、,at,1);h,w=freqz(bz,az);f=w*fss/(2*pi); figure(3)plot(f,20*log10(abs(h);xlabel(频率(HZ);ylabel(幅度(db);title(带通切比雪夫数字滤波器幅频特性);axis(0,12000,-150,10); x=filter(bz,az,y2);X=fft(x,15000);figure(4);subplot(2,1,1);plot(abs(Y2(1:7500); title(滤波前信号频谱);subplot(2,1,2);plot(abs(X(1:7500); title(滤波后信号频谱); figure(5

52、);subplot(2,1,1);plot(y2);title(滤波前信号波形);subplot(2,1,2);plot(x); title(滤波后信号波形);axis(0,18000,-0.2,0.2);sound(x,fs,bits); 6. 切比雪夫型低通滤波器y,fs,bits=wavread(e:over.wav,1 16384);f=fs*(0:7500)/15000;t=0:1/16384:(size(y)-1)/16384; %将所加噪声信号的点数调整到与原始信号相同Au=0.03;d=Au*cos(2*pi*5000*t);y2=y+d;sound(y2,16384); %播

53、放加噪声后的语音信号 Y=fft(y,15000);f=fs*(0:7500)/15000;Y2=fft(y2,15000); 语音信号的数字滤波处理第 32 页 共 36 页figure(1)subplot(2,1,1);plot(y);title(原始语音信号);subplot(2,1,2);plot(y2);title(加噪后的信号);xlabel(时间 n);ylabel(幅值); figure(2)subplot(2,1,1);plot(abs(Y(1:7500);title(原始语音信号频谱);xlabel(Hz);ylabel(幅值);subplot(2,1,2);plot(abs(Y2(1:7500);title(加噪后的信号频谱);a