数字信号处理课程设计报告

频率

1.5 2 2.5x104D=2MATLAB编写程序如下，并绘出二倍抽取后的时域波形和频谱：%对原始信号分别进行D倍的直接抽取x1=x(1:2:end);X1=fft(x1);N1=length(x1);n=0:N1-1;f1=n*fs/N1;figure(2);subplot(2,1,1)plot(x1);title('二倍抽取后信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)stem(f1,abs(X1),'.');title('二倍抽取后的信号频谱幅值');xlabel('频率f');ylabel('幅值A');sound(x1,fs/2);%回放抽取后的语音信号二倍抽取后信号时域波形10.5A值 0幅-0.5-10 1000

3000

4000 5000 时间轴n

7000 8000 9000二倍抽取后的信号频谱幅值200150A100幅5000 0.5

频率

1.5 2 2.5x104这里采用直接抽取的方式，而不选择MATLAB的decimate函数。因为decimate函数会先对序列进行抗混叠滤波再抽取，就失去了后续自己设计模拟低通滤波器的意义，所以这里采用直接抽取，后续再自行滤除。分析：可以看到整数倍抽取后原序列的频谱带宽扩展，已抽样序列x(n)和抽取序列y(n)的频谱差别在频率尺度上不同。信号的整数2倍抽取，每个1个抽取一个，降低了采样率，会引起频谱混叠现象，这样就无法恢复原信号，所以后续要进行抗混叠滤波。还是可以清晰地听到“away个单词。3）按内插因子I=2进行内插，将采样率提高2倍利用MATLAB编写程序如下，并绘出二倍抽取后二倍插值的时域波形和频谱：%对信号进行L倍直接插值x2=zeros(1,2*length(x));x2(1:2:end)=x;X2=fft(x2);N2=length(x2);n=0:N2-1;f2=n*fs/N2;figure(3);subplot(2,1,1)plot(x2);title('二倍插值后信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)stem(f2,abs(X2),'.');title('二倍插值后的信号频谱幅值');xlabel('频率f');ylabel('幅值A');sound(x2,fs);二倍插值后信号时域波形10.5A值 0幅-0.5-10 0.5 1 1.5 时间轴n

2.5 3 3.5x104二倍插值后的信号频谱幅值400300A200幅10000 0.5

频率

1.5 2 2.5x104同抽取一样，这里采用直接内插的方式，而不选择MATLAB的interp函数。因为interp函数会先对序列进行滤波再内插。这里采用直接内插，后续再自行滤除。分析：整数2倍内插，在相邻两个样点之间等间隔插入1个0值点，然后进行低通滤波，即可求得2倍内插的结果。往往内插结果，以一定间隔对信号进行采样然后再进行傅里叶变换，频谱出现的变化很明显可以看到镜像频谱。回放语音信号，感觉语音变得很清晰，听起来跟之前抽取后的信号相比似乎更尖锐些，高频部分更清晰，听得更清楚。4）设计汉宁窗低通滤波器恢复出语音信号编写MATLAB代码，绘出滤波器时域和频域波形wp=0.75*pi;ws=0.125*pi;rp=0.25;as=50;Bt=wp-ws;N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;hn=fir1(N-1,wc,'low',hanning(N));%绘出滤波器时域和频域波形滤波器单位脉冲响应0.6

0.40.20-0.21 2 3 4

6 7 8 9 10 11i滤波器频率响应波形)H02

0-50-1000 0.1 0.2 0.3 0.4 0.5/

0.6 0.7 0.8 0.9 1%对抽取插值后的输出信号进行低通滤波h=conv(hn,x2);H=fft(h)figure(5)N3=length(h);n=0:N3-1;f3=n*fs/N3;%绘出模拟滤波后信号时域波形和频谱幅值sound(h,fs*2)%回放语音信号0.5A值 幅

模拟滤波后信号时域波形-0.50 0.5 1 1.5 时间轴n

2.5 3 3.5x104模拟低通滤波后的信号频谱幅值原始信号时域波形400原始信号时域波形30A值.0幅值100幅-050 0.5 1 1.5 2 2.5-1400

2000

频率f x1046000 8000 10000 12000 14000 16000 18000时间轴n原始信号频谱幅值300A200幅10000 0.5

频率

1.5 2 2.5x104分析：与原始信号频谱幅值相对比，可以看出滤波后频谱接近。通过低通滤波，滤除了镜像频谱，使得抽取之后的抽样率仍满足抽样定理要求时，才不会产生混叠失真。频谱与原始信号频谱相近。回放语音信号，可以明显感觉到恢复的语音听起来很自然，声音接近原始语音信号。设计内容读入一段语音信号（或音乐信号）在语音信号中分别加入以下几种噪声：（1）（2）单频噪色（正弦干扰（3）多频噪声（多正弦干扰；进行对比，分析并体会含噪语音信号频域和时域波形的改变3）IIRFIR用所设计的滤波器对被污染的语音信号进行滤波；的信号进行对比，分析信号的变化；设计思想信号加噪给信号分别加入白噪声、单频噪声和多频噪声。这里选用添加高斯白噪声，MATLABwgn()awgn()噪声则分别通过设定不同频率的正弦函数，添加到原信号即可，这里为了后续更方便滤除噪声，正弦噪声的频率设置比较高。噪声滤除滤波设计原理：滤波器的作用是利用离散时间系统的特性，对输入信号波形(或频谱)进行加工与处理，或者说利用数字方法，按预定的要求对信号进行变换。可以将滤波器理解为一个算法或者程序，将代表输入信号的数字时间序列，转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。从性能上说，IIR滤波器传输函数的极点可位于单位圆内的任何地方，因此滤波器却可以得到严格的线性相位输出，但由于FIRIIR滤波器高5~10倍，但是成本较高，信号延时也较大；所以如若按相同的选择性和相同的线性要求来说的话，则IIR加滤波器的节数和复杂性。整体来看，IIR与FIR各有优缺点。IIR滤波器能达到同样效果，且阶数少，延迟也较低，但是会有稳定性和非线性相位问题；FIR滤波器虽然没有稳定性问题，且线性相位，但阶数多，延迟大。IIR滤波器进行设计。这里我根据之前绘出的信号（加噪）频域波形特性，分别对白噪声和正弦噪声设计了两个低通滤波器，滤波器的参数设置如下：%设计巴特沃斯低通滤波器1fp1=1000; 通带截止频fs1=1600; 阻带截止频rp=0.1; 通带最大衰减rs=40; 带最小衰减%设计巴特沃斯低通滤波器2Fp2=1000; 通带截止频Fs2=1500; 阻带截止频率rp=0.1;rs=40;

带最小衰减wp=2*Fp2/fs;ws=2*Fs2/fs; 计算数字滤波器设计指标设计步骤1) 首先读入一段语音信号“hello.wvMATLAB域波形和频谱。原始信号时域波形10.5A值 0幅-0.5-10 5000 10000 15000时间轴n原始信号频谱幅值400300A200幅10000 2000 4000 6000频率f/Hz

