語音信號濾波去噪——使用雙線性變換法設

1、 語音信號濾波去噪使用雙線性變換法設計的切比雪夫I型濾波器 第 頁 共19頁語音信號濾波去噪使用雙線性變換法設計的并聯(lián)型切比雪夫I型濾波器學生姓名：歐曉燕 指導老師：吳志敏摘 要 本課程設計是采用雙線性變換法設計的切比雪夫I型濾波器對雙音頻信號濾波去噪。在網上下載一段雙音頻信號，在MATLAB集成環(huán)境下，首先用wavread函數(shù)求出雙音頻信號的相關參數(shù)，對雙音頻信號進行讀取和加噪；然后再給定相應技術指標，設計一個滿足指標的切比雪夫I型濾波器，對該雙音頻信號進行濾波去噪處理，并繪制對比圖，比較濾波前后的波形和頻譜并進行分析；最后通過回放雙音頻信號，對比濾波前后的信號變換。本課程設計成功的對雙音頻

2、信號進行濾波去噪，初步完成了設計指標。關鍵詞 雙音頻信號； 濾波設計； MATLAB；切比雪夫I型濾波器1 引 言用麥克風采集一段8000Hz，8k的雙音頻信號，繪制波形并觀察其頻譜，給定通帶截止頻率為2000Hz，阻帶截止頻率為2150Hz，通帶波紋為1dB，阻帶波紋為35dB，用雙線性變換法設計的一個滿足上述指標的切比雪夫I型IIR濾波器，對該雙音頻信號進行濾波去噪處理。1.1 課程設計目的數(shù)字信號處理課程設計是在學生完成數(shù)字信號處理和MATLAB的結合后的基本實驗以后開設的。本課程設計的目的是為了讓學生綜合數(shù)字信號處理和MATLAB并實現(xiàn)一個較為完整的小型濾波系統(tǒng)。這一點與驗證性的基本實

3、驗有本質性的區(qū)別。開設課程設計環(huán)節(jié)的主要目的是通過系統(tǒng)設計、軟件仿真、程序安排與調試、寫實習報告等步驟，使學生初步掌握工程設計的具體步驟和方法，提高分析問題和解決問題的能力，提高實際應用水平。1.2 課程設計的要求（1）學會 MATLAB 的使用，掌握 MATLAB 的程序設計方法；（2）濾波器指標必須符合工程實際，根據(jù)模擬濾波器的性能指標，確定數(shù)字濾波器指標；（3）采用雙線性變換法，設計滿足上述性能指標要求的ChebyshevI型數(shù)字低通濾波器；（4）設計完后應檢查其頻率響應曲線是否滿足指標；（5）處理結果和分析結論應該一致，而且應符合理論；（6）獨立完成課程設計并按要求編寫課程設計報告書；

4、1.3 設計平臺本次課程設計是在MATLAB軟件平臺上進行的。MATLAB是矩陣實驗室（MATRIX LABORATORY）的簡稱，是美國MATHWORKS公司推出的具有強大數(shù)值分析、矩陣運算、圖形繪制和數(shù)據(jù)處理等功能的軟件，現(xiàn)已廣泛應用到教學、科研、工程設計等領域2。隨著MATLAB軟件信號處理工具箱的推出，MATLAB已成為信息處理，特別是數(shù)字信號處理DSP應用中分析和設計的主要工具。就MATLAB信號處理中的濾波器設計而言，在很大程度上能快速有效地實現(xiàn)濾波器的分析、設計及仿真，大大節(jié)約了設計時間，相對傳統(tǒng)設計而言，簡化了濾波器設計的難度。2 設計原理用麥克風采集一段雙音頻信號，繪制波形并

5、觀察其頻譜，給定相應技術指標，用雙線性變換法設計的切比雪夫I型IIR濾波器，對該雙音頻信號進行濾波去噪處理，比較濾波前后的波形和頻譜并進行分析。2.1 IIR濾波器從離散時間來看，若系統(tǒng)的單位抽樣(沖激)響應延伸到無窮長，稱之為“無限長單位沖激響應系統(tǒng)”，簡稱為IIR系統(tǒng)。無限長單位沖激響應（IIR）濾波器有以下幾個特點:(1) 系統(tǒng)的單位沖激響應h(n)是無限長；(2) 系統(tǒng)函數(shù)H(z)在有限z平面（0<<）；(3) 結構上存在著輸出到輸入的反饋，也就是結構上是遞歸型的。IIR濾波器采用遞歸型結構，即結構上帶有反饋環(huán)路。同一種系統(tǒng)函數(shù)H（z）可以有多種不同的結構，基本網絡結構有直

