數(shù)字信號處理課程設計報告

頻率

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函數(shù)。因為decimate函數(shù)會先對序列進行抗混疊濾波再抽取，就失去了后續(xù)自己設計模擬低通濾波器的意義，所以這里采用直接抽取，后續(xù)再自行濾除。分析：可以看到整數(shù)倍抽取后原序列的頻譜帶寬擴展，已抽樣序列x(n)和抽取序列y(n)的頻譜差別在頻率尺度上不同。信號的整數(shù)2倍抽取，每個1個抽取一個，降低了采樣率，會引起頻譜混疊現(xiàn)象，這樣就無法恢復原信號，所以后續(xù)要進行抗混疊濾波。還是可以清晰地聽到“away個單詞。3）按內(nèi)插因子I=2進行內(nèi)插，將采樣率提高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同抽取一樣，這里采用直接內(nèi)插的方式，而不選擇MATLAB的interp函數(shù)。因為interp函數(shù)會先對序列進行濾波再內(nèi)插。這里采用直接內(nèi)插，后續(xù)再自行濾除。分析：整數(shù)2倍內(nèi)插，在相鄰兩個樣點之間等間隔插入1個0值點，然后進行低通濾波，即可求得2倍內(nèi)插的結果。往往內(nèi)插結果，以一定間隔對信號進行采樣然后再進行傅里葉變換，頻譜出現(xiàn)的變化很明顯可以看到鏡像頻譜?；胤耪Z音信號，感覺語音變得很清晰，聽起來跟之前抽取后的信號相比似乎更尖銳些，高頻部分更清晰，聽得更清楚。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分析：與原始信號頻譜幅值相對比，可以看出濾波后頻譜接近。通過低通濾波，濾除了鏡像頻譜，使得抽取之后的抽樣率仍滿足抽樣定理要求時，才不會產(chǎn)生混疊失真。頻譜與原始信號頻譜相近。回放語音信號，可以明顯感覺到恢復的語音聽起來很自然，聲音接近原始語音信號。設計內(nèi)容讀入一段語音信號（或音樂信號）在語音信號中分別加入以下幾種噪聲：（1）（2）單頻噪色（正弦干擾（3）多頻噪聲（多正弦干擾；進行對比，分析并體會含噪語音信號頻域和時域波形的改變3）IIRFIR用所設計的濾波器對被污染的語音信號進行濾波；的信號進行對比，分析信號的變化；設計思想信號加噪給信號分別加入白噪聲、單頻噪聲和多頻噪聲。這里選用添加高斯白噪聲，MATLABwgn()awgn()噪聲則分別通過設定不同頻率的正弦函數(shù)，添加到原信號即可，這里為了后續(xù)更方便濾除噪聲，正弦噪聲的頻率設置比較高。噪聲濾除濾波設計原理：濾波器的作用是利用離散時間系統(tǒng)的特性，對輸入信號波形(或頻譜)進行加工與處理，或者說利用數(shù)字方法，按預定的要求對信號進行變換。可以將濾波器理解為一個算法或者程序，將代表輸入信號的數(shù)字時間序列，轉(zhuǎn)化為代表輸出信號的數(shù)字時間序列，并在轉(zhuǎn)化過程中，使信號按預定的形式變化。從性能上說，IIR濾波器傳輸函數(shù)的極點可位于單位圓內(nèi)的任何地方，因此濾波器卻可以得到嚴格的線性相位輸出，但由于FIRIIR濾波器高5~10倍，但是成本較高，信號延時也較大；所以如若按相同的選擇性和相同的線性要求來說的話，則IIR加濾波器的節(jié)數(shù)和復雜性。整體來看，IIR與FIR各有優(yōu)缺點。IIR濾波器能達到同樣效果，且階數(shù)少，延遲也較低，但是會有穩(wěn)定性和非線性相位問題；FIR濾波器雖然沒有穩(wěn)定性問題，且線性相位，但階數(shù)多，延遲大。IIR濾波器進行設計。這里我根據(jù)之前繪出的信號（加噪）頻域波形特性，分別對白噪聲和正弦噪聲設計了兩個低通濾波器，濾波器的參數(shù)設置如下：%設計巴特沃斯低通濾波器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; 計算數(shù)字濾波器設計指標設計步驟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處出現(xiàn)了一條明顯的峰譜線，可推定為尖銳正弦噪聲的譜線。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，可以明顯與原語音區(qū)分開來，也方便后續(xù)的濾除。播放語音信號，聽1500Hz1600Hz1700Hz噪聲譜線。IIRFIR在設計思想中，提到了濾波設計的原理。從性能上說，IIRFIRIIR5~10倍，但是成本較高，信號延時也較大；這里由于對階數(shù)和相位沒有特殊的要求，通過觀察信號頻譜的特性，我在這里選用IIR濾波器，分別對白噪聲和正弦噪聲設計了兩個巴特沃斯低通濾波器。關鍵部分代碼如下：%設計巴特沃斯低通濾波器1fp1=1000;fs1=1600;rp=0.1; 通帶最大衰減rs=40; 帶最小衰減wp=2*fp1/fs;ws=2*fs1/fs;%計算數(shù)字濾波器設計指標[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以上幅度較小，但由于高斯白噪聲在整個信號周期內(nèi)都存在，故用低通濾波器沒有辦法完全將其濾除，播放語音時還是會聽到一些輕微的噪聲，總體噪聲減小。利用巴特沃斯低通濾波器對疊加單頻噪聲信號進行濾波，畫出利用巴特沃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除單頻噪聲后語音基本一致。三、設計過程中遇到的問題和解決方案概念不清晰或淡忘。對于信號的整數(shù)周期截取、DFT等概念，一開始做課程設計的時候?qū)λ鼈兌加幸稽c模糊，通過上網(wǎng)查找資料，詢問老師和溫顧書本的知識以后，才記起具體的概念，對這些概念有了進一步的理解以后我才開始進行設計。包括多采樣率語音信號處理和信號加噪及噪聲濾除之前課堂上都沒有講到過，我們都是通過自己查找資料，結合之前所學的濾波器設計才有所掌握，受益匪淺。MATLAB封裝decimateinterp函數(shù)，但是得到內(nèi)插和抽取的波形后，可以明顯的看到信號已進行了濾波（抗混疊濾波和鏡像濾波。這樣后續(xù)自己設計濾波器MATLAB就相當于把那兩個函數(shù)重寫了一次，可以觀察到明顯的不同。在對信號進行加噪聲時，一開始，我寫好了代碼，運行也沒有問題，即可得到我想要的結果。濾波器的參數(shù)設置問題。這個問題讓我頭疼了一段時間，因為它不是對頻域波形的觀察，不斷調(diào)整設置參數(shù)，才得到了大致想要的結果。WienerWienerMMSELsA-MMSE、先驗信噪比、聽覺掩蔽效應等不同的方法可以較好地濾除高斯白噪聲，可是最后在設計檢查的時候，很遺憾，我最后也沒有能用其他有效地方法更好地濾除白噪聲。但是這個問題，還可以進一步研究，不是設計結束了就停止了，我在后來的時間里還是通過資料學到了和維納濾波器濾除高斯白噪聲的方法，套用網(wǎng)上原有的代碼，也得到了自己的結果。四、設計心得與體會這一次的課程設計為期總的算起來時間不算很長，但確實讓我感覺到自己學習了很多的知識。數(shù)字信號處理這門課是上學期學的，但是課堂上我們知識學習了一些理論的知識，和題目的計算，對實際操作還是很不了解，通過上學期的兩次實驗我感覺到學了些有用的東西，通過自己繪出信號波形，自己設計濾波器，感受到了這門學科的價值。但是這一次的課程設計，就不光光是停留在那個層面了。雖然說總體感覺題目也并沒有很難，但是我覺得遠比在書本上學到的知識有用，有趣。就再通過自己設計一個濾波器來把噪聲濾除，這種感覺讓人十分地有成就感，也讓我覺得特別的有趣，似乎稍微將這門學科實際化了一步，貼近實際生活中的應用，這對我們學習這門學科也有很大的幫助。DSP會有收獲，有收獲，那實踐就是有價值的。五、參考文獻報：自然科學版，2012.8六、程序源代碼清單DFT和IDFT函數(shù)的MATLAB代碼如下：%DFT函數(shù)function[Xk]=dft(xn,N)%Xk0<=n<=N-1DFT系數(shù)數(shù)組%xn=N點有限持續(xù)時間序列%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系數(shù)數(shù)組%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;%計算數(shù)字濾波器設計指標[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('