8000 100002）MATLAB（下只列出部分关键代码）%向原始语音信号中添加高斯白噪声snr=30;x1=awgn(x,snr);X1=fft(x1,N);N1=length(x1);n=0:N1-1;f1=n*fs/N1;%绘制原始信号和叠加高斯白噪声信号后的时域频域波形sound(x1,fs,nbits);%回放语音信号10.5A值 0幅-0.5-10

原始信号时域波形A值幅5000 10000 15000时间轴n

40030020010000

原始信号频谱幅值5000 10000频率f/Hz叠加高斯白噪声信号时域波形10.5A值 0幅-0.5-1

叠加高斯白噪声信号频谱幅值400300A值200幅10000 5000 10000 时间轴n

0 5000 10000频率f/Hz分析：播放语音信号，可以明显听到语音中夹杂有“嘶嘶”嘈杂的噪声。幅值一定程度增加，主要是高频段频谱变得模糊。3）在语音信号中加入单频噪声，编写MATLAB代码，绘出叠加单频噪声的时域信号和频谱，并与原始信号波形进行了对比。（下只列出部分关键代码）%向原始语音信号中添加单频噪声信号x2=x+dt';X2=fft(x2,N);n=0:N2-1;f2=n*fs/N2;%绘制原始信号和叠加单频噪声信号后的时域频域波形sound(x2,fs,nbits);%回放语音信号10.5A值 幅-0.5-1