6、接型、直接型、級聯(lián)型、并聯(lián)型四種，都具有反饋回路。同時，IIR數(shù)字濾波器在設計上可以借助成熟的模擬濾波器的成果，巴特沃斯(Butterworth)濾波器、切比雪夫(Chebyshev)濾波器、橢圓(Cauer)濾波器、貝塞爾(Bessel)濾波器等，這些典型的濾波器各有特點。有現(xiàn)成的設計數(shù)據(jù)或圖表可查，在設計一個IIR數(shù)字濾波器時，我們根據(jù)指標先寫出模擬濾波器的公式，然后通過一定的變換，將模擬濾波器的公式轉換成數(shù)字濾波器的公式。2.2 切比雪夫I型濾器切比雪夫濾波器（又譯車比雪夫濾波器）是在通帶或阻帶上頻率響應幅度等波紋波動的濾波器。在通帶波動的為“I型切比雪夫濾波器”，在阻帶波動的為“II型

7、切比雪夫濾波器”。切比雪夫濾波器在過渡帶比巴特沃斯濾波器的衰減快，但頻率響應的幅頻特性不如后者平坦。切比雪夫濾波器和理想濾波器的頻率響應曲線之間的誤差最小，但是在通頻帶內存在幅度波動。這種濾波器來自切比雪夫多項式，因此得名，用以紀念俄羅斯數(shù)學家巴夫尼提·列波維其·切比雪夫 巴特沃茲濾波器在通帶內幅度特性是單調下降的，如果階次一定，則在靠近截止 處，幅度下降很多，或者說，為了使通帶內的衰減足夠小，需要的階次N很高，為了克服這一缺點，采用切比雪夫多項式來逼近所希望的 。切比雪夫濾波器的 在通帶范圍內是等幅起伏的，所以在同樣的通常內衰減要求下，其階數(shù)較巴特沃茲濾波器要小。切比雪夫

8、濾波器的振幅平方函數(shù)為 （2-1） 式中   c有效通帶截止頻率    與通帶波紋有關的參量， 大,波紋大 0< <1    VN（x）N階切比雪夫多項式 （2-2）|x|1時,|VN(x)|1 |x|>1時, |x|, VN(x) 切比雪夫濾波器的振幅平方特性如圖所示,通帶內的變化范圍為 1(max) (min)時，|x|>1,隨 ， 0 (迅速趨于零)當 =0時， （2-3） N為偶數(shù)，cos2( )=1，得到min， （2-4）N為奇數(shù)，cos2( ，得到max， （2-5）切比雪夫濾

9、波器的振幅平方特性如圖2.1所示。圖2.1 切比雪夫濾波器的振幅平方特性2.3 雙線性變換法雙線性變換法是使數(shù)字信號濾波器的頻率響應與模擬濾波器的頻率響應相似的一種變換方法。為了客服多值映射這一缺點，我們首先把整個s平面壓縮變換到某一中介的s1平面的一條橫帶里（寬度為,即從到),其次再通過上面討論過的標準變換關系將此橫帶變換到這個z平面上去，這樣就使s平面與z平面式一一對應的關系，消除了多值變換性，也就消除了頻譜混疊現(xiàn)象。將s平面整個軸壓縮變換到s1平面軸上的到一段，可以采用以下變換關系： （2-6） 這樣，變換到，變到可將（6）式寫成 （2-7）解析延拓到整個s平面和s1平面，令，則得 （2

10、-8）再將s1平面通過以下標準變換關系映射到z平面： （2-9）從而得到s平面和z平面的單值映射的關系為 （2-10） （2-11）一般來說，為了使模擬濾波器的某一頻率與數(shù)字濾波器的任一頻率有對應的關系，可以引入待定常數(shù)c，使（6）式和（7）式變換成 （2-12） （2-13）仍將代入（13）式，可得 （2-14） （2-15） （14）式和（15）式是s平面與z平面之間的單值映射關系，這種變換就稱為雙線性變+換。3 設計步驟3.1設計流程圖語音信號濾波去噪使用雙線性變換法設計的切比雪夫I型濾波器的設計流程如圖3.1所示：開始結束回放音頻信號比較濾波前后雙音頻信號的波形及頻譜設計好的濾波器進行

