语音信号的数字滤你波处理.doc

《语音信号的数字滤你波处理.doc》由会员分享，可在线阅读，更多相关《语音信号的数字滤你波处理.doc（45页珍藏版）》请在三一办公上搜索。

1、 语音信号的数字滤波处理（3）目录一 绪论.21.1设计目的.21.2 设计内容.21.3 设计题目.3二 设计指标要求.3三 预习题.43.1线性卷积.43.2FFT快速卷积.63.3圆周卷积.83.4采样定理演示.11四 课程设计原理.134.1巴特沃斯低通滤波器 .134.2用窗函数法设计FIR滤波器.14五 设计的步骤和过程.155.1设计方案.155.2设计步骤.16六 用巴特沃斯设计IIR滤波器.306.1巴特沃斯低通滤波器.196.2巴特沃斯带通滤波器.206.3巴特沃斯高通滤波器.23七 用hamming窗设计FIR滤波器.257.1blackman高通滤波器.25 7.2bl

2、ackman低通滤波器.267.3blackman带通滤波器.28八 心得体会.30参考文献. 31附录A巴特沃斯低通滤波器程序. .32附录B巴特沃斯高通滤波器程序. 34附录C巴特沃斯带通滤波器程序. 37附录D：blackman高通滤波器程序.40附录E：blackman低通滤波器程序 42附录F：blackman带通滤波器程序 44 一 绪论 在本次数字信号处理实习过程中，需要我们用自己设计的数字滤波器分别对被不同噪声污染的信号进行滤波。 数字滤波器可分为IIR和FIR两大类。而且数字滤波是数字信号处理的重要内容，对于IIR数字滤波器的设计，我们需要借助模拟原型滤波器，再将模拟滤波器转

3、化为数字滤波器，文中采用的设计方法是脉冲响应不变法、双向性变换法和完全函数设计法；对于FIR数字滤波器的设计，可以根据所给定的频率特性直接设计，文中采用的设计方法是窗函数法。根据IIR滤波器和FIR滤波器的特点，在MATLAB坏境下分别用双线性变换法设计IIR和用窗函数设计FIR数字滤波器，并对采集的语音信号进行分析，最后给出了IIR和FIR对语音滤波的效果。 数字滤波器可分为IIR和FIR两大类。而且数字滤波是数字信号处理的重要内容，对于IIR数字滤波器的设计，我们需要借助模拟原型滤波器，再将模拟滤波器转化为数字滤波器，文中采用的设计方法是脉冲响应不变法、双向性变换法和完全函数设计法；对于F

4、IR数字滤波器的设计，可以根据所给定的频率特性直接设计，文中采用的设计方法是窗函数法。根据IIR滤波器和FIR滤波器的特点，在MATLAB坏境下分别用双线性变换法设计IIR和用窗函数设计FIR数字滤波器，并对采集的语音信号进行分析，最后给出了IIR和FIR对语音滤波的效果。 Digital filter can be divided into two categories: IIR and FIR. And the digital filter is an important part of digital signal processing, for the design of IIR di

5、gital filter, we need the help of analog prototype filter, and then the analog filter transformation digital filter and the design method is pulse response method, two-way transformation method and full function design method; for the design of FIR digital filter can according to the given the frequ

6、ency characteristics of the direct design, the design method is the window function method. According to the characteristics of IIR filter and FIR filter, in the bad environment of the MATLAB respectively with bilinear transform design method of IIR and using window function design FIR digital filte

7、r and on the acquisition of the speech signal analysis, finally gives the IIR and FIR filtering of the speech.1.1设计目的此次课程设计的目的是为了让我们综合运用数字信号处理和MATLAB去设计语音信号去噪的数字滤波器，使学生加强自身通过实践发现问题，探讨问题从而解决问题的能力。不但能够使学生对这门课程有更加深入的了解与运用而且能够让学生学会自主设计的具体步骤和方法，提高自我创新的能力，提高实际应用水平。1.2 设计内容（1）设计卷积运算的演示程序可输入任意两个序列x1(n)、x2(n

8、)，指定x1(n)为自己的学号, x1(n)=2,0,1,15,7,0,5,0,2,1,2例如x2(n)= 1, 2.43, 6.17,12.93,22.17,32.25,40.88, 45.87, 45.87, 40.88, 132.25, 22.17, 12.93, 6.17, 2.43,1.0000。分别动态演示两个序列进行线性卷积x1(n)x2(n)和圆周卷积x1(n)x2(n)的过程；要求分别动态演示翻转、移位、乘积、求和的过程；圆周卷积默认使用两个序列中的最大长度，但卷积前可以指定卷积长度N用以进行混叠分析；改变圆周卷积长度N，根据实验结果分析2类卷积的关系。 在计算机操作系统下选