原始信号时域波形A值幅0 5000 10000 15000时间轴n

40030020010000

原始信号频谱幅值5000 10000频率f/Hz10.5A值 幅-0.5-1

单频噪声信号波形

叠加单频噪声信号频谱幅值800600A值400幅20000 5000 10000 时间轴n

0 5000 10000频率f/Hz1500Hz1500Hz处出现了一条明显的峰谱线，可推定为尖锐正弦噪声的谱线。4）MATLAB（下只列出部分关键代码）%向原始语音信号中添加多频噪声信号dt1=0.2*sin(2*pi*1600*(1:size(x))/fs);dt2=0.2*sin(2*pi*1700*(1:size(x))/fs);x3=x+dt'+dt1'+dt2';X3=fft(x3,N);N3=length(x3);n=0:N3-1;f3=n*fs/N3;%绘制原始信号和叠加多频噪声信号后的时域频域波形sound(x3,fs,nbits);原始信号时域波形 原始信号频谱幅值4001300A值0A值200幅幅-11000 5000 10000 1500000 5000 10000时间轴n频率f/Hz多频噪声信号波形叠加多频噪声信号频谱幅值15001A值0A值1000幅幅500-10 5000 10000 1500000 5000 10000时间轴n频率f/Hz1500Hz、1600Hz和1700Hz，可以明显与原语音区分开来，也方便后续的滤除。播放语音信号，听1500Hz1600Hz1700Hz噪声谱线。IIRFIR在设计思想中，提到了滤波设计的原理。从性能上说，IIRFIRIIR5~10倍，但是成本较高，信号延时也较大；这里由于对阶数和相位没有特殊的要求，通过观察信号频谱的特性，我在这里选用IIR滤波器，分别对白噪声和正弦噪声设计了两个巴特沃斯低通滤波器。关键部分代码如下：%设计巴特沃斯低通滤波器1fp1=1000;fs1=1600;rp=0.1; 通带最大衰减rs=40; 带最小衰减wp=2*fp1/fs;ws=2*fs1/fs;%计算数字滤波器设计指标[N2,wc]=buttord(wp,ws,rp,rs);[b,a]=butter(N2,wc,'low');%设计巴特沃斯低通滤波器2Fp2=1000; 通带截止频Fs2=1500; 阻带截止频率rp=0.1;rs=40;

带最小衰减wp=2*Fp2/fs;ws=2*Fs2/fs; [N1,wc]=buttord(wp,ws,rp,rs);[b,a]=butter(N1,wc,'low');%用巴特沃斯低通滤波器滤除单频噪声%绘出两个低通滤波器的频谱巴特沃斯低通滤波器频域响应0d（-50度幅-1000 0.1

0.3

0.5 0.6 /

0.8 0.9 1巴特沃斯低通滤波器频域响应0d（-50度幅-1000 0.2 0.4 0.6 0.8 1/