11、濾波處理驗證并進行頻譜分析語音信號的采集（wavread函數(shù)），畫時域圖快速傅里葉變換，并且畫頻譜圖設定濾波器性能指標，通帶截止頻率fb=2000，阻帶截止頻率fc=2150，通帶波紋Rp=1，阻帶波紋As=35雙線性變換法設計切比雪夫I型濾波器圖3.1 雙線性變換法切比雪夫I型濾波器對雙音頻信號去噪流程圖3.2語音信號的采集點擊windows系統(tǒng)桌面的“開始”按鈕，點擊開始菜單欄里的“附件”，選擇“錄音機”選項，點擊錄音機“文件”選項，進入“聲音選定”設置，把屬性一欄設置成“8000Hz，8位，單聲道，7KB/秒”(見圖3.2)。點擊確定，然后開始雙音頻信號的采集，采集時間為23秒左右為最佳

12、。采集的聲音文件以“.wav”格式存儲（見圖3.3）。圖3.2 采集聲音的參數(shù)設置圖3.3 采集聲音3.3語音信號的頻譜分析在MATLAB中編輯m函數(shù)，使用wavread函數(shù)讀取采集的聲音文件（.wav)將它賦值給某一向量，再對其進行采樣,然后使用plot語句畫出相關的頻譜圖形。（1）Wavread函數(shù)調用格式：y,fs,nbits=wavread(file)功能說明：采樣值放在向量x中，fs表示采樣頻率（Hz），bits表示采樣位數(shù)。（2）快速傅里葉變換算法FFT計算DFT的函數(shù)fft,其調用格式如下：Xk=fft(x,n)參數(shù)x為被變換的時域序列向量，N是DFT變換區(qū)間長度，當n大于x的長

13、度時，fft函數(shù)自動在x后面補零。，當n小于xn的長度時，fft函數(shù)計算x的前n個元素，忽略其后面的元素。在本次課程設計中，我們利用fft函數(shù)對雙音頻信號進行快速傅里葉變換，就可以得到信號的頻譜特性。（3）聲音采樣文件讀取的程序（文件名：xiaomiao.wav)；雙音頻信號的提取：x,fs,bits=wavread('D:MATLAB7workxiaomiao.wav'); % 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位

14、數(shù)回放N=length(x); % 計算信號x的長度fn=2200; % 單頻噪聲頻率，此參數(shù)可改t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數(shù)除以采樣頻率x=x(:,1)' y=x+sin(fn*2*pi*t);sound(y,fs,bits); % 應該可以明顯聽出有尖銳的單頻嘯叫聲X=abs(fft(x); Y=abs(fft(y); % 對原始信號和加噪信號進行fft變換,取幅度譜X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率范圍得到

15、的原始語音信號和加上噪音后的語音信號的時域波形和頻譜圖如圖3.4所示。圖3.4 原始雙音頻信號和加噪后信號的時域和頻譜圖3.4濾波器設計 設計指標：通帶截止頻率為2000Hz，阻帶截止頻率為2150Hz，通帶波紋為1dB，阻帶波紋為35dB，用雙線性變換法設計的一個滿足上述指標的切比雪夫I型濾波器雙線性變換法設計切比雪夫I型濾波器fp=fn-200;fc=fn-50; %定義通帶和阻帶截止頻率Rp=1;As=35; % 定義通帶波紋和阻帶衰減wp=fp/fs*2*pi; ws=fc/fs*2*pi; %計算對應的數(shù)字頻率T=1;Fs = 1/T; %定義采樣間隔Omegap=(2/T)*tan

16、(wp/2);Omegas=(2/T)*tan(ws/2); %截止頻率預畸變cs,ds = afd_chb1(Omegap,Omegas,Rp,As); %選擇濾波器最小階數(shù)b,a = bilinear(cs,ds,Fs);C,B,A = dir2cas(b,a) %雙線性變換法實現(xiàn)模擬濾波器到數(shù)字濾波器的轉換db,mag,pha,grd,w=freqz_m(b,a); %繪制數(shù)字濾波器頻率響應幅度圖delta=1,zeros(1,99);ha=filter(b,a,delta);得到切比雪夫濾波器的幅度和相位譜如圖3.5所示。圖3.5 切比雪夫濾波器的幅度和相位譜3.5雙音頻的濾波程序及濾