9、一段声音文件(XP系统在“C:WINDOWSMedia”)，读取文件取10ms的声音数据产生时域序列x1(n),序列内容自定义。利用x2(n)= 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。利用FFT实现快速卷积，验证时域卷积定理，并与直接卷积进行效率对比（验证时采用matlab子函数）。（2）编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应对下面连续信号进行采样：，A为幅度因子，a为衰减因子，为模拟角频率，

10、其中n为学号(例如，王墨同学n=23)，要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示不同采样频率下对应信号的频谱。1.3设计题目利用Windows下的录音机或其他软件，进行语音信号的采集（*.wav）；语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；污染信号的频谱分析，画出被污染的语音信号时域波形和频谱； 根据有关的频谱特性，采用间接法设计IIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个IIR滤波器）a. 模拟滤波器类型：巴特沃思滤波器（低通、带通

11、、高通）b总体要求：Matlab原程序+仿真波形+技术指标 根据有关的频谱特性，采用直接法设计FIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个FIR滤波器）a. 滤波器类型：hamming窗（低通、带通、高通）b总体要求：Matlab原程序+仿真波形+技术指标+窗函数 用自己设计的这些滤波器分别对被不同噪声污染的信号进行滤波用自己设计的这些滤波器分别对被不同噪声污染的信号进行滤波；分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号。二 设计指标要求1. 掌握数字信号处理的基本概念，基本理论和基本方法。2. 熟悉离散信号和系统的

12、时域特性。3. 掌握序列快速傅里叶变换方法。4. 学会MATLAB的使用，掌握MATLAB的程序设计方法。5. 掌握利用MATLAB对语音信号进行频谱分析。6. 掌握滤波器的网络结构。三 预习题3.1 线性卷积x1(n)=2,0,1,2,5,7,0,5,0,2,2,5,x2(n)= 2,4,5,6,2,7,4,8,9,6,4,8,3,5,7,8 ，其线性卷积动态演示Matlab程序如下：clf;n1=1:40;n2=1:40;x1=zeros(1,12),2,0,1,2,5,7,0,5,0,2,2,5,zeros(1,16); x2=zeros(1,12),2,4,5,6,2,7,4,8,9,