1.2 1.4 1.6 1.8 2用所设计的滤波器对被污染的语音信号进行滤波，编写MATLAB的变化，并回放语音信号，比较变化。%设计巴特沃斯低通滤波器1%用巴特沃斯低通滤波器滤除高斯白噪声z1=filter(b,a,x1);%滤除白噪声Nz1=length(z1);Z1=fft(z1,Nz1);%设计巴特沃斯低通滤波器2%用巴特沃斯低通滤波器滤除单频噪声z2=filter(b,a,x2);Nz2=length(z2); 算滤波后的长度Z2=fft(z2,Nz2)DFT变换%滤除多频噪声z3=filter(b,a,x3);Nz3=length(z3); 算滤波后的长度Z3=fft(z3,Nz3)DFT变换1并与叠加噪声后和叠加噪声前的波形进行对比：10.5A值 幅-110.5A值 幅-1

原始信号时域波形0 5000 10000 时间轴n加高斯白噪声信号的时域波形

400300A值幅1000400300A值幅1000

原始信号频谱幅值0 2000 4000 6000 8000频率f/Hz加高斯白噪声信号的频谱0 5000 10000 时间轴n加高斯白噪声信号滤波后的波形10.5A值 0幅-0.5-1

0 2000 4000 6000 800010000频率f/Hz加高斯白噪声信号滤波后的频谱400300A值200幅10000 5000 时间轴n

15000

0 2000 4000 6000 800010000频率f/Hz观察分析：使用巴特沃斯低通滤波器保留了信号的低频部分，从时域波形上可看出与叠加噪声后相比，波形波动减小，接近原始信号波形，从频谱上看，频率2000Hz以上幅度较小，但由于高斯白噪声在整个信号周期内都存在，故用低通滤波器没有办法完全将其滤除，播放语音时还是会听到一些轻微的噪声，总体噪声减小。利用巴特沃斯低通滤波器对叠加单频噪声信号进行滤波，画出利用巴特沃2噪声前的波形进行对比：原始信号时域波形1

400

原始信号频谱幅值0.5A值 0幅-0.550001000050001000015000200040006000 800010000010.5A值 0幅-0.550001000050001000015000200040006000 8000100000

时间轴n加正弦噪声信号的波形

300A200幅10000800600A400幅20000

频率f/Hz加正弦噪声信号的频谱时间轴n加正弦噪声信号低通滤波后的波形频率f/Hz加正弦噪声语音信号滤波后的频谱14000.5300A值0A值200幅幅-0.51005000100001500050001000015000200040006000 8000100000时间轴n

00频率f/Hz播放语音时没有了原先尖锐的噪声，恢复得到原语音非常接近。利用巴特沃斯低通滤波器对叠加多频噪声信号进行滤波，画出利用巴特沃2噪声前的波形进行对比：原始信号时域波形1A值0幅-1

400300A值幅100

原始信号频谱幅值0 5000 10000时间轴n加多弦噪声信号的波形1A值0

15000

00 2000 4000 6000 800010000频率f/Hz加正弦噪声信号的频谱1500A 1000值幅 -10 5000 10000 15000

50000

2000 4000 6000 800010000时间轴n