17、波效果圖濾波程序：y_fil=filter(b,a,y); % 用設計好的濾波器對y進行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半畫出濾波前后雙音頻信號的時域和頻譜圖如圖3.6所示。圖3.6 濾波前后雙音頻信號的時域和頻譜圖3.6結果分析在MATLAB中，經sound函數(shù)，對經過切比雪夫I型濾波器之后的信號進行回放，可以聽出濾波之后的信號比原始信號更清晰一些，清除了環(huán)境噪音。通過以下語句來進行語音信號回放比較：sound（x，fs，bits） %播放原始雙音頻信號sound（y，fs，bits）%播放加噪后的雙音頻信號sound（

18、y_fil，fs，bits） %播放經過濾波處理后的雙音頻信號所得結果證明了切比雪夫I型濾波器去噪設計成功。4 出現(xiàn)的問題及解決方法 在這次的課程設計中，由于理論知識的不踏實以及其他各種原因，我們遇到了不少問題。（1）在進行雙音頻信號提取時，進過多次錄取才得到理想的雙音頻信號，在得到理想的波形時，通過多次嘗試，和查找書籍及同學討論，最后猜得到理想的雙音頻信號的時域圖和頻譜圖（2）在運用Matlab設計濾波器時，當編輯完前面兩條程序時無法放出聲音，后來發(fā)現(xiàn)我們應當把采集的雙音頻信號wav文件放到Matlab的work文件夾中。（3）還要在濾波器性能曲線的As處畫一根豎線，這樣更方便看出結果。（4

19、）所有的時間波形橫坐標都要化為時間，濾波前后頻譜的橫坐標應是頻率，這樣在觀察通帶截止頻率和阻帶截止頻率時更加精確，誤差較小。5 課程設計心得體會采用MATLAB設計濾波器，使原來非常繁瑣復雜的程序設計變成了簡單的函數(shù)調用，為濾波器的設和實現(xiàn)開辟了廣闊的天地，尤其是MATLAB工具箱使各個領域的研究人員可以直觀方便地進行科學研究與工程應用。其中的信號處理工具箱、圖像處理工具箱、小波工具箱等更是為數(shù)字濾波研究的蓬勃發(fā)展提供了可能。MATLAB 信號處理工具箱為濾波器設計及分析提供了非常優(yōu)秀的輔助設計工具, 在設計數(shù)字濾波器時, 善于應用MATLAB進行輔助設計, 能夠大大提高設計效率。兩周的課程設

20、計讓我對MATLAB軟件的使用更加的熟練，對切比雪夫濾波器的濾波原理有了更深刻的認識。6 參考文獻1 程佩青數(shù)字信號處理教程M北京：清華大學出版社，20022 薛年喜主編MATLAB 在數(shù)字信號處理中的應用M北京：清華大學出版社，20023 維納K恩格爾，約翰G普羅克斯.數(shù)字信號處理M.西安交通大學出版社，20024 董長虹等. MATLAB信號處理與應用M.北京：國防工業(yè)出版社，20055 美 M.H.海因斯 著，張建華等 譯.數(shù)字信號處理M.北京：科學出版社，20026 張葛祥，李 娜. MATLAB仿真技術與應用M.北京：清華大學出版社，2007附錄一: 源程序設計清單原始雙音頻信號的采

21、集及加噪：x,fs,bits=wavread('D:MATLAB7workxiaomiao.wav'); % 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放N=length(x); % 計算信號x的長度fn=2200; % 單頻噪聲頻率，此參數(shù)可改t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數(shù)除以采樣頻率x=x(:,1)' y=x+sin(fn*2*pi*t); sound(y,fs,bits);

