目錄
1 clc; 2 clear; 3 close all; 4 5 [x1,fs,bits]=wavread('qq.wav'); %讀取語音信號的數據，賦給變量x1 6 L=length(x1); 7 sound(x1,fs); %播放語音信號 8 9 N=L; 10 k=[0:N-1]; 11 y1=fft(x1,N); %對信號做L個點FFT變換 12 figure(1); 13 subplot(2,1,1); 14 plot(k(1:20000),abs(y1(1:20000))); %做原始語音信號的FFT頻譜圖 15 title('原始語音信號FFT頻譜'); 16 17 f=fs*(0:N-1)/N; 18 subplot(2,1,2); 19 plot(f(1:20000),abs(y1(1:20000))); 20 title('原始語音信號頻譜'); 21 xlabel('f Hz');ylabel('fuzhi');
說明：程序第17行的解釋見“四、DFT的周期性”公式四。
說明：上圖是FFT變換的頻譜圖（數字頻率k），下圖是轉換成模擬頻率下的頻譜圖。
1 %給原始的語音信號加上一個高頻余弦噪聲，頻率為5kHz。畫出加噪后的語音信號時域和頻譜圖，與原始信號對比，可以很明顯的看出區別。 2 t=[0:1/fs:(L-1)/fs]; %將所加噪聲信號的點數調整到與原始信號相同，構造采樣時間點（模擬時間） 3 Au=0.03; %噪聲幅度 4 d=[Au*cos(2*pi*8000*t)]'; %噪聲為8kHz（2*pi*8000/(2*pi)）的余弦信號(模擬時間) 5 x2=x1+d; 6 7 figure(2); 8 subplot(2,1,1); 9 plot(x1); %做原始語音信號的時域圖形 10 title('原始語音信號');xlabel('time n');ylabel('fuzhi'); 11 12 sound(x2,fs); %播放加噪聲后的語音信號 13 subplot(2,1,2); 14 plot(x2); %做原始語音信號的時域圖形 15 title('加噪后的信號');xlabel('time n');ylabel('fuzhi'); 16 17 y2=fft(x2,N); %對信號做L個點FFT變換 18 figure(3); 19 subplot(2,1,1); 20 plot(k(1:20000),abs(y2(1:20000))); %做原始語音信號的FFT頻譜圖 21 title('加躁語音信號FFT頻譜'); 22 23 subplot(2,1,2); 24 plot(f(1:20000),abs(y2(1:20000))); 25 title('加躁語音信號頻譜'); 26 xlabel('f Hz');ylabel('fuzhi');
說明：程序前5行是構造高頻雜波信號，高頻噪聲信號也是按照fs為采樣頻率的，而且采樣點數也是L（與原始信號等長）。
說明：通過頻譜圖就可以看到所加噪聲信號是8Khz。
1 wp=2*pi*4000; %通帶邊界角頻率 2 ws=2*pi*5000; %阻帶邊界角頻率 3 Rp=1; %通帶最大衰減 4 Rs=15; %阻帶最小衰減 5 [NN,Wn]=buttord(wp,ws,Rp,Rs,'s'); %選擇濾波器的最小階數 6 [Z,P,K]=buttap(NN); %創建butterworth模擬濾波器 7 [Bap,Aap]=zp2tf(Z,P,K); 8 [b,a]=lp2lp(Bap,Aap,Wn); 9 [bz,az]=bilinear(b,a,fs); %用雙線性變換法實現模擬濾波器到數字濾波器的轉換 10 11 [H,W]=freqz(bz,az); %繪制頻率響應曲線 12 figure(4); 13 plot(W*fs/(2*pi),abs(H)); 14 grid; 15 xlabel('頻率／Hz'); 16 ylabel('頻率響應幅度'); 17 title('Butterworth');
說明：
（1）前4行是設置波特沃斯濾波器的指標
（2）第5行是根據前4行的參數計算濾波器的階數NN和3db截止頻率Wn。
（3）第6行根據階數NN計算巴特沃斯歸一化模擬低通原型濾波器系統函數的零、極點和增益因子。
（4）第9行將模擬濾波器轉換到數字濾波器，而bz、az分別是分子和分母的系數。
到此，有了bz、az，濾波器就構造出來了。
（5）第11行計算濾波器的頻率響應，freqz函數的使用見freqz()--matlab函數。
（6）第13行的f（模擬頻率）=W*fs/(2*pi)解釋見“2.2.2采樣后的離散傅里葉頻譜”。
說明：可以看到巴特沃斯濾波器的特性是與之前設置的參數匹配的。
x3=filter(bz,az,x2); figure(5); subplot(2,1,1); plot(t,x2); %畫出濾波前的時域圖 title('濾波前的時域波形'); subplot(2,1,2); plot(t,x3); %畫出濾波后的時域圖 title('濾波后的時域波形'); sound(x3,fs); %播放濾波后的信號 y2=fft(x2,N); %對信號做L個點FFT變換 figure(6); subplot(2,1,1); plot(f(1:20000),abs(y2(1:20000))); title('加躁語音信號頻譜'); xlabel('f Hz');ylabel('fuzhi'); y2=fft(x3,N); %對信號做L個點FFT變換 subplot(2,1,2); plot(f(1:20000),abs(y2(1:20000))); title('去噪語音信號頻譜'); xlabel('f Hz');ylabel('fuzhi');
說明：程序的第一行“x3=filter(bz,az,x2);”就是調用上邊構造的巴特沃斯低通濾波器對含躁信號x2濾波，其實就是數值運算。
說明：通過對比可以看到，8Khz的高頻噪聲確實是被濾掉了。
附：C調音符與頻率對照表
西電《數字信號處理》第三版