数字信号处理课程设计基于FIR的有噪声语音信号处理.doc

1 课程设计的主要目的和内容及要求课程设计目的：1. 通过本次课程设计，综合运用数字信号处理技术课程和其他有关先修课程的理论和生产实际知识去分析和解决具体问题，并使所学知识得到进一步巩固、深化和发展。2. 初步培养学生对工程设计的独立工作能力，学习设计的一般方法。3. 通过课程设计树立正确的设计思想，提高学生分析问题、解决问题的能力。4. 进行设计基本技能的训练，如查阅设计资料和手册、程序的设计、调试等。本课题设计内容：1. 利用Windows下的录音机或其他软件，选择Windows系统的“叮······”（ding.wav），并对该信号进行采样；2. 语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；3. 产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；4. 污染信号的频谱分析，画出被污染的语音信号时域波形和频谱； 5. 根据有关的频谱特性，采用Kaiser窗设计FIR数字低通，高通，带通，带阻滤波器并画出4种滤波器的幅频、相频图；6. 用自己设计的4种滤波器分别对被不同噪声污染的信号进行滤波；7. 分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；8. 回放语音信号。设计要求：1. 熟悉离散信号和系统的时域特性。2. 掌握数字信号处理的基本概念，基本理论和基本方法。3. 掌握序列快速傅里叶变换方法。4. 学会MATLAB的使用，掌握MATLAB的程序设计方法。5. 掌握利用MATLAB对语音信号进行频谱分析。6. 掌握MATLAB设计FIR数字滤波器的方法和对信号进行滤波的方法。2 课程设计的总体方案2.1设计框图加入噪声信号noise4加入噪声信号noise3加入噪声信号noise2加入噪声信号noise1信号采样并画频域图画信号的时域图 获取原始信号x xxxx x x x对信号采样画频域图画信号的时域图对信号采样画频域图画信号的时域图对信号采样画频域图画信号的时域图对信号采样画频域图画信号的时域图低通滤波器带阻滤波器带通滤波器高通滤波器画y4的时域图并采样画频域图画y2的时域图并采样画频域图画y1的时域图并采样画频域图画y3的时域图并采样画频域图 低通滤波器幅频与相频图带通滤波器幅频与相频图带阻滤波器幅频与相频图高通滤波器幅频与相频图2.2总体设计概述本次课程设计所采用的主要软件为Matlab。设计中用到的函数主要有wavread,sound,fft,subplot,plot,fir1,freqz,filter。设计的主要内容为数字滤波器的设计。从功能上分可把滤波器分为低通，高通，带阻，带通滤波器。同时数字滤波器的设计又可分为无限脉冲响应（IIR）数字滤波器设计和有限脉冲响应（FIR）数字滤波器。IIR数字滤波器设计方法主要可分为脉冲响应不变法和双线性变换法，其中主要有巴特沃斯滤波器，切比雪夫I型滤波器，切比雪夫II型滤波器，椭圆滤波器以及贝塞尔滤波器等五种模拟滤波器做为原型；而对于FIR数字滤波器设计方法主要有窗函数法，频率采样法和等波纹最佳逼近法，其中窗函数法中主要有六种典型的窗函数，它们分别是矩形窗，三角形窗，汉宁（Hanning）窗，哈明（Hamming）窗，布莱克曼（Blackman）窗以及凯塞贝塞尔窗（Kaiser-Basel Window）。本课程设计采用的是FIR数字滤波器中的窗函数法，选用的是凯塞贝塞尔窗.主要的设计思路为首先获取原始信号x。其次在画出x的时域图和经采样后的频域图后分别加入四种噪声信号noise1,noise2,noise3,noise4分别获得X1,X2,X3,X4四种带噪声信号。然后分别画出X1,X2,X3,X4的时域图并进行采样画出频域图。再次，利用Matlab中的函数fir1,Kaiser,freqz,filter分别设计出低通，高通，带阻，带通四种数字滤波器并分别画出它们的幅频特性和相频特性。最后分别让X1，X2，X3，X4四种信号分别通过四种滤波器得到四种滤波后的信号y1,y2,y3,y4，分别画出y1,y2,y3,y4四种信号的时域图和经采样后的频域图与X1,X2,X3,X4进行比较。3 设计的详细原理设计主要可分为原始信号部分，噪声信号部分，低通滤波部分，高通滤波部分，带阻滤波部分，带通滤波部分等六个部分。3.1原始信号部分本次使用的原始信号是windows系统的“叮······”（ding.wav）语音信号，可利用Windows下的录音机（wavread）或其他软件进行获取。利用sound函数进行试听。其中wavread的主要参数有采样频率（Fs），采样信号（x），采样位数（bits），本次设计中Fs10000，bits16从而可获得时间长为1s的语音信号x。sound函数调用格式为sound(x,Fs,bits) 。对x信号进行采样y=fft(x,m)，其中m=length(x)。最后利用plot函数分别画出x信号的时域图和频域图subplot(211);plot(t,x);title('原始信号波形');xlabel('time(s)');subplot(212);plot(f,abs(y);title('原始信号频谱');xlabel('frequency(hz)');便可完成原始信号部分的设计。同时通过图像可观察x的特性，其频率在800Hz左右。3.2噪声信号部分本次课程设计中主要有四种噪声信号，它们分别是高频噪声信号noise1=0.01*sin(2*pi*6000*t) 低频噪声信号noise2=0.01*sin(2*pi*200*t) 带阻噪声信号noise3=0.01*sin(2*pi*4000*t)+0.01*sin(2*pi*5000*t) 带通噪声信号noise4=0.01*sin(2*pi*5000*t)+0.01*sin(2*pi*200*t) 。其中noise1为一个频率为5000Hz的高频信号可作为低通滤波中x的噪声输入得到X1信号，noise2为频率为200Hz的低频信号可作为高通滤波中x的噪声输入得到X2信号，noise3为频率分别为4000Hz和5000Hz的带阻信号可作为带阻滤波中x的噪声输入得到X3信号，noise4为频率分别为200Hz和5000Hz的带通信号可作为带通滤波中x的噪声输入得到X4信号。3.3低通滤波部分低通滤波部分可分两小部分：噪声部分和低通滤波器设计部分。噪声部分低通滤波的噪声信号为noise1，使X1=x+noise1'获得噪声污染后的信号X1，并利用sound回放信号X1与原始信号x进行比较。然后利用plot函数画X1的时域图和经采样后的频域图，与原信号时域图和频域图进行比较。低通滤波器设计部分所谓低通滤波器是当信号频率（w）大于通带边界频率（wp）时将无法通过此滤波器，而当w<wp时信号就能通过。低通滤波器设计的主要参数有通带边界频率wp1，阻带截止频率ws1，3dB通带截止频率wc1。其中wp12*pi*fp1/Fs，ws12*pi*fs1/Fs（其中fp1，fs1为模拟频率）wc1(wp1+ws1)/2。本次设计中fp12000，fs13000。选定参数后利用b1=fir1(N,wc1,kaiser(N+1) ，hn1 w1=freqz(b1,1,512)语句实现低通滤波器的设计，其中N为滤波器的阶数，本次设计中假定N=38。在通过x1get=filter(b1,1,X1); 语句进行对X1信号的滤波得到信号x1get，再利用Y1get=fft(x1get,m);语句进行x1get的采样。同时运用subplot(323);plot(w1/pi,20*log10(abs(hn1);title('低通滤波器幅频');xlabel('w'); subplot(324);plot(w1/pi,unwrap(angle(hn1);title('低通滤波器相频');xlabel('w');语句绘出滤波器的幅频和相频图，并同时画出x1get的时域图和频域图与X1信号进行比较。3.4高通滤波部分高通滤波部分可分两小部分：噪声部分和高通滤波器设计部分。噪声部分高通滤波的噪声信号为noise2，使X2=x+noise1'获得噪声污染后的信号X2，并利用sound回放信号X2与原始信号x进行比较。然后利用plot函数画X2的时域图和经采样后的频域图，与原信号时域图和频域图进行比较。高通滤波器设计部分所谓高通滤波器是当信号频率（w）小于通带边界频率（wp）时将无法通过此滤波器，而当w>wp时信号便可通过。高通滤波器设计的主要参数有通带边界频率wp2，阻带截止频率ws2，3dB通带截止频率wc2。其中wp22*pi*fp2/Fs，ws22*pi*fs2/Fs（其中fp2，fs2为模拟频率）wc2(wp2+ws2)/2。本次设计中fp2500，fs2300。选定参数后利用b2=fir1(N,wc2,'high',kaiser(N+1) ，hn2 w2=freqz(b2,1,512)语句实现低通滤波器的设计,其中N为滤波器的阶数，本次设计中假定N=38。在通过x2get=filter(b2,1,X2); 语句进行对X2信号的滤波得到信号x2get，再利用Y2get=fft(x2get,m);语句进行x2get的采样。同时运用subplot(323);plot(w2/pi,20*log10(abs(hn2);title('高通滤波器幅频');xlabel('w');subplot(324);plot(w2/pi,unwrap(angle(hn2);title('高通滤波器相频');xlabel('w');语句绘出滤波器的幅频和相频图，并同时画出x2get的时域图和频域图与X2信号进行比较。3.5带阻滤波部分带阻滤波部分可分两小部分：噪声部分和带阻滤波器设计部分。噪声部分带阻滤波的噪声信号为noise3，使X3=x+noise3'获得噪声污染后的信号X3，并利用sound回放信号X3与原始信号x进行比较。然后利用plot函数画X3的时域图和经采样后的频域图，与原信号时域图和频域图进行比较。带阻滤波器设计部分所谓带阻滤波器是当信号频率（w）大于通带边界频率（wpl）小于通带边界频率（wph）时将无法通过此滤波器，而当w>wph或w<wpl时信号便可通过。带阻滤波器设计的主要参数有通带边界频率wp3l，wp3h，阻带截止频率ws3l，ws3h，3dB通带截止频率wc3。其中wp3=2*fpl3/Fs 2*fph3/Fs;ws3=2*fsl3/Fs 2*fsh3/Fs;（其中fp3，fs3为模拟频率）wc3(wp3+ws3)/2。本次设计中fpl32000，fph38000，fsl33000，fsh37000。选定参数后利用b3=fir1(N,wc3,'stop',kaiser(N+1)，hn3 wp3=freqz(b3,1,512)语句实现低通滤波器的设计,其中N为滤波器的阶数，本次设计中假定N=38。在通过x3get=filter(b3,1,X3); 语句进行对X3信号的滤波得到信号x3get，再利用Y3get=fft(x3get,m);语句进行x3get的采样。同时运用subplot(323);plot(w2/pi,20*log10(abs(hn3);title('带阻滤波器幅频');xlabel('w');subplot(324);plot(w2/pi,unwrap(angle(hn3);title('带阻滤波器相频');xlabel('w');语句绘出滤波器的幅频和相频图，并同时画出x3get的时域图和频域图与X3信号进行比较。3.6带通滤波部分带通滤波部分可分两小部分：噪声部分和带通滤波器设计部分。噪声部分带通滤波的噪声信号为noise4，使X4=x+noise4'获得噪声污染后的信号X4，并利用sound回放信号X4与原始信号x进行比较。然后利用plot函数画X4的时域图和经采样后的频域图，与原信号时域图和频域图进行比较。带通滤波器设计部分所谓带通滤波器是当信号频率（w）小于通带边界频率（wpl）大于通带边界频率（wph）时将无法通过此滤波器，而当wpl<w<wph时信号便可通过。带通滤波器设计的主要参数有通带边界频率wp4l，wp4h，阻带截止频率ws4l，ws4h，3dB通带截止频率wc4。其中wp4=2*fpl4/Fs 2*fph4/Fs;ws4=2*fsl4/Fs 2*fsh4/Fs;（其中fp4，fs4为模拟频率）wc4(wp4+ws4)/2。本次设计中fpl4500，fph42000，fsl4700，fsh41800。选定参数后利用b4=fir1(N,wc4,kaiser(N+1)，hn4 wp4=freqz(b4,1,512)语句实现低通滤波器的设计,其中N为滤波器的阶数，本次设计中假定N=38。在通过x4get=filter(b4,1,X4); 语句进行对X4信号的滤波得到信号x4get，再利用Y4get=fft(x4get,m);语句进行x4get的采样。同时运用subplot(323);plot(wp4/pi,20*log10(abs(hn4);title('带通滤波器幅频');xlabel('w');subplot(324);plot(wp4/pi,unwrap(angle(hn4);title('带通滤波器相频');xlabel('w');语句绘出滤波器的幅频和相频图，并同时画出x4get的时域图和频域图与X4信号进行比较。4 设计的步骤和过程本次课程设计的步骤为原始信号部分噪声信号部分低通滤波部分高通滤波部分带阻滤波部分带通滤波部分完成设计并运行4.1原始信号部分的设计程序如下：clf;%原始音频抽样Fs=10000;bits=16;T=1;ts=1/Fs;N=38;x,Fs,bits=wavread('ding.wav'); %利用Windows下的录音机，录制一段自己的语音信号，时间为1s x=x(:,1); m=length(x);sound(x,Fs,bits);y=fft(x,m); %对录制的信号进行傅立叶变换f=(Fs/m)*1:m;t=1:m/Fs;subplot(211);plot(t,x);title('原始信号波形');xlabel('time(s)');subplot(212);plot(f,abs(y);title('原始信号频谱');xlabel('frequency(hz)');4.2噪声信号部分的设计程序如下：%噪声信号noise1=0.01*sin(2*pi*5000*t); %高频噪声noise2=0.01*sin(2*pi*200*t); %低频噪声noise3=0.01*sin(2*pi*4000*t)+0.01*sin(2*pi*5000*t); %带阻噪声noise4=0.01*sin(2*pi*5000*t)+0.01*sin(2*pi*200*t); %带通噪声4.3低通滤波部分的设计程序如下：% 低通滤波 fp1=2000; fs1=3000; wp1=2*pi*fp1/Fs; ws1=2*pi*fs1/Fs; X1=x+noise1'%加了高频噪声的信号 sound(X1,Fs,bits);%回放噪声信号 Y11=fft(X1,m);%对加了高频噪声的信号进行傅立叶变换 figure; subplot(321);plot(t,X1);title('加了高频噪声的信号时域图');xlabel('time(s)'); subplot(322);plot(f,abs(Y11);title('加了高频噪声的信号频谱图');xlabel('frequency(hz)'); wc1=(wp1+ws1)/2; b1=fir1(N,wc1,kaiser(N+1);%调用kaiser计算低通FIRDF的b1（n） hn1 w1=freqz(b1,1,512); subplot(323);plot(w1/pi,20*log10(abs(hn1);title('低通滤波器幅频');xlabel('w'); subplot(324);plot(w1/pi,unwrap(angle(hn1);title('低通滤波器相频');xlabel('w'); x1get=filter(b1,1,X1); Y1get=fft(x1get,m); sound(x1get,Fs,bits); subplot(325);plot(t,x1get);title('加了高频噪声的信号滤波后时域图');xlabel('time(s)');subplot(326);plot(f,abs(Y1get);title('加了高频噪声的信号滤波后频谱图');xlabel('frequency(hz)');4.4高通滤波部分的设计程序如下：%高通滤波fp2=500;fs2=300;wp2=2*fp2/Fs;ws2=2*fs2/Fs;X2=x+noise2'%加了低频噪声的信号sound(X2,Fs,bits);%回放噪声信号Y21=fft(X2,m);%对加了低频噪声的信号进行傅立叶变换figure;subplot(321);plot(t,X2);title('加了低频噪声的信号时域图');xlabel('time(s)');subplot(322);plot(f,abs(Y21);title('加了低频噪声的信号频谱图');xlabel('frequency(hz)');wc2=(wp2+ws2)/2;b2=fir1(N,wc2,'high',kaiser(N+1);%调用kaiser计算高通FIRDF的b2（n）hn2 w2=freqz(b2,1,512);subplot(323);plot(w2/pi,20*log10(abs(hn2);title('高通滤波器幅频');xlabel('w');subplot(324);plot(w2/pi,unwrap(angle(hn2);title('高通滤波器相频');xlabel('w');x2get=filter(b2,1,X2);Y2get=fft(x2get,m);sound(x2get,Fs,bits);subplot(325);plot(t,x2get);title('加了低频噪声的信号滤波后时域图');xlabel('time(s)');subplot(326);plot(f,abs(Y2get);title('加了低频噪声的信号滤波后频谱图');xlabel('frequency(hz)');4.5带阻滤波部分的设计程序如下：%带阻滤波fpl3=2000;fph3=8000;fsl3=3000;fsh3=7000;wp3=2*fpl3/Fs 2*fph3/Fs;ws3=2*fsl3/Fs 2*fsh3/Fs;X3=x+noise3'%加了带阻噪声的信号sound(X3,Fs,bits);%回放噪声信号Y31=fft(X3,m);%对加了带阻噪声的信号进行傅立叶变换figure;subplot(321);plot(t,X3);title('加了带阻噪声的信号时域图');xlabel('time(s)');subplot(322);plot(f,abs(Y31);title('加了带阻噪声的信号频谱图');xlabel('frequency(hz)');wc3=(wp3+ws3)/2;b3=fir1(N,wc3,'stop',kaiser(N+1);%调用kaiser计算带阻FIRDF的b3（n）hn3 wp3=freqz(b3,1,512);subplot(323);plot(w2/pi,20*log10(abs(hn3);title('带阻滤波器幅频');xlabel('w');subplot(324);plot(w2/pi,unwrap(angle(hn3);title('带阻滤波器相频');xlabel('w');x3get=filter(b3,1,X3);Y3get=fft(x3get,m);sound(x3get,Fs,bits);subplot(325);plot(t,x3get);title('加了带阻噪声的信号滤波后时域图');xlabel('time(s)');subplot(326);plot(f,abs(Y3get);title('加了带阻噪声的信号滤波后频谱图');xlabel('frequency(hz)');4.6带通滤波部分的设计程序如下：%带通滤波fpl4=500;fph4=2000;fsl4=700;fsh4=1800;wp4=2*fpl4/Fs 2*fph4/Fs;ws4=2*fsl4/Fs 2*fsh4/Fs;X4=x+noise4'%加了带通噪声的信号sound(X4,Fs,bits);%回放噪声信号Y41=fft(X4,m);%对加了带通噪声的信号进行傅立叶变换figure;subplot(321);plot(t,X4);title('加了带通噪声的信号时域图');xlabel('time(s)');subplot(322);plot(f,abs(Y41);title('加了带通噪声的信号频谱图');xlabel('frequency(hz)');wc4=(wp4+ws4)/2;b4=fir1(N,wc4,kaiser(N+1);%调用kaiser计算带通FIRDF的b4（n）hn4 wp4=freqz(b4,1,512);subplot(323);plot(wp4/pi,20*log10(abs(hn4);title('带通滤波器幅频');xlabel('w');subplot(324);plot(wp4/pi,unwrap(angle(hn4);title('带通滤波器相频');xlabel('w');x4get=filter(b4,1,X3);Y4get=fft(x4get,m);sound(x4get,Fs,bits);subplot(325);plot(t,x4get);title('加了带通噪声的信号滤波后时域图');xlabel('time(s)');subplot(326);plot(f,abs(Y4get);title('加了带通噪声的信号滤波后频谱图');xlabel('frequency(hz)');至此，设计的六个主体部分都基本完成，可完成数字滤波器设计，也完成X1、X2、X3、X4四个有噪声语音信号的滤波过程。5设计程序的调试及运行结果5.1设计程序调试运行程序，观察滤波后的信号频谱中的噪声成分是否去除。若滤波后信号频谱中的噪声成分没有去除干净，则观察程序进行调试，调试后再观察，直至噪声成分完全去除；若滤波后信号频谱中的噪声成分已经去除干净，则记录运行结果并进行结果分析。5.2程序运行结果图1.原始信号的时域图与频域图由上图观察可知原始信号的频率为800hz左右，同时又由于周期性的原因，图上才有两个频率成分。 图2.低通滤波部分运行结果图由噪声信号和滤波后的信号的时域图和频谱图可看出噪声信号中的噪声成分已成功的去除，因而设计成功由滤波器的幅频和相频图可看出滤波器的设计也符合要求。图3.高通滤波部分运行结果图由噪声信号和滤波后的信号的时域图和频谱图可看出噪声信号中的噪声成分已成功的去除，因而设计成功由滤波器的幅频和相频图可看出滤波器的设计也符合要求。图4.带阻滤波部分运行结果图由噪声信号和滤波后的信号的时域图和频谱图可看出噪声信号中的噪声成分已成功的去除，因而设计成功由滤波器的幅频和相频图可看出滤波器的设计也符合要求。图5.带通部分运行结果图由噪声信号和滤波后的信号的时域图和频谱图可看出噪声信号中的噪声成分已成功的去除，因而设计成功由滤波器的幅频和相频图可看出滤波器的设计也符合要求。5 课程设计心得体会时间过得很快，是一个说了上百回的句子，可再用上一千回，它还是对的。刚刚接到课程设计的任务，想着长达两周的时间，总以为时间有的是。只等到要交时原来自己还是蒙然未懂，一切的一切又是模糊不清。但在这过程中自己体会到了一些，一点乃至半点，这两周便是有收获。作课程设计已不是第一回了，有了经验但同时也让自己有些懒散。本次设计的课题时数字滤波器的设计。从第一个周一的选题开始就有些不以为然，感觉做完了实验，对这些滤波器的了解应该比较深。因而也就随便选了一个课题，因为个人觉得所有的课题都差不多。选完题当然就是做了。拿出书，打开电脑，开启MATLAB，接着举手，到最后便无奈的放下手，设计的第一步就让自己感到棘手。原始信号的获得，老师说过可用windows系统下的录音机或软件，可是两个自己都没用过。察书，没有，自己借的书有的只是纯理论上的东西，与实际运用差了十万八千里，没办法最后用了一个最有效但一开始没想到的办法上网。还真别说，有效啊，确实有效，查到了wavread和sound两个成对的函数。心里想着，这网络被人说成是毒遍天下，是杀人无形的兵刃，但看来用好了倒也是个百科书林。原始信号有了，自己也得了教训，便开始计划开来。自己选的课题是运用凯塞窗函数设计一个FIR数字滤波器。于是，翻书看完有关FIR滤波器的相关介绍，重点看了凯塞窗。然后查资料找出与本次设计有关的MATLAB知识。做完这些自然是开始设计，总体框架，具体步骤，都在脑袋瓜里转了几转。动工，这时才深刻体会到“纸上学来终觉浅，用到实际才叫真！”开始设计从噪声选取，到低通，高通，带阻，带通滤波器的设计，其中是错漏百出。当然，作为老手自然不会再去犯什么程序设计上的低级错误，主要还是对滤波器的认知不够深刻。有的只是照抄书本，最后是连自己都为自己所犯的错误惊愕不已。其中许多错误差点让自己从头来过，一度让自己有放弃的打算。最后还是坚持了下来，自己大叹不容易。作完了，自己心里有些高兴。至于报告，自己还不是手到擒来。可是，到真正开始写报告才发觉，报告才是检验自己掌握程度的真正难关。报告是对理论的理解，这让刚刚大叹理论的我又开始缩手缩脚。看来自己还是老毛病，理论基础不过关。搅尽脑子，搜肠刮肚，抓耳挠勺，终于让自己掌握的有限的理论知识化为了一篇报告。设计结束了，自己也再次用到了时间过得很快的措辞。这回是真的感概，不再是少年不知愁滋味，作为一名大三的学生，另一个代名词便是“准大学毕业生”。一谈到毕业，以前总是报之一笑，然后诙谐的挥挥手，嘴里嘟哝“毕业，早着！”一步步的“早着”，毕业也就不早了。通过此次的课程设计让自己意识到自己已可以做些实际东西，有了基础，但是，自己的一系列的错误，一系列的踌躇也让自己清醒的认识到自己已无时间让自己去挥霍了。初浅的理论知识，不再适合今后的就业。作为专业课的数字信号处理的课程设计，是作为对自己专业能力检测的一个重要手段。认真，细心，专业，理论知识，实际应用这些说了无数遍的名词，真真切切的让自己在这次设计中体会到了。设计的时间过了，但它带来的影响我想是不灭的！要有一个良好的就业前景。拥有扎实的理论知识，认真细致的处事态度，随机应变的实际应用方法，不拘一格的思路都将使得自己能有好的未来。参考资料1 高西全,丁玉美,阔永红.数字信号处理-原理、实现及应用.M.北京:电子工业出版社,20062 张延华,姚林泉,郭玮.数字信号处理-基础与应用. M.北京:机械工业出版社,20053 王宏. MATLAB6.5及其在信号处理中的应用.M. 北京:清华大学出版社,2004附录MATLAB程序：clf;%原始音频抽样Fs=10000;bits=16;T=1;ts=1/Fs;N=38;x,Fs,bits=wavread('ding.wav'); %利用Windows下的录音机，录制一段自己的语音信号，时间为1s x=x(:,1); m=length(x);sound(x,Fs,bits);y=fft(x,m); %对录制的信号进行傅立叶变换f=(Fs/m)*1:m;t=1:m/Fs;figure;subplot(211);plot(t,x);title('原始信号波形');xlabel('time(s)');subplot(212);plot(f,abs(y);title('原始信号频谱');xlabel('frequency(hz)');%噪声信号noise1=0.01*sin(2*pi*6000*t); %高频噪声noise2=0.01*sin(2*pi*200*t); %低频噪声noise3=0.01*sin(2*pi*4000*t)+0.01*sin(2*pi*5000*t); %带阻噪声noise4=0.01*sin(2*pi*5000*t)+0.01*sin(2*pi*200*t); %带通噪声% 低通滤波 fp1=2000; fs1=3000; wp1=2*pi*fp1/Fs; ws1=2*pi*fs1/Fs; X1=x+noise1'%加了高频噪声的信号 sound(X1,Fs,bits);%回放噪声信号 Y11=fft(X1,m);%对加了高频噪声的信号进行傅立叶变换 figure; subplot(321);plot(t,X1);title('加了高频噪声的信号时域图');xlabel('time(s)'); subplot(322);plot(f,abs(Y11);title('加了高频噪声的信号频谱图');xlabel('frequency(hz)'); wc1=(wp1+ws1)/2; b1=fir1(N,wc1,kaiser(N+1);%调用kaiser计算低通FIRDF的b1（n） hn1 w1=freqz(b1,1,512); subplot(323);plot(w1/pi,20*log10(abs(hn1);title('低通滤波器幅频');xlabel('w'); subplot(324);plot(w1/pi,unwrap(angle(hn1);title('低通滤波器相频');xlabel('w'); x1get=filter(b1,1,X1); Y1get=fft(x1get,m); sound(x1get,Fs,bits); subplot(325);plot(t,x1get);title('加了高频噪声的信号滤波后时域图');xlabel('time(s)');subplot(326);plot(f,abs(Y1get);title('加了高频噪声的信号滤波后频谱图');xlabel('frequency(hz)');%高通滤波fp2=500;fs2=300;wp2=2*fp2/Fs;ws2=2*fs2/Fs;X2=x+noise2'%加了低频噪声的信号sound(X2,Fs,bits);%回放噪声信号Y21=fft(X2,m);%对加了低频噪声的信号进行傅立叶变换figure;subplot(321);plot(t,X2);title('加了低频噪声的信号时域图');xlabel('time(s)');subplot(322);plot(f,abs(Y21);title('加了低频噪声的信号频谱图');xlabel('frequency(hz)');wc2=(wp2+ws2)/2;b2=fir1(N,wc2,'high',kaiser(N+1);%调用kaiser计算高通FIRDF的b2（n）hn2 w2=freqz(b2,1,512);subplot(323);plot(w2/pi,20*log10(abs(hn2);title('高通滤波器幅频');xlabel('w');subplot(324);plot(w2/pi,unwrap(angle(hn2);title('高通滤波器相频');xlabel('w');x2get=filter(b2,1,X2);Y2get=fft(x2get,m);sound(x2get,Fs,bits);subplot(325);plot(t,x2get);title('加了低频噪声的信号滤波后时域图');xlabel('time(s)');subplot(326);plot(f,abs(Y2get);title('加了低频噪声的信号滤波后频谱图');xlabel('frequency(hz)');%带阻滤波fpl3=2000;fph3=8000;fsl3=3000;fsh3=7000;wp3=2*fpl3/Fs 2*fph3/Fs;ws3=2*fsl3/Fs 2*fsh3/Fs;X3=x+noise3'%加了带阻噪声的信号sound(X3,Fs,bits);%回放噪声信号Y31=fft(X3,m);%对加了带阻噪声的信号进行傅立叶变换figure;subplot(321);plot(t,X3);title('加了带阻噪声的信号时域图');xlabel('time(s)');subplot(322);plot(f,abs(Y31);title('加了带阻噪声的信号频谱图');xlabel('frequency(hz)');wc3=(wp3+ws3)/2;b3=fir1(N,wc3,'stop',kaiser(N+1);%调用kaiser计算带阻FIRDF的b3（n）hn3 wp3=freqz(b3,1,512);subplot(323);plot(w2/pi,20*log10(abs(hn3);title('带阻滤波器幅频');xlabel('w');subplot(324