13、6,4,8,3,5,7,8,zeros(1,12);x11=zeros(size(x1);subplot(3,1,1)stem(n1,x1)subplot(3,1,2)stem(n2 ,x2)pause(2);for i=1:24 x11(i)=x1(25-i);endfor i=25:40 x11(i)=0;endsubplot(3,1,1)stem(n1,x11);pause(2)x3=zeros(size(x11);for i=1:28x11=zeros(1,1),x11(1:(40-1);x3(i)=sum(x11.*x2);subplot(3,1,1)stem(n1,x11)subp

14、lot(3,1,3)x4=zeros(1,12),x3(1:i),zeros(1,28-i);stem(n1,x4)pause(2)end 运行结果：图3.1 线性卷积3.2 FFT快速卷积x2(n)= 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，其FFT快速卷积Matlab程序如下：clf;x1=wavread(c:WindowsMediading.wav,10);x1=x1;x1=x1(1,:);x2=1 2.43 6.17 12.93

15、22.17 32.25 40.88 45.87 45.87 40.88 32.25 22.17 12.93 6.17 2.43 1;subplot(7,1,1);stem(x1);ylabel(x1(n);title(x1(n);subplot(7,1,2);stem(x2);ylabel(x2(n);title(x2(n);y=conv(x1,x2);subplot(7,1,3);stem(y);ylabel(y);title(x1(n)与x2(n)的卷积);N1=length(x1);N2=length(x2);N=N1+N2-1;X1=fft(x1,N);X2=fft(x2,N);sub

16、plot(7,1,4);stem(X1);ylabel(X1);title(x1(n)的N点DFT);subplot(7,1,5);stem(X2);ylabel(X2);title(x2(n)的N点DFT);Y1=X1.*X2;subplot(7,1,6);stem(Y1);ylabel(Y1);title(X1与X2相乘的结果);Y2=ifft(Y1);subplot(7,1,7);stem(Y2);ylabel(Y2);title(Y1的IDFT结果);运行结果： 图3.2 FFT快速卷积3.3圆周卷积x1(n)=2,0,1,2,5,7,0,5,0,2,2,5，x2(n)= 2,4,5,

17、6,2,7,4,8,9,6,4,8,3,5,7,8，其圆周卷积Matlab程序如下：clf;x=zeros(1,12),2,0,1,2,5,7,0,5,0,2,2,5,zeros(1,16);y=zeros(1,12),2,4,5,6,2,7,4,8,9,6,4,8,3,5,7,8,zeros(1,12);subplot(4,1,1)stem(x);subplot(4,1,2)stem(y);N1=length(x);N2=length(y);N=N1;if N2N1 N=N2;endM=input(请输入圆周卷积长度n);if M=N N=M;endif N=N1 x=x;else x=x,

18、zeros(1,N-N1);endX=zeros(N,N);x1=zeros(1,N);pause(2)for n=1:N if n=1 x1(n)=x(n); else x1(n)=x(N-n+2); endendX(1,1:N)=x1;for i=2:Npause(0.3) for n=1:N if n=1 X(i,n)=X(i-1,N); else X(i,n)=X(i-1,n-1); endendsubplot(4,1,3);stem(X(i,1:N);endif N2Ny=y,zeros(1,N-N2);endY=y;W=X*Y;W=W;n=1:N;subplot(4,1,4);st

19、em(n,W);当圆周卷积长度为40，50时，结果如图3.3、3.4： 图3.3 圆周卷积（N=40）图3.4 圆周卷积（N=50）3.4采样定理演示连续信号：，A为幅度因子，a为衰减因子，为模拟角频率，其中n为25，其采样Matlab程序如下：fs=input(请输入采样频率:);n=input(请输入两位学号:);t=0:1/fs:1;t1=0:0.0005:1;A=n;a=sqrt(2)*n*pi;w=a;x=A*exp(-a*t).*sin(w*t);x1=A*exp(-a*t1).*sin(w*t1);subplot(5,1,1)plot(x1);xlabel(时间(s);ylabe

20、l(幅度);xn=ones(1,length(t);subplot(512)stem(xn);xlabel(时间(s);ylabel(幅度);yn=xn.*x;subplot(5,1,3)stem(yn);xlabel(时间(s);ylabel(幅度);y=fft(yn,fs);f=fs*(0:fs/2)/fs;subplot(5,1,4)plot(f,y(1:fs/2+1);xlabel(频率(HZ);ylabel(幅度);yn=ifft(y);subplot(5,1,5)plot(yn);xlabel(时间(s);ylabel(幅度); 当fs=80 时运行结果如图3.5：图3.5 采样定

21、理演示（fs=80） 当fs= 40时运行结果如图3.6：图3.6 采样定理演示（fs=40）四 课程设计原理4.1 巴特沃斯低通滤波器 (1)特性:巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。 一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率 单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形

22、状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。 (2)传递函数：巴特沃斯低通滤波器可用如下振幅的平方对频率的公式表示：其中,= 滤波器的阶数= 截止频率 = 振幅下降为 -3分贝时的 频率= 通频带边缘频率在通频带边缘的数值4.2 用窗函数法设计FIR滤波器 根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的

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

24、权处理，得到 如果要求线性相位特性， 则h(n)还必须满足： 根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。五 设计方案及步骤5.1 设计方案 1.要求利用windows下的录音机（开始程序附件娱乐录音机，文件属性立即转换8000Hz，8位，单声道）录制一段自己

25、的话音，或者采用Windows自带的声音文件（默认为22050Hz），时间控制在几秒左右。然后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解采样频率、采样位数等概念。 wavread函数调用格式： y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。 y,fs,nbits=wavread(file)，采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。 y=wavread(file,N)，读取前N点的采样值放在向量y中。 y=wavread(file,N1,

26、N2)，读取从N1点到N2点的采样值放在向量y中。 2语音信号的频谱分析 要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。 3设计数字滤波器和画出频率响应 根据语音信号的特点给出有关滤波器的性能指标，例如：1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=2800 Hz，fp=3000 Hz，As=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，f

27、c1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB。 4.回放语音信号 在Matlab中,函数sound可以对声音进行回放。其调用格式：sound(x,fs,bits);可以感觉滤波前后的声音有变化。 5.方案设计流程图如图3.5.1所示。（图为巴特沃斯滤波器，blackman窗滤波器与巴特沃斯类似）图5.1 方案设计流程图5.2设计步骤（1）语音信号的采集Matlab程序如下:运行结果如图5.2x,FS,bits=wavread(C:WINDOWSMediair_inter.wav);x=x(:,1);figure(1);subplot(2,1,1);plot(x)

28、;sound(x,FS,bits); title(语音信号时域波形图)y=fft(x,3260);f=(FS/1630)*1:1630; subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);title(语音信号频谱图);运行结果如图5.2所示：图5.2 语音信号（2）对语音信号加噪Matlab程序如下：t=0:length(x)-1;zs=0.002*cos(2*pi*10000*t/22050);zs0=0.002*cos(2*pi*10000*t/22050000);figure(2);subplot(2,1,1)plot(zs0)title(噪声信号波

29、形);zs1=fft(zs,1200);sound(zs,FS,bits); subplot(2,1,2)plot(f(1:600),abs(zs1(1:600);title(噪声信号频谱);x1=x+zs;sound(x1,FS,bits); y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title(加入噪声后的信号波形);subplot(2,1,2);plot(f(1:600),abs(y1(1:600);title(加入噪声后的信号频谱);运行结果如图5.3所示：（a）（b）图5.3 对语音信号加噪六 用巴特沃斯设计IIR滤波器6.1

30、巴特沃斯低通滤波器Matlab程序见附录A巴特沃斯低通滤波器程序，运行程序波形图分析如下：图6.1.1图6.1.2图6.1.3图6.1.4图6.1.56.2. 巴特沃斯带通滤波器Matlab程序见附录B巴特沃斯带通滤波器程序，运行程序波形图分析如下：图6.2.1图6.2.2图6.2.3图6.2.4图6.2.56.3 巴特沃斯高通滤波器Matlab程序见附录C巴特沃斯高通滤波器程序，运行程序波形图分析如下：图6.3.1图6.3.2图6.3.3图6.3.4图6.3.5七 blackman窗设计FIR滤波器7.1 blackman高通滤波器调试程序见附录Dblackman窗高通滤波器程序,运行程序波

31、形图分析：图7.1.1图7.1.2图7.1.3图7.1.47.2 blackman低通滤波器调试程序见附录blackman窗低通滤波器程序，运行程序波形图分析如下：图7.2.1图7.2.2图7.2.3图7.2.47.3 blackman带通滤波器调试程序见附录blackman窗带通滤波器程序，运行程序波形图分析如下：图7.3.1图7.3.2图7.3.3图7.3.4五 心得体会 这次的课程设计是我们对数字信号处理这门课程知识综合应用的实践训练，也是是我们进入社会，从事自身工作前一个必不少的过程学习就像建房子，必须打好基础，一块砖一片瓦的累积起来，通过这次课程设计，我深深体会到这句话的含义只有在每

32、次的课程设计的实践当中去发现自己的问题，提高自己的能力才能有更大的收获。 在这次的课程设计中，我遇到了很多问题，发现自己以前的专业知识学习的不是很牢固，由于数字信号处理学的不是很好，很多相关知识的运用不是很熟悉，最主要的事MATLAB 的很多相关指令的运用不熟练，导致课程设计做起来有点困难，但是通过查阅资料以及询问老师最终还是解决了很多的难题，通过这次数字信号处理的课程设计，不但巩固了以前学的知识，还学习了新内容，收获颇丰。原来滤去语音中的噪声并不是像我想象中的那么遥远那么困难，而且更加巩固了MATLAB软件的使用，在从这次的课程设计中体会到，我觉得无论什么事都要自己认认真真去做了，才能使自己

33、成长，没有学习没有实践就不可能有提升能力的机会，希望这次的经历能让我们在以后的学习生活中不断成长，在以后的学习中更能与实际相结合，为以后的工作打下良好的基础。 从书本到实验，从理论到实践，在这段日子里，可以说得是苦多于甜，但是可以学到很多很多的东西，同时不仅可以巩固了以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能真正为社会服务，从而提高自己的实际动手能力和独立思考的能力。在设计的过程中遇到问题，可以说得是困难重重，但可喜的是最终都得到了解决。

34、 这次课程设计的每个实验细节和每个数据，都离不开老师的细心指导，和同学的热心帮助，让我明白了一个人的力量是渺小的，我们不能离群索居，只有大家互相帮助，相互鼓励，才能有更大的提升。参考文献（1）数字信号处理（第二版），丁玉美等，西安电子科技大学出版社；（2）数字信号处理试验指导书王创新、文卉编 长沙理工大学印刷（内部使用）（3）数字信号处理及其MATLAB实现，陈怀琛等译，电子工业出版社；（4）MATLAB及在电子信息课程中的应用，陈怀琛等，电子工业出版社。（5）数字信号处理A.V.奥本海姆，R.W.谢弗著，北京：科学出版社（6）数字信号处理理论、算法与实现（第二版）胡广书编著，北京：电子工业出

35、版社附录A：巴特沃斯低通滤波器程序Fs=22050;x,FS,bits=wavread(C:WINDOWSMediair_begin.wav);x=x(:,1);figure(1);subplot(2,1,1)plot(x);sound(x,FS,bits);title(语音信号时域波形图)y=fft(x,3260);f=(FS/1630)*1:1630; subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);title(语音信号频谱图);t=0:length(x)-1;zs=0.002*cos(2*pi*10000*t/22050);zs0=0.002*co

36、s(2*pi*10000*t/22050000);figure(2);subplot(2,1,1)plot(zs0);title(噪声信号波形);zs1=fft(zs,1200); sound(x,FS,bits);subplot(2,1,2)plot(f(1:600),abs(zs1(1:600); title(噪声信号频谱);x1=x+zs;sound(x1,FS,bits);y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title(加入噪声后的信号波形);subplot(2,1,2);plot(f(1:600),abs(y1(1:60

37、0);title(加入噪声后的信号频谱);fp=3000,fs=3500; Fs=22050;rp=1;rs=10;wp=2*pi*fp/FS;ws=2*pi*fs/FS;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=buttord(wap,was,rp,rs,s);B,A=butter(N,wc,s);Bz,Az=bilinear(B,A,Fs1);figure(4);h,w=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h);title(巴特沃斯低通滤波器);xlabel(频率(HZ);ylabel(耗损（dB）);gr

38、id on;yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(yd);title(滤波后的信号波形);ydd=fft(yd,1200);subplot(2,1,2);plot(f(1:600),abs(ydd(1:600); title(滤波后的信号频谱);sound(yd,FS,bits);附录B：巴特沃斯高通滤波器程序Fs=22050; x,FS,bits=wavread(C:WINDOWSMediair_begin.wav); x=x(:,1); figure(1); subplot(2,1,1); plot(x);sound(x,FS,

39、bits); title(语音信号时域波形图) y=fft(x,3260); f=(FS/1630)*1:1630; subplot(2,1,2); plot(f(1:1630),abs(y(1:1630); title(语音信号频谱图); t=0:length(x)-1; zs0=0.002*cos(2*pi*100*t/22050); figure(2); subplot(2,1,1); plot(zs0);title(噪声信号波形); zs1=fft(zs0,1200); sound(zs0,FS,bits); subplot(2,1,2); plot(f(1:600),abs(zs1(

40、1:600); title(噪声信号频谱); x1=x+zs0; sound(x1,FS,bits); y1=fft(x1,1200); figure(3); subplot(2,1,1);plot(x1); title(加入噪声后的信号波形); subplot(2,1,2); plot(f(1:600),abs(y1(1:600); title(加入噪声后的信号频谱); fp=600;fs=400;Fs=22050; rp=1;rs=10; wp=2*pi*fp/Fs; ws=2*pi*fs/Fs; T=1;Fs1=1; wap=2*tan(wp/2); was=2*tan(ws/2); N

41、,wc=buttord(wap,was,rp,rs,s); B,A=butter(N,wc,high,s); Bz,Az=bilinear(B,A,Fs1); figure(4); h,w=freqz(Bz,Az,512,Fs1*22050); plot(w,abs(h); title(巴特沃斯高通滤波器); xlabel(频率HZ);ylabel(耗损dB); grid on; yd=filter(Bz,Az,x1); figure(5); subplot(2,1,1);plot(yd); title(滤波后信号波形); ydd=fft(yd,1200); subplot(2,1,2);plot(f(1:600),abs(ydd(1:600); title(滤波后信号频谱); sound(yd,FS,bits);附录C：巴特沃斯带通滤波器程序Fs=22050; x,FS,bits=wavread(C:WINDOWSMedi