0时间轴n加多弦噪声信号低通滤波后的波形频率f/Hz时间轴n加多弦噪声信号低通滤波后的波形频率f/Hz加多弦噪声语音信号滤波后的频谱4001300A值0A值200幅幅-11000 5000 10000 15000频率f/Hz2除单频噪声后语音基本一致。三、设计过程中遇到的问题和解决方案概念不清晰或淡忘。对于信号的整数周期截取、DFT等概念，一开始做课程设计的时候对它们都有一点模糊，通过上网查找资料，询问老师和温顾书本的知识以后，才记起具体的概念，对这些概念有了进一步的理解以后我才开始进行设计。包括多采样率语音信号处理和信号加噪及噪声滤除之前课堂上都没有讲到过，我们都是通过自己查找资料，结合之前所学的滤波器设计才有所掌握，受益匪浅。MATLAB封装decimateinterp函数，但是得到内插和抽取的波形后，可以明显的看到信号已进行了滤波（抗混叠滤波和镜像滤波。这样后续自己设计滤波器MATLAB就相当于把那两个函数重写了一次，可以观察到明显的不同。在对信号进行加噪声时，一开始，我写好了代码，运行也没有问题，即可得到我想要的结果。滤波器的参数设置问题。这个问题让我头疼了一段时间，因为它不是对频域波形的观察，不断调整设置参数，才得到了大致想要的结果。WienerWienerMMSELsA-MMSE、先验信噪比、听觉掩蔽效应等不同的方法可以较好地滤除高斯白噪声，可是最后在设计检查的时候，很遗憾，我最后也没有能用其他有效地方法更好地滤除白噪声。但是这个问题，还可以进一步研究，不是设计结束了就停止了，我在后来的时间里还是通过资料学到了和维纳滤波器滤除高斯白噪声的方法，套用网上原有的代码，也得到了自己的结果。四、设计心得与体会这一次的课程设计为期总的算起来时间不算很长，但确实让我感觉到自己学习了很多的知识。数字信号处理这门课是上学期学的，但是课堂上我们知识学习了一些理论的知识，和题目的计算，对实际操作还是很不了解，通过上学期的两次实验我感觉到学了些有用的东西，通过自己绘出信号波形，自己设计滤波器，感受到了这门学科的价值。但是这一次的课程设计，就不光光是停留在那个层面了。虽然说总体感觉题目也并没有很难，但是我觉得远比在书本上学到的知识有用，有趣。就再通过自己设计一个滤波器来把噪声滤除，这种感觉让人十分地有成就感，也让我觉得特别的有趣，似乎稍微将这门学科实际化了一步，贴近实际生活中的应用，这对我们学习这门学科也有很大的帮助。DSP会有收获，有收获，那实践就是有价值的。五、参考文献报：自然科学版，2012.8六、程序源代码清单DFT和IDFT函数的MATLAB代码如下：%DFT函数function[Xk]=dft(xn,N)%Xk0<=n<=N-1DFT系数数组%xn=N点有限持续时间序列%N=DFT的长度k=[0:1:N-1];WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk;%计算逆离散傅里叶变换function[xn]=idft(Xk,N)%xn=0<=n<=N-1间的N点有限时间序列%Xk=0<=k<=N-1DFT系数数组%N=DFT的长度n=0:1:N-1的行向量k=0:1:N-1的行向量WN=exp(-j*2*pi/N)因子nk=n'*k; nkNNWNnk=WN.^(-nk);矩阵xn=(Xk*WNnk)/N;的行向量%计算N=32时的幅频特性N=32;n=0:1:N-1;xn=cos(7/16*pi.*n)+0.5*cos(9/16*pi.*n)+0.75*cos(0.5*pi.*n);subplot(2,1,1);stem(n,xn,'.');title('信号x(n)');xlabel('n');ylabel('x(n)');axis([0,32,-4,4])X=dft(xn,N)magX=abs(X(1:1:(N/2+1)));k=0:1:(N/2);w=2*pi/N*k;subplot(2,1,2);stem(w/pi,magX,'.');xlabel('\omega/\pi');ylabel('幅度');title('N=32点DFT幅频特性图')%N=30的信号的幅频特性N=30;n=0:1:N-1;xn=cos(7/16*pi.*n)+0.5*cos(9/16*pi.*n)+0.75*cos(0.5*pi.*n);subplot(2,1,1);stem(n,xn,'.');title('信号x(n)');xlabel('n');ylabel('x(n)');axis([0,32,-4,4])X=dft(xn,N)magX=abs(X(1:1:(N/2+1)));k=0:1:(N/2);w=2*pi/N*k;subplot(2,1,2);stem(w/pi,magX,'.');xlabel('\omega/\pi');ylabel('幅度');title('N=30点DFT幅频特性图')%N=40的信号的幅频特性图：N=40;n=0:1:N-1;xn=cos(7/16*pi.*n)+0.5*cos(9/16*pi.*n)+0.75*cos(0.5*pi.*n);subplot(2,1,1);stem(n,xn,'.');title('信号x(n)');xlabel('n');ylabel('x(n)');axis([0,32,-4,4])X=dft(xn,N);magX=abs(X(1:1:(N/2+1)));k=0:1:(N/2);w=2*pi/N*k;subplot(2,1,2);stem(w/pi,magX,'.');xlabel('\omega/\pi');ylabel('幅度');title('N=40点DFT幅频特性图');%多采样率语音信号处理%读取语音信号幅值给y。fs为采样频率[x,fs,bits]=wavread('C:\Users\Administrator\Desktop\away.wav');sound(x,fs);X=fft(x);N=length(x);n=0:N-1;f=n*fs/N;figure(1);subplot(2,1,1)plot(x);title('原始信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)stem(f,abs(X),'.');title('原始信号频谱幅值');xlabel('频率f');ylabel('幅值A');%对原始信号分别进行D倍的直接抽取x1=x(1:2:end);X1=fft(x1);N1=length(x1);n=0:N1-1;f1=n*fs/N1;figure(2);subplot(2,1,1)plot(x1);title('二倍抽取后信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)stem(f1,abs(X1),'.');title('二倍抽取后的信号频谱幅值');xlabel('频率f');ylabel('幅值A');sound(x1,fs/2);%回放抽取后的语音信号%对信号进行L倍直接插值x2=zeros(1,2*length(x));x2(1:2:end)=x;X2=fft(x2);N2=length(x2);n=0:N2-1;f2=n*fs/N2;figure(3);subplot(2,1,1)plot(x2);title('二倍插值后信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)stem(f2,abs(X2),'.');title('二倍插值后的信号频谱幅值');xlabel('频率f');ylabel('幅值A');sound(x2,fs);%设计汉宁窗，低通滤波器wp=0.75*pi;ws=0.125*pi;rp=0.25;as=50;Bt=wp-ws;N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;hn=fir1(N-1,wc,'low',hanning(N));figure(4)subplot(2,1,1);stem(hn,'.');title('滤波器单位脉冲响应');xlabel('i');ylabel('h(i)');w=(0:1023)*2/1024;subplot(2,1,2);滤波器频率响应波形');xlabel('\omega/\pi');ylabel('20lg|Hg(\omega)|');axis([0,1,-100,20]);gridon%低通滤波h=conv(hn,x2);H=fft(h)figure(5)N3=length(h);n=0:N3-1;f3=n*fs/N3;subplot(2,1,1)plot(h);title('模拟滤波后信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)stem(f3,abs(H),'.');title('模拟低通滤波后的信号频谱幅值');xlabel('频率f');ylabel('幅值A');sound(h,fs*2)%第三问[x,fs,nbits]=wavread('C:\Users\Administrator\Desktop\hello.wav');sound(x,fs,nbits);%播放语音信号%原始语音信号N=length(x);n=0:N-1;f=n*fs/N;figure(1);subplot(2,1,1)plot(x);grid;title('原始信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(2,1,2)X=fft(x,N);stem(f,abs(X),'.');grid;title('原始信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([0112500400])%向原始语音信号中添加高斯白噪声snr=30;x1=awgn(x,snr);X1=fft(x1,N);N1=length(x1);n=0:N1-1;f1=n*fs/N1;figure(2);subplot(2,2,1)plot(x);grid;title('原始信号时域波形');xlabel('时间轴n');ylabel('幅值A');axis([015000-1.21.2])subplot(2,2,2)X=fft(x,N);stem(f,abs(X),'.');grid;title('原始信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([0112500400])subplot(2,2,3);时间轴n');ylabel('幅值A');axis([015000-1.21.2])subplot(2,2,4);stem(f1,abs(X1),'.');grid;title('叠加高斯白噪声信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([0112500400])sound(x1,fs,nbits);%向原始语音信号中添加单频噪声信号x2=x+dt';X2=fft(x2,N);n=0:N2-1;f2=n*fs/N2;figure(3)subplot(2,2,1)plot(x);grid;title('原始信号时域波形');xlabel('时间轴n');ylabel('幅值A');axis([015000-1.21.2])subplot(2,2,2)X=fft(x,N);stem(f,abs(X),'.');grid;title('原始信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([0112500400])subplot(2,2,3)plot(x2);grid;title('时间轴幅值A');axis([015000-1.21.2])subplot(2,2,4)频率f/Hz');ylabel('幅值A');axis([0112500800])sound(x2,fs,nbits);%向原始语音信号中添加多频噪声信号x3=x+dt'+dt1'+dt2';X3=fft(x3,N);N3=length(x3);n=0:N3-1;f3=n*fs/N3;figure(4)subplot(2,2,1)plot(x);grid;title('原始信号时域波形');xlabel('时间轴n');ylabel('幅值A');axis([015000-1.51.5])subplot(2,2,2)X=fft(x,N);stem(f,abs(X),'.');grid;title('原始信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([0112500400])subplot(2,2,3)plot(x3);grid;title('时间轴幅值A');axis([015000-1.51.5])subplot(2,2,4)stem(f3,abs(X3),'.');grid;title('叠加多频噪声信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([01125001500])sound(x3,fs,nbits);%设计巴特沃斯低通滤波器1fp1=1000;fs1=1600;rp=0.1; 通带最大衰减rs=40; 带最小衰减wp=2*fp1/fs;ws=2*fs1/fs;%计算数字滤波器设计指标[N2,wc]=buttord(wp,ws,rp,rs);[b,a]=butter(N2,wc,'low');w=0:0.01*pi:pi;[h2,w]=freqz(b,a,w);h2=20*log10(abs(h2));%用巴特沃斯带通滤波器滤除高斯白噪声z1=filter(b,a,x1);%滤除白噪声Nz1=length(z1);Z1=fft(z1,Nz1);%设计巴特沃斯低通滤波2Fp2=1000; 通带截止频Fs2=1500; 阻带截止频率rp=0.1;rs=40;

带最小衰减wp=2*Fp2/fs;ws=2*Fs2/fs; [N1,wc]=buttord(wp,ws,rp,rs);[b,a]=butter(N1,wc,'low');w=0:0.01*pi:pi;[h1,w]=freqz(b,a,w);%计算高通滤波器频率响应h1=20*log10(abs(h1));%用巴特沃斯高通滤波器滤除单频噪声z2=filter(b,a,x2);Nz2=length(z2); 算滤波后的长度Z2=fft(z2,Nz2);%对滤除噪声后的信号进行DFT变换%滤除多频噪声z3=filter(b,a,x3);Nz3=length(z3); 算滤波后的长度Z3=fft(z3,Nz3);%对滤除噪声后的信号进行DFT变换%两个低通滤波器的频谱figure(5)subplot(211)plot(w,-abs(h1),'b');axis([0,1,-120,20])xlabel('\omega/\pi');ylabel('幅度（dB)')title('巴特沃斯低通滤波器1频域响应')subplot(212)plot(w,-abs(h2),'b');axis([0,2,-120,20])xlabel('\omega/\pi');ylabel('幅度（dB)')title('巴特沃斯低通滤波器2频域响应')figure(6)subplot(3,2,1)plot(x);grid;title('原始信号时域波形');xlabel('时间轴n');ylabel('幅值A');subplot(3,2,2)X=fft(x,N);stem(f,abs(X),'.');grid;title('原始信号频谱幅值');xlabel('频率f/Hz');ylabel('幅值A');axis([0112500400])subplot(3,2,3);plot(x2);grid;title('加正弦噪声信号的波形');xlabel('时间轴n');ylabel('