22、 % 應該可以明顯聽出有尖銳的單頻嘯叫聲X=abs(fft(x); Y=abs(fft(y); % 對原始信號和加噪信號進行fft變換,取幅度譜X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分Warning: Integer operands are required for colon operator when used as index.Warning: Integer operands are required for colon operator when used as index.deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-de

23、ltaf; % 計算頻譜頻率范圍subplot(2,2,1);plot(t,x);xlabel('時間(單位:s)');ylabel('幅度');title('原始雙音頻信號');gridsubplot(2,2,2);plot(f,X);xlabel('頻率(單位:Hz)');ylabel('幅度譜');title('原始雙音頻信號幅度譜圖'); axis(0,6*103,0,6000);gridsubplot(2,2,3);plot(t,y);xlabel('時間(單位:s)');

24、ylabel('幅度');title('加噪后的雙音頻信號');gridsubplot(2,2,4);plot(f,Y);xlabel('頻率(單位:Hz)');ylabel('幅度譜');title('加噪后的雙音頻信號幅度譜圖'); axis(0,6*103,0,6000);grid使用雙線性變換法設計切比雪夫I型濾波器：fp=fn-200;fc=fn-50; %定義通帶和阻帶截止頻率Rp=1;As=35; % 定義通帶波紋和阻帶衰減wp=fp/fs*2*pi; ws=fc/fs*2*pi; %計算對應的數(shù)字頻

25、率T=1;Fs = 1/T; %定義采樣間隔Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2); %截止頻率預畸變cs,ds = afd_chb1(Omegap,Omegas,Rp,As); %選擇濾波器最小階數(shù)Warning: Function call afd_chb1 invokes inexact match d:MATLAB7workAFD_CHB1.M.* 切比雪夫-1 濾波器階次 = 14 Warning: Function call u_chb1ap invokes inexact match d:MATLAB7workU_CHB1AP.

26、M.> In AFD_CHB1 at 29b,a = bilinear(cs,ds,Fs);C,B,A = dir2cas(b,a) %雙線性變換法實現(xiàn)模擬濾波器到數(shù)字濾波器的轉換Warning: Function call dir2cas invokes inexact match d:MATLAB7workDIR2CAS.M.C =3.3307e-016B = 1.0000 -0.0321 0.4021 1.0000 -0.2091 0.0238 1.0000 -0.2206 5.8342 1.0000 -0.3400 0.1487 1.0000 -0.5438 0.7109 1.0

27、000 -2.1153 8.5100 1.0000 -12.5391 364.1406A = 1.0000 -1.9050 0.9776 1.0000 -1.9116 0.9540 1.0000 -1.9147 0.9950 1.0000 -1.9185 0.9754 1.0000 -1.9315 0.9425 1.0000 -1.9354 0.9570 1.0000 -1.9430 0.9444db,mag,pha,grd,w=freqz_m(b,a); %繪制數(shù)字濾波器頻率響應幅度圖Warning: Function call freqz_m invokes inexact match d

28、:MATLAB7workFREQZ_M.M. delta=1,zeros(1,99);ha=filter(b,a,delta);figureSubplot(221);plot(w/pi,db); title('切比雪夫I 濾波器幅度(dB)');xlabel('w');ylabel('dB');axis(0,1,-150,10);gridSubplot(222);plot(w/pi,mag);title('切比雪夫I濾波器幅度響應');xlabel('w');ylabel('幅值 |H|'); ax

29、is(0,1,0,1);gridSubplot(223);plot(w/pi,pha);title('切比雪夫I 濾波器相位響應'); xlabel('w');ylabel('相位（單位： ）');axis(0,1,-4,4);grid對加噪信號進行濾波處理并畫出加噪前后的信號時域和頻譜：y_fil=filter(b,a,y); % 用設計好的濾波器對y進行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半Warning: Integer operands are required for

30、colon operator when used as index.figuresubplot(3,2,1);plot(t,x);xlabel('時間(單位:s)');ylabel('幅度');title('原始雙音頻信號');gridsubplot(3,2,2);plot(f,X);xlabel('頻率(單位:Hz)');ylabel('幅度譜');title('原始雙音頻信號幅度譜圖');gridaxis(0,6*103,0,6000)subplot(3,2,3);plot(t,y);xlabe

31、l('時間(單位:s)');ylabel('幅度');title('加噪后的雙音頻信號');gridsubplot(3,2,4);plot(f,Y);xlabel('頻率(單位:Hz)');ylabel('幅度譜');title('加噪后的雙音頻信號幅度譜圖');gridaxis(0,6*103,0,6000)subplot(3,2,5);plot(t,y_fil);xlabel('時間(單位:s)');ylabel('幅度');title('濾波后的雙音頻信

32、號');gridaxis(0,6,-1,1)subplot(3,2,6);plot(f,Y_fil);xlabel('頻率(單位:Hz)');ylabel('幅度譜');title('濾波后的雙音頻信號幅度譜圖');gridaxis(0,6*103,0,6000)附錄二：afd_chb1函數(shù)function b,a = afd_chb1(Wp,Ws,Rp,As);% 切比雪夫-1型模擬低通濾波器設計% -% b,a = afd_chb1(Wp,Ws,Rp,As);% b = Ha(s) 分子的系數(shù)% a = Ha(s) 分母的系數(shù)% Wp = 以弧度/秒為單位的通帶邊緣頻率; Wp >