数字信号实习报告正文.doc

语音信号的数字滤波处理(八) 目 录 1 设计目的及内容.1 1.1 课程设计背景1 1.2 课程设计目的11.3 设计内容1 1.3.1 预习题部分.1 1.3.2 设计题部分.22 课程设计原理.22.1 频谱分析原理.22.2 IIR设计原理. . . .32.3 FIR设计原理. .33 设计步骤及过程 . . . . .5 3.1 语音信号的处理.5 3.2 语音信号的频谱分析.6 3.3 设计数字滤波器和画出频率响应.6 3.4 回放语音信号.64 程序的调试和运行结果. . .6 4.1 采样定理的演示. .6 4.2 间接法设计IIR数字滤波器. .74.2.1 切比雪夫低通滤波器的设计.74.2.2 切比雪夫高通滤波器的设计.84.2.3 切比雪夫带通滤波器的设计.104.3 直接法设计FIR数字滤波器. .114.3.1 Blackman窗低通滤波器的设计.114.3.2 Blackman窗高通滤波器的设计.134.3.3 Blackman窗带通滤波器的设计.145 结论 .15参考文献.17附录.18 1 设计目的及内容1.1 课程设计背景MATLAB是一个数据分析和处理功能十分强大的工程实用软件，他的数据采集工具箱为实现数据的输入和输出提供了十分方便的函数和命令。MATLAB是MathWorks公司于1982年推出的一款功能强大、易于使用的高效数值计算和可视化软件，它为进行算法开发、数据计算、信号分析与可视化提供了交互式应用开发环境，主要包括基本数学计算、编程环境(M语言)、数据可视化、GUIDE等。数字滤波器在信号的过滤、检测和参数估计等方面起着重要的作用。信号往往夹杂着噪声及无用信号成分，必须将这些干扰成分滤除。数字滤波器对信号进行筛选，可通过特定频段的信号。现代滤波器的作用是从含有噪声的信号中估计出信号的某些特征或信号本身，一旦信号被估计出，那么估计出来的信号与原信号相比会有更高的信噪比。1.2 课程设计目的(1)掌握数字信号处理的基本概念，基本理论和基本方法。(2)熟悉离散信号和系统的时域特性。(3)掌握序列快速傅里叶变换方法。(4)学会MATLAB的使用，掌握MATLAB的程序设计方法。(5)掌握利用MATLAB对语音信号进行频谱分析。(6)掌握滤波器的网络结构。(7)掌握MATLAB设计IIR、FIR数字滤波器的方法和对信号进行滤波的方法。1.3 设计内容1.3.1 预习题部分对下面连续信号进行采样：学号，A为幅度因子，a为衰减因子，为模拟角频率，其中n为学号(例如，贺娜同学n=201)；要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示不同采样频率下对应信号的频谱。1.3.2 设计题部分（1）语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；（2）产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；（3）污染信号的频谱分析，画出被污染的语音信号时域波形和频谱；（4）根据有关的频谱特性，采用间接法设计IIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个IIR滤波器）a模拟滤波器类型：切比雪夫滤波器（低通、带通、高通）b总体要求：Matlab原程序+仿真波形+技术指标（5）根据有关的频谱特性，采用直接法设计FIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个FIR滤波器）a滤波器类型：布莱克曼(Blackman)窗（低通、带通、高通）b总体要求：Matlab原程序+仿真波形+技术指标+窗函数（6）用自己设计的这些滤波器分别对被不同噪声污染的信号进行滤波；（7）分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；（8）回放语音信号。2 课程设计原理2.1 频谱分析原理时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。对于给定的时域信号y，可以通过Fourier变换得到频域信息Y。Y可按下式计算 （3-1）式中，N为样本容量，为采样间隔。采样信号的频谱是一个连续的频谱，不可能计算出所有的点的值，故采用离散Fourier变换(DFT)，即 （3-2）式中，。但上式的计算效率很低，因为有大量的指数(等价于三角函数)运算，故实际中多采用快速Fourier变换(FFT)。其原理即是将重复的三角函数算计的中间结果保存起来，以减少重复三角函数计算带来的时间浪费。由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。2.2 IIR设计原理语音信号的频谱范围主要为1000HZ左右，因此，在设计低通滤波器时，应把噪声频谱设定在2000HZ以上，这样，通过低通滤波器，即可滤除噪声信号从而还原语音信号；在设计高通滤波器时，应把噪声设定在500HZ以内，以通过高通滤波器滤除低频的噪声信号，从而还原相对频较高的语音信号；在设计带通滤波器时，可把噪声设计在低于500HZ或高于2000HZ的频谱上，已通过带通滤波器还原带通范围内的语音信号。(1) 切比雪夫滤波器原理切比雪夫滤波器的幅频特性具有等波纹特性。它有两种形式：振幅特性在通带内是等波纹的，在阻带内是单调下降的切比雪夫I型；振幅特性在阻带内是等波纹的，在通带内是单调下降的切比雪夫II型。(2) 双线性变换法工作原理 双线性变换中数字域频率和模拟频率之间的非线性关系限制了它的应用范围，只有当非线性失真是允许的或能被忽略时，才能采用双线性变换法，通常低通、高通、带通和带阻等滤波器等具有分段恒定的频率特性，可以采用预畸变的方法来补偿频率畸变，因此可以采用双线性变换设计方法。(3) 脉冲响应不变法工作原理 冲激响应不变法遵循的准则是使数字滤波器的单位取样响应与参照的模拟滤波器的脉冲响应的取样值完全一样，即h(n)=ha(nT),其中T为取样周期。实际是由模拟滤波器转换成为数字滤波器，就是要建立模拟系统函数Ha（S）与数字系统函数H（z）之间的关系。脉冲响应不变法是从S平面映射到z平面，这种映射不是简单的代数映射，而是S平面的每一条宽为2/T的横带重复地映射到整个z平面。 2.3 FIR设计原理 有限脉冲响应数字滤波器具有一般FIR系统的一切特征，与IIR系统的特征相比在单位脉冲响应h(n)、差分方程形式、系统函数H(z)以及系统结构方面有明显的区别。另外，可以设定条件来保证FIR滤波器具有线性相位特征，但IIR滤波器很难实现线性相位。设一个因果有限脉冲响应数字滤波器的单位脉冲响应为h(n),n=0N-1,则其相位的系统函数H(z)如下： （1）上式两边求Z反变换可以得到滤波器系统的差分方程描述形式： （2）显然，上式右边也就是单位脉冲响应h(n)与输入信号x(n)的线性卷积h(n)*x(n)。上式说明，FIR数字滤波器额当前输出信号值仅仅与当前输入信号值和N-1个以前的输入信号值有关。由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。但是这些特性是以牺牲线性相位频率特性为代价的，即用Butterworth、切比雪夫和椭圆法设计的数字滤波器逼近理想的滤波器的幅度频率特性，得到的滤波器往往是非线性的。在许多电子系统中，对幅度频率特性和线性相位特性都有较高的要求，所以IIR滤波器在这些系统中往往难以胜任。有限长单位冲激响应(FIR)数字滤波器具有以下优良的特点：(1) 可在设计任意幅度频率特性滤波器的同时，保证精确、严格的线性相位特性。(2) FIR数字滤波器的单位冲激响应h(n)是有限长的，可以用一个固定的系统来实现，因而FIR数字滤波器可以做成因果稳定系统。 (3) 允许设计多通带(多阻带)系统。 窗函数法就是设计FIR数字滤波器的最简单的方法。它在设计FIR数字滤波器中有很重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数4有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、海明窗(Hamming window)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window)。在MATLAB中，实现凯塞窗的函数为kaiser，调用格式为： w=kaiser (N,beta)其中beta为窗函数的参数。各种窗函数的性能比较可见表1，在设计FIR滤波器的过程中可以根据要求选择合适的窗函数：表1 各种窗函数的性能比较窗 函 数第一旁瓣相对于主瓣衰减/dB主 瓣 宽阻带最小衰减/dB矩形窗134/N21三角窗258/N25汉宁窗318/N44海明窗418/N53布拉克曼窗5712/N74凯塞窗可调可调可调切比雪夫窗可调可调可调3 设计步骤及过程3.1 语音信号的采集 利用windows下的录音机（开始程序附件娱乐录音机，文件属性立即转换8000Hz，8位，单声道）录制一段自己的话音，或者采用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,N2)，读取从N1点到N2点的采样值放在向量y中。3.2 语音信号的频谱分析 要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。3.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，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB。3.4 回放语音信号 在Matlab中，利用函数sound可以对声音进行回放，感受滤波前后的声音有变化。其调用格式：sound(x,fs,bits)。4 程序的调试和运行结果4.1 采样定理的演示采样序列的频谱是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为400Hz时频谱混叠很小，当采样频率为50Hz时频谱混叠很严重。要使采样信号可以恢复到原信号，采样频率必须满足时域采样定理，从而验证采样定理。图4.1 采样定理的演示4.2 间接法设计IIR数字滤波器4.2.1 切比雪夫低通滤波器的设计将语音信号和噪音信号叠加可以得到含噪声信号，噪声信号的频率设为7000HZ，分析其频谱特性之后用切比雪夫低通滤波器进行滤波。用MATLAB仿真，仿真处理如图所示。图4.2 语音信号的时域波形和频谱图 图4.3 噪声信号的时域波形和频谱图 图4.4 被污染后信号的时域波形和频谱图 图4.5 低通滤波器的幅频特性图图4.6 滤波后信号的时域波形和频谱图4.2.2 切比雪夫高通滤波器的设计将语音信号和噪音信号叠加可以得到含噪声信号，噪声信号的频率设为7000HZ，分析其频谱特性之后用切比雪夫高通滤波器进行滤波。用MATLAB仿真，仿真处理如图所示。 图4.7 语音信号的时域波形和频谱图 图4.8 噪声信号的时域波形和频谱图 图4.9 被污染后信号的时域波形和频谱图 图4.10 高通滤波器的幅频特性图图4.11 滤波后信号的时域波形和频谱图4.2.3 切比雪夫带通滤波器的设计将语音信号和噪音信号叠加可以得到含噪声信号，设两个噪声信号，其频率分别为100HZ和2200HZ，分析其频谱特性之后用切比雪夫带通滤波器进行滤波。用MATLAB仿真，仿真处理如图所示。 图4.12 语音信号的时域波形和频谱图 图4.13 噪声信号的时域波形和频谱图 图4.14 被污染后信号的时域波形和频谱图 图4.15 带通滤波器的幅频特性图图4.16 滤波后信号的时域波形和频谱图4.3 直接法设计FIR数字滤波器4.3.1 Blackman窗低通滤波器的设计将语音信号和噪音信号叠加可以得到含噪声信号，噪声信号的频率设为7000HZ，分析其频谱特性之后用Blackman低通滤波器进行滤波。用MATLAB仿真，仿真处理如图所示。 图4.17 语音信号的时域波形和频谱图 图4.18 噪声信号的时域波形和频谱图 图4.19 被污染后信号的时域波形和频谱图 图4.20 blackman低通滤波器的幅频特性图图4.21 滤波后信号的时域波形和频谱图4.3.2 Blackman窗高通滤波器的设计将语音信号和噪音信号叠加可以得到含噪声信号，噪声信号的频率设为100HZ，分析其频谱特性之后用Blackman高通滤波器进行滤波。用MATLAB仿真，仿真处理如图所示。 图4.22 语音信号的时域波形和频谱图 图4.23 噪声信号的时域波形和频谱图 图4.24 被污染后信号的时域波形和频谱图 图4.25 blackman高通滤波器的幅频特性图图4.26 滤波后信号的时域波形和频谱图4.3.3 Blackman窗带通滤波器的设计将语音信号和噪音信号叠加可以得到含噪声信号，噪声信号的频率设为7000HZ，分析其频谱特性之后用切比雪夫带通滤波器进行滤波。用MATLAB仿真，仿真处理如图所示。 图4.27 语音信号的时域波形和频谱图 图4.28 噪声信号的时域波形和频谱图 图4.29 被污染后信号的时域波形和频谱图 图4.30 blackman带通滤波器的幅频特性图图4.31 滤波后信号的时域波形和频谱图5 结论 通过本次为期一周的数字信号处理的课程设计，使我熟悉和巩固了切比雪夫滤波器和Blackman滤波器的设计方法和原理，能实现滤波器设计的有关经典算法，更重要的是熟练掌握使用MATLAB语言设计各种要求的数字滤波器。这一部分的内容相对来说是有些难度了，做起来花费的精力也多了一些，不过对数字信号处理内容的掌握上又加深了一层。这次的课程设计也让我体会到了能读懂程序和自己编写程序的差距是很大的，理论与实际相结合是很重要的。在之前没有做这次系统的课设之前 ，我懂得的只是书本上的理论知识，但经过这次的课程设计让我对这些滤波器有了更深入的了解，同时我也意识到只有理论知识是远远不够的，要把所学的理论知识与实践相结合起来，经过大量的亲身实践后才能够得到更深的掌握，最终达到融会贯通的效果。参考文献1 高西全，丁玉美数字信号处理（第三版）M西安：西安电子科技大学出版社，2008：167-168，211-2142 王创新，文卉数字信号处理试验指导书长沙理工大学印刷（内部使用）3 陈怀琛数字信号处理及其MATLAB实现北京：电子工业出版社4 陈怀琛MATLAB及在电子信息课程中的应用北京：电子工业出版社5 奥本海姆，谢弗著数字信号处理北京：科学出版社6 胡广书数字信号处理理论、算法与实现（第二版）北京：电子工业出版社附录 源程序%采样定理的演示t=0:0.005:1 xa=108*exp(-108*sqrt(2)*pi*t).*sin(108*sqrt(2)*pi*t);subplot(4,1,1);plot(t,xa);xlabel( 时间u(t)（s）);ylabel( 幅度xa(t); title(信号波形)t=0.1; n=0:t:1; subplot(4,1,2);stem(n);xlabel( 采样点（N）);ylabel(幅度xa(t); title(采样脉冲)y1=fft(xa,400); subplot(4,1,3);plot(abs(y1);xlabel( 频率(HZ）);ylabel( 幅度xa(t); title(采样频率为400时的频谱)y1=fft(xa,50); subplot(4,1,4);plot(abs(y1);xlabel(频率（HZ）);ylabel(幅度xa(t);title(采样频率为50时的频谱)%切比雪夫低通滤波器clf;x,Fs,bits=wavread(H:数字信号实习yuyin.wav);t=(0:length(x)-1)/44100;x=x(:,1); figure(1);subplot(2,1,1);plot(t,x); sound(x,Fs,bits);xlabel(时间（s）);ylabel(幅度(db);title(语音信号时域波形图); %采样后语音信号的时域波形图y=fft(x,3260); f=(Fs/3260)*1:1630;subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(语音信号频谱图); %采样后语音信号的频谱图d=0.1*cos(2*pi*7000*t);%产生污染信号figure(2);subplot(2,1,1)plot(t,d);xlabel(时间（s）);ylabel(幅度(db);title(噪声信号波形);%污染信号波形d1=fft(d,3260);%sound(d,Fs,bits);subplot(2,1,2)plot(f(1:1000),abs(d1(1:1000);xlabel(频率（HZ）);ylabel(幅度(db);title(噪声信号频谱);%污染信号频谱x1=x+d;%给语音信号加污染信号y1=fft(x1,3260); figure(3); subplot(2,1,1);plot(t,x1); xlabel(时间（s）);ylabel(幅度(db);title(被污染后的信号波形);%被污染的语音信号时域波形图subplot(2,1,2); plot(f(1:1024),abs(y1(1:1024);xlabel(频率（HZ）);ylabel(幅度(db);title(被污染后的信号频谱);%被污染的语音信号频谱图fp=1000;fs=1500;rp=1;%设置低通滤波器通带内的波纹为1dbas=10;%设置低通滤波器阻带内的波纹为1dbwp=2*pi*1000/Fs;ws=2*pi*1500/Fs;N,wc=cheb1ord(wp,ws,rp,as,s); %计算切比雪1型模拟低通滤波器阶数和通带边界频率B,A=cheby1(N,rp,wc,s); %计算切比雪1型模拟低通滤波器系统函数系数Bz,Az=bilinear(B,A,Fs1);figure(4);h,w=freqz(Bz,Az,length(x),Fs);plot(w,abs(h);grid onxlabel(频率（HZ）);ylabel(幅度(db);title(低通滤波器幅频特性);yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(t,yd);title(滤波后信号波形);xlabel(时间（s）);ylabel(幅度(db);ydd=fft(yd,3260);subplot(2,1,2);plot(f(1:1630),abs(ydd(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(滤波后信号频谱);sound(yd,Fs,bits);%切比雪夫高通滤波器clf;x,Fs,bits=wavread(H:数字信号实习yuyin.wav);t=(0:length(x)-1)/44100;x=x(:,1); figure(1);subplot(2,1,1);plot(t,x); sound(x,Fs,bits);xlabel(时间（s）);ylabel(幅度(db);title(语音信号时域波形图); %采样后语音信号的时域波形图y=fft(x,3260); f=(Fs/3260)*1:1630;subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(语音信号频谱图); %采样后语音信号的频谱图d=0.1*cos(2*pi*100*t);%产生污染信号figure(2);subplot(2,1,1)plot(t,d);xlabel(时间（s）);ylabel(幅度(db);title(噪声信号波形);%污染信号波形d1=fft(d,3260);%sound(d,Fs,bits);subplot(2,1,2)plot(f(1:1000),abs(d1(1:1000);xlabel(频率（HZ）);ylabel(幅度(db);title(噪声信号频谱);%污染信号频谱x1=x+d;%给语音信号加污染信号y1=fft(x1,3260); figure(3); subplot(2,1,1);plot(t,x1); xlabel(时间（s）);ylabel(幅度(db);title(被污染后的信号波形);%被污染的语音信号时域波形图subplot(2,1,2); plot(f(1:1024),abs(y1(1:1024);xlabel(频率（HZ）);ylabel(幅度(db);title(被污染后的信号频谱);%被污染的语音信号频谱图fp=1000;fs=500;rp=1;%设置高通滤波器通带内的波纹为1dbrs=100;%设置高通滤波器阻带内的波纹为10dbwp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=cheb1ord(wap,was,rp,rs,s); %计算切比雪1型模拟高通滤波器阶数和通带边界频率B,A=cheby1(N,rp,wc,high,s); %计算切比雪1型模拟高通滤波器系统函数系数Bz,Az=bilinear(B,A,Fs1);figure(4);h,w=freqz(Bz,Az,length(x),Fs);plot(w,abs(h);grid onxlabel(频率（HZ）);ylabel(幅度(db);title(高通滤波器幅频特性);yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(t,yd);title(滤波后信号波形);xlabel(时间（s）);ylabel(幅度(db);ydd=fft(yd,3260);subplot(2,1,2);plot(f(1:1630),abs(ydd(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(滤波后信号频谱);sound(yd,Fs,bits);%切比雪夫带通滤波器clf;x,Fs,bits=wavread(H:数字信号实习yuyin.wav);t=(0:length(x)-1)/44100;x=x(:,1); figure(1);subplot(2,1,1);plot(t,x); sound(x,Fs,bits);xlabel(时间（s）);ylabel(幅度(db);title(语音信号时域波形图); %采样后语音信号的时域波形图y=fft(x,3260); f=(Fs/3260)*1:1630;subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(语音信号频谱图); %采样后语音信号的频谱图d=0.1*cos(2*pi*100*t);%产生污染信号ddx=0.1*cos(2*pi*2200*t);%产生污染信号dxfigure(2);subplot(2,2,1)plot(t,d);xlabel(时间（s）);ylabel(幅度(db);title(噪声d信号波形);%污染信号d波形d1=fft(d,3260);subplot(2,2,2)plot(f(1:1000),abs(d1(1:1000);xlabel(频率（HZ）);ylabel(幅度(db);title(噪声信号d频谱);%污染信号d频谱subplot(2,2,3)plot(t,dx);xlabel(时间（s）);ylabel(幅度(db);title(噪声信号dx波形);%污染信号dx波形d2=fft(dx,3260);subplot(2,2,4)plot(f(1:1000),abs(d2(1:1000);xlabel(频率（HZ）);ylabel(幅度(db);title(噪声信号dx频谱);%污染信号dx频谱x1=x+d+dx;%给语音信号加污染信号y1=fft(x1,3260); figure(3); subplot(2,1,1);plot(t,x1); xlabel(时间（s）);ylabel(幅度(db);title(被污染后的信号波形);%被污染的语音信号时域波形图subplot(2,1,2); plot(f(1:1024),abs(y1(1:1024);xlabel(频率（HZ）);ylabel(幅度(db);title(被污染后的信号频谱);%被污染的语音信号频谱图fp1=1000;fp2=1500;fs1=500;fs2=2000;rp=1;%设置带通滤波器通带内的波纹为1dbrs=10;%设置带通滤波器阻带内的波纹为10dbwp=2*pi*fp1,fp2/Fs;ws=2*pi*fs1,fs2/Fs;Fs1=Fs/Fs;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=cheb1ord(wap,was,rp,rs,s); %计算切比雪1型模拟带通滤波器阶数和通带边界频率B,A=cheby1(N,rp,wc,s); %计算切比雪1型模拟带通滤波器系统函数系数Bz,Az=bilinear(B,A,Fs1);figure(4);h,w=freqz(Bz,Az,length(x),Fs);plot(w,abs(h);grid onxlabel(频率（HZ）);ylabel(幅度(db);title(带通滤波器幅频特性);yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(t,yd);title(滤波后信号波形);xlabel(时间（s）);ylabel(幅度(db);ydd=fft(yd,3260);subplot(2,1,2);plot(f(1:1630),abs(ydd(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(滤波后信号频谱);sound(yd,Fs,bits);%Blackman低通滤波器clf;x,Fs,bits=wavread(H:数字信号实习yuyin.wav);t=(0:length(x)-1)/Fs;x=x(:,1); figure(1);subplot(2,1,1);plot(t,x); sound(x,Fs,bits);xlabel(时间（s）);ylabel(幅度(db);title(语音信号时域波形图); %采样后语音信号的时域波形图y=fft(x,3260); f=(Fs/3260)*1:1630;subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);xlabel(频率（HZ）);ylabel(幅度(db);title(语音信号频谱图); %采样后语音信号的频谱图d=0.1*cos(2*pi*7000*t);%产生污染信号figure(2);subplot(2,1,1)plot(t,d);xlabel(时间（s）);ylabel(幅度(db);title(噪声信号波形);%污染信号波形d1=fft(d,3260);subplot(2,1,2)plot(f(1:1000),abs(d1(1:1000);xlabel(频率（HZ）);ylabel(幅度(db);title(噪声信号频谱);%污染信号频谱x1=x+d;%给语音信号加污染信号y1=fft(x1,3260); figure(3); subplot(2,1,1);plot(t,x1); xlabel(时间（s）);ylabel(幅度(db);title(被污染