彈撥音樂濾波去噪——使用GAUSSWIN設計的FIR濾波器.doc

季政府 彈撥音樂濾波去噪-使用GAUSSWIN設計的FIR濾波器 第19頁 共18頁 彈撥音樂濾波去噪使用GAUSSWIN設計的FIR濾波器學生姓名：季政府 指導老師：胡雙紅摘 要 本課程設計主要內(nèi)容是設計利用窗口設計法選擇GAUSSWIN設計一個FIR濾波器，對一段含噪彈撥音樂信號進行濾波去噪處理并根據(jù)濾波前后的波形和頻譜分析濾波性能。本課程設計仿真平臺為MATLAB7.0，開發(fā)工具是M語言編程。首先從網(wǎng)站上下載一段彈撥樂器演奏音樂，繪制波形并觀察其頻譜特點，然后加入一個帶外單頻噪聲，用高斯窗設計一個滿足指標的FIR濾波器，最后對該音樂信號進行濾波去噪處理，比較濾波前后波形和頻譜并進行分析。由分析結(jié)果可知，濾波器后的彈撥音樂信號與原始信號基本一致，即設計的FIR濾波器能夠去除信號中所加單頻噪聲，達到了設計目的。關(guān)鍵詞 課程設計；濾波去噪；FIR濾波器；GAUSSWIN；MATLAB7.01 引 言本課程設計主要是將一段彈撥音樂信號加入噪聲，然后用某種函數(shù)法設計出的FIR濾波器對加入噪聲后的彈撥音樂信號進行濾波去噪處理，處理時采用的是利用窗口設計法選擇GAUSSWIN設計的FIR濾波器，通過課程設計了解 FIR 濾波器設計的原理和步驟，掌握用 Matlab 語言設計濾波器的方法，觀察音樂信號濾波前后的時域波形的比較，加深對濾波器作用的理解1。1.1課程設計的目的數(shù)字信號處理（Digital Signal Proccessing，簡稱DSP）是一門涉及許多學科而廣泛應用于許多領(lǐng)域的新型學科。20世紀60年代后，隨著計算機和信息技術(shù)的飛速發(fā)展，數(shù)字信號處理技術(shù)應運而生并并得到迅速的發(fā)展。在過去的二十多年時間里，數(shù)字信號處理已經(jīng)在通信領(lǐng)域得到極為廣泛的應用。數(shù)字信號處理是利用計算機或?qū)Ｓ锰幚碓O備，以數(shù)字形式對信號進行采集、變換、濾波、估值、增強、壓縮、識別等處理，以得到符合人們的需要的信號形式。在本次課程設計中，最主要的設計是設計FIR濾波器，F(xiàn)IR濾波器的設計方法主要分為兩類，第一類是基于逼近理想濾波器特性的方法，包括窗函數(shù)法、頻域采樣法和等波紋最佳逼近法；第二類是最優(yōu)設計法。本次的課程設計主要采用的是第一類設計方法，是利用GAUSSWIN法設計FIR濾波器對一段彈撥音樂信號進行濾波去噪，通過這一過程，對濾波前后波形進行對比分析得到結(jié)論。此課程設計比較簡單，主要是將書本中的知識應用到現(xiàn)實中，并且根據(jù)自己對設計題目的理解，運用軟件編寫出程序?qū)崿F(xiàn)這一設計，也是我們對數(shù)字信號處理的原理進行驗證的一個過程。對此，也可以加深我們對所學知識的理解，培養(yǎng)我們的動手能力。1.2課程設計的要求（1）濾波器指標必須符合工程實際。（2）設計完后應檢查其頻率響應曲線是否滿足指標。（3）處理結(jié)果和分析結(jié)論應該一致，而且應符合理論。（4）獨立完成課程設計并按要求編寫課程設計報告書。1.3設計平臺MATLABMATLAB是由美國Math Works公司20世紀80年代中期推出的數(shù)學軟件。MATLAB是“Matric Laboratory”的縮寫，意及“矩陣實驗室”，優(yōu)秀的數(shù)值計算能力和卓越的數(shù)據(jù)可視化能力使其很快在數(shù)學軟件中脫穎而出。Matlab已經(jīng)發(fā)展成為多學科、多種工作平臺的功能強大的大型軟件。在歐美的高校和研究機構(gòu)中，MATLAB是一種非常流行的計算機語言，許多重要的學術(shù)刊物上發(fā)表的論文均是用MATLAB來分析計算以及繪制出各種圖形。MATLAB是一完整的并可擴展的計算機環(huán)境，是一種進行科學和工程計算的交互式程序語言。它的基本數(shù)據(jù)單元是不需要指定維數(shù)的矩陣，它可直接用于表達數(shù)學的算式和技術(shù)概念，而普通的高級語言只能對一個個具體的數(shù)據(jù)單元進行操作。因此，解決同樣的數(shù)值計算問題，使用MATLAB要比使用Basic、Fortran和C語言等提高效率許多倍。許多人贊譽它為萬能的數(shù)學“演算紙”。MATLAB采用開放式的環(huán)境，你可以讀到它的算法，并能改變當前的函數(shù)或增添你自己編寫的函數(shù)2。MATLAB包含的內(nèi)容非常豐富，功能強大，可以概括為以下幾個方面：(1)可以在多種操作系統(tǒng)下運行，如DOS、Windows 95/98/2000/2000/NT、Compaq Alpha、LinuxSun Solaris等。(2)有超過500種的數(shù)學、統(tǒng)計、科學及工程方面的函數(shù)，使用簡單快捷，并且有很強的用戶自定義函數(shù)的能力。(3)有強大的圖形繪制和可視化功能，可以進行視覺數(shù)據(jù)處理和分析，進行圖形、圖像的顯示及編輯，能夠繪制二維、三維圖形，使用戶可以制作高質(zhì)量的圖形，從而寫出圖文并茂的文章。(4)有從外部文件及外部硬件設備讀入數(shù)據(jù)的能力。(5)有豐富的工具箱toolbox。各個領(lǐng)域的專家學者將眾多學科領(lǐng)域中常用的算法編寫為一個個子程序，即m文件，這些m文件包含在一個個工具箱中。其工具箱可以分為兩大類，即功能性工具箱和科學性工具箱。功能性工具箱主要用來擴充MATLAB的符號計算、圖形可視化、建模仿真、文字處理等功能以及與硬件實時交互的功能。學科性工具箱是按學科領(lǐng)域來分類的，如信號處理、控制、通信、神經(jīng)網(wǎng)絡圖像處理、系統(tǒng)辨識、魯棒控制、模糊邏輯、小波等工具箱。MATLAB中的信號處理工具箱內(nèi)容豐富，使用簡便。在數(shù)字信號處理中常用的算法，如FFT，卷積，相關(guān)，濾波器設計，參數(shù)模型等，幾乎都只用一條語句即可以調(diào)用。數(shù)字信號處理所常用的函數(shù)有波形的產(chǎn)生、濾波器的分析和設計、傅里葉變換、Z變換等3。2 設計原理2.1 FIR濾波器 濾波器根據(jù)其沖激響應函數(shù)的時域特性，可分為2種，即無限長沖激響應（IIR）濾波器和有限長沖激響應（FIR）濾波器。FIR和IIR的濾波原理都是進行卷積，就是對數(shù)入信號進行某種計算。FIR用處就在于對數(shù)字信號進行必要的處理，得到所需的輸出信號。FIR系統(tǒng)有自己突出的優(yōu)點：系統(tǒng)總是穩(wěn)定的；易實現(xiàn)線性相位；允許設計多通帶（或多阻帶）濾波器，后兩項都是IIR系統(tǒng)不易實現(xiàn)的。FIR數(shù)字濾波器的設計方法有多種，如窗函數(shù)設計法、頻率采樣法和Chebyshev逼近法等。隨著Matlab軟件尤其是Matlab的信號處理工作箱的不斷完善，不僅數(shù)字濾波器的計算機輔助設計有了可能，而且還可以使設計達到最優(yōu)化。FIR數(shù)字濾波器設計的基本步驟如下： （1）確定技術(shù)指標 在設計一個濾波器之前，必須首先根據(jù)工程實際的需要確定濾波器的技術(shù)指標。在很多實際應用中，數(shù)字濾波器常被用來實現(xiàn)選頻操作。因此，指標的形式一般在頻域中給出幅度和相位響應。幅度指標主要以2種方式給出。第一種是絕對指標。他提供對幅度響應函數(shù)的要求，一般應用于FIR濾波器的設計。第二種指標是相對指標。他以分貝值的形式給出要求。本文中濾波器的設計就以線性相位FIR濾波器的設計為例。（2）逼近 確定了技術(shù)指標后，就可以建立一個目標的數(shù)字濾波器模型（通常采用理想的數(shù)字濾波器模型）。之后，利用數(shù)字濾波器的設計方法（窗函數(shù)法、頻率采樣法等），設計出一個實際濾波器模型來逼近給定的目標。 （3）性能分析和計算機仿真 上兩步的結(jié)果是得到以差分或系統(tǒng)函數(shù)或沖激響應描述的濾波器。根據(jù)這個描述就可以分析其頻率特性和相位特性，以驗證設計結(jié)果是否滿足指標要求；或者利用計算機仿真實現(xiàn)設計的濾波器，再分析濾波結(jié)果來判斷4。2.2窗口設計法窗口法設計的基本想法是要選取某一種合適的理想頻率選擇性濾波器（這種濾波器總是有一個非因果，無限長的脈沖響應），然后將它的脈沖響應截斷（或加窗）以得到一個線性相位和因果的FIR濾波器。因此，這種方法的重點在于選擇某種恰當?shù)拇昂瘮?shù)和一個合適的理想濾波器。現(xiàn)用代表一理想頻率選擇性濾波器，它在整個通帶內(nèi)有單位幅度增益和線性相位特性，而阻帶內(nèi)有零響應5。 用窗口設計法基本步驟如下： （1）構(gòu)造希望逼近的頻率響應函數(shù)。以低通線性相位FIRDF設計為例，一般選擇為線性理想低通濾波器，即 （2-1）其中也稱為截止頻率，稱為樣本延遲。（2） 這個濾波器的脈沖響應應該具有無限長，求出。對進行IFT得到 （2-2）注意，是關(guān)于對稱的，這一點對于線性相位FIR濾波器來說是有用的。（3）為了從得到一個FIRl濾波器必須在兩邊將它截斷。為了得到一個長度為M的因果且線性相位的FIR濾波器，就必須有和 （2-3）這種運算叫做“加窗”，一般來說，可以當做是和某一窗函數(shù)相乘而得到的即 （2-4）式中，稱為窗口函數(shù)，其長度為。窗口法設計的基本思想：對于給定的濾波器技術(shù)要求，選擇濾波器長度M和具有最窄主瓣寬度和盡可能小的旁瓣衰減的某個窗函數(shù)6。常見的窗函數(shù)性能表如下表2-1所示。表2-1 常見的窗函數(shù)性能表名稱濾波器過渡帶寬最小阻帶衰減名稱濾波器過渡帶寬最小阻帶衰減矩形1.8/M21dBPARZENWIN6.6/M56db巴特利特6.1/M25dBFLATTOPWIN19.6/M108db漢寧6.2/M44dBGAUSSWIN5.8/M60db漢明6.6/M51dBBARTHANNWIN3.6/M40db布萊克曼11/M74dBBLACKMANHARRIS16.1/M109dbBOHMANWIN5.8/M51.5dbCHEBWIN15.2/M113dbNUTTALLWIN15.4/M108dbTUKEYWIN2.4/M22db2.3 GAUSSIAN窗高斯窗是一種指數(shù)窗。高斯窗譜無負的旁瓣，第一旁瓣衰減達-55dB。高斯幅度譜的主瓣較寬，故而頻率分辨力低。高斯窗函數(shù)常被用來截短一些非周期信號，如指數(shù)衰減信號等。GAUSSWIN的時域形式可以表示為： (2-5)其中k與的取值范圍必須滿足0N,2以下是畫高斯圖和幅度響應圖的主程序：M=100;n=0:M-1;gauss=gausswin(M);hd_bs=ideal_lp(0.5*pi,M);h_bs=hd_bs.*gauss;db,mag,pha,grd,w=freqz_m(h_bs,1); 高斯窗及其幅度響應圖如圖2-1所示。圖2-1高斯窗及其幅度響應圖3設計步驟3.1設計流程圖濾波器設計流程圖如圖3-1所示。在網(wǎng)上下載一段彈撥音樂加入單頻噪聲對音樂信號進行頻譜分析，畫出時域和頻域波形圖 用GAUSSWIN設計FIR濾波器和畫出其頻譜響應波器用濾波器對信號進行濾波比較濾波前后語音信號的波形及頻譜回放音樂信號開始結(jié)束圖3-1 濾波器設計流程圖3.2 下載并截取彈撥音樂信號從網(wǎng)上下載一段彈撥音樂，從中截取一段格式為.wav的語音信號，時間為45S，并加入噪聲。在MALAB平臺上，觀察原始語音信號與加入噪聲后的時域和頻譜圖。原程序如下： x,fs,bits=wavread(高山流水2.wav); %讀取音樂信號的數(shù)據(jù)sound(x,fs,bits);N=length(x); % 計算信號x的長度fn=5000; % 單頻噪聲頻率t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數(shù)除以采樣頻率x=x(:,1);x=x;y=x+ 0.1*sin(fn*2*pi*t); %加單頻噪聲sound(y,fs,bits); % 應該可以明顯聽出有尖銳的單頻嘯叫聲X=abs(fft(x); Y=abs(fft(y); % 對原始信號和加噪信號進行fft變換X=X(1:length(X)/2); Y=Y(1:length(Y)/2); %截取前半部分deltaf=fs/2/length(X); %計算頻譜的譜線間隔 f=0:deltaf:fs/2-deltaf; %計算頻譜頻率范圍figure (1);原始信號與含噪信號的時域圖和頻域圖如圖3-2所示。圖3-2 原始信號和含噪的信號時域圖和頻譜圖3.3 濾波器設計截取好原始信號，接下來的工作是設計一個GAUSSIAN 窗濾波器，利用公式：Rp=-(min(db(1:1:wp/delta_w+1) （3-1）As=-round(max(db(ws/delta_w+1:1:501) （3-2）分別得到通帶波紋及最小阻帶衰減，并輸入相應濾波器（高斯窗）的頻率參數(shù)，接著用freqz求得濾波器的相應頻率特性。源程序如下：fpd=4700;fsd=4850;fsu=5150;fpu=5300;Rp=1;As=60 %帶阻濾波器設計指標fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); %計算上下帶中心頻率和頻率間隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(10.8*pi/dw)+1; %計算GAUSSWIN窗設計該濾波器時需要的階數(shù)n=0:M-1; %定義時間范圍gauss=gausswin(M);hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 調(diào)用自編函數(shù)計算理想帶阻濾波器的脈沖響應h_bs=gauss.*hd_bs; % 用窗口法計算實際濾波器脈沖響應db,mag,pha,grd,w=freqz_m(h_bs,1); %調(diào)用自編函數(shù)計算濾波器的頻率特性figure (2);濾波器幅度響應圖、相位響應圖、脈沖響應圖如圖3-3所示。圖3-3 濾波器幅度和相位響應圖3.4信號濾波處理設計好濾波器后，我們要對語音信號進行濾波，對比濾波前和濾波后的語音信號。原程序如下：y_fil=fftfilt(h_bs,y); %用設計好的濾波器對y進行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:length(Y_fil)/2); %計算頻譜取前一半y_fil=fftfilt(h_bs,y); %用設計好的濾波器對y進行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:length(Y_fil)/2); %計算頻譜取前一半figure (3);彈撥音樂原始信號、加噪后信號與濾波后信號的時域、頻域圖如圖3-4所示。圖3-4濾波前和濾波后的對比3.5結(jié)果分析通過觀察上圖，濾波后的語音信號發(fā)生了衰減，說明濾波器起到了濾波作用，同時通過頻譜對比，可以看出濾波器濾掉了一部分頻率范圍內(nèi)的信號。分別聽原始語音和濾波后的語音信號，發(fā)現(xiàn)濾波后的語音信號噪聲減小了，同時原始信號強度稍有減弱，基本達到了濾波的效果。4出現(xiàn)的問題及解決方法在本次課程設計中我們設計的濾波器對音樂信號的濾噪處理中，我遇到了許多細節(jié)疏忽和一些比較難處理的問題，當遇到問題時，我們不但運用自己學到的知識將其解決，還得到了同學老師的幫助，將其不斷完善和修改。設計中出現(xiàn)的問題及解決方法如下：1、 在調(diào)用音樂時，用了語句x,fs,bits=wavread(高山流水2.wav);出現(xiàn)錯誤，后來得知我下載的高山流水高山流水音樂是MP3格式，不是wav格式，后來用千千靜聽軟件將其改為wav格式之后就正確了。2、 在運用Matlab設計濾波器時，當編輯完前面兩條程序時無法發(fā)出聲音，后來發(fā)現(xiàn)我們應當把音樂信號wav文件放到Matlab的work文件夾中。3、 在對加入噪聲后的音樂信號用FIR濾波器進行去噪處理時，所編寫的函數(shù)出現(xiàn)錯誤，繪制不出圖形，分析錯誤主要是因為調(diào)用FIR濾波器是所用的函數(shù)改為fftfilt函數(shù)。4、 在進行音樂信號的提取時沒設采樣頻率，結(jié)果頻率沒達到要求，通過向老師詢問，在windows自帶錄音機軟件中改變音樂信號頻率的屬性，使時域圖達到了要求。5、 繪制出濾波后的的波形，發(fā)現(xiàn)FIR濾波器沒有濾掉單頻噪聲。通過自己仔細的檢查，是單頻噪聲的頻率改動后，F(xiàn)IR濾波器沒有改動。所以單頻噪聲的頻率也應該先定義，F(xiàn)IR濾波器的截止頻率應該以單頻信號的頻率為中心，這樣重新運行后，結(jié)果正確。6、 在采用 gausswin 窗設計的 FIR 濾波器時得不到理想的濾波器，把通帶截止頻率與阻帶起始頻率之間的差值設置的太小或者太大，后來通過適當?shù)倪x擇參數(shù) fp 和 fs，繪制出來的圖形效果比較明顯，基本符合設計指標。這樣通過 MATLAB 運算出來的濾波器的阻帶波紋達到了濾波器設 計的要求，得到了比較理想的濾波器。7、 還要在濾波器性能曲線的wc處畫一根豎線，這樣更方便看出結(jié)果，其中wc處線的確定還需計算出wb/pi的值。8、 所有的時間波形橫坐標都要化為時間，濾波前后頻譜的橫坐標應是頻率，這樣在觀察同代截止頻率和阻帶截止頻率時更加精確，誤差較小。對利用 gausswin函數(shù)設計 FIR 濾波器的設計步驟很生疏,不知如何計算 M 的值，自己上網(wǎng)查閱資料，或向圖書館借閱資料，掌握利用 gausswin 窗函數(shù)設計 FIR 濾 波器的方法和步驟。通過同學的幫助，利用 matlab 編寫關(guān)于 gausswin 的 M 文件得到 text5 函數(shù)，再通過 4 個不同參數(shù)代入 text5 函數(shù)，計算出過渡帶寬。9、 在對濾波器進行設計的過程中，由于程序較多，總會或多或少出現(xiàn)些錯誤，但是通過請教同學，在程序包中找到了freqz_m.m和ideal_lp.m這兩個程序包，從而解決了對濾波器的設計問題。10、 最后當觀察濾波前和濾波后的時域和頻域波形對比圖時，發(fā)現(xiàn)濾波后的信號頻譜仍有一小段噪聲沒有濾除，但是與未濾波前的噪聲相比，這一小段噪聲可以忽略，通過聽取前后兩段語音信號也沒有出現(xiàn)任何問題。5 結(jié)束語至此，本次課程設計就告一個段落了。這次課程設計中，是對上學期的專業(yè)課數(shù)字信號處理和MATLAB的一次復習，和通過實踐加深了對所學知識的理解，所以這次課程設計對每位同學來說都是受益頗多的。學習的過程是相互討論共同進步的，多多討論課題中遇到的問題，可以鞏固我們的知識掌握能力，增加熟練運用度。從本次課程設計的中心來看，課題是希望將數(shù)字信號處理技術(shù)應用于某一實際領(lǐng)域，這里就是指對音樂的處理。作為存儲于計算機中的音樂信號，其本身就是離散化了的向量，我們只需將這些離散的量提取出來，就可以對其進行處理了。在這里，用到了處理數(shù)字信號的強有力工具MATLAB，通過MATLAB里幾個命令函數(shù)的調(diào)用，很輕易的在實際化音樂與數(shù)字信號的理論之間搭了一座橋。課題的特色在于它將音樂看作了一個向量，于是音樂數(shù)字化了，則可以完全利用數(shù)字信號處理的知識來解決。我們可以像給一般信號做頻譜分析一樣，來給音樂信號做頻譜分析，也可以較容易的用數(shù)字濾波器來對音樂進行濾波處理。改變參數(shù)，理論結(jié)合實際，分析各參數(shù)對圖形的影響，從而加深對各個參數(shù)的理解。在完成這次課程設計過程中學到了許多東西，進一步理解了濾波器設計方法和各參數(shù)意義，通過分析信號時域和頻域的關(guān)系等，加深了對濾波性能的理解，而且學會了使用Matlab一些基本函數(shù)，增加了進一步學習Matlab軟件的興趣。同時，通過本次課程設計，鍛煉了我的動手能力，和提高了我分析問題，解決問題的能力。同學們在一起思考問題，通過親自動手做一個FIR濾波器，發(fā)現(xiàn)了許多潛在的問題，而這些問題是在平常被我們所忽略的，甚至認為是不成問題的。讓我再次感受到我們應當把所學知識和實踐相結(jié)合，才能夠超越水平，取得更大的發(fā)展。在這次課程設計中，胡雙紅老師給了我們很大的幫助和啟發(fā)，隨時隨地為我們答疑，引導我們思考，而不是一味的接受。在這里向所有所有老師和同學給予的熱心支持和幫助表示感謝，謝謝你們。參考文獻1 程佩青數(shù)字信號處理教程M北京：清華大學出版社，20022 薛年喜主編MATLAB 在數(shù)字信號處理中的應用北京：清華大學出版社3 維納K恩格爾，約翰G普羅克斯.數(shù)字信號處理.西安交通大學出版社，20024 陳其宗.數(shù)字信號處理頻譜計算與濾波器設計M.北京：電子工業(yè)出版社，20025 劉敏，魏玲Matlab通信仿真與應用M北京：國防工業(yè)出版社，20016 百科ROBOT,zivenwong.窗函數(shù).百度百科：附錄1：彈撥音樂信號濾波去噪設計源程序清單% 程序名稱：voicenoisefilter.m% 程序功能：采用基于gausswin的窗口設計法，設計FIR濾波器對含噪語音進行濾波去噪處理。% 程序作者：季政府/% 最后修改日期：2012-3-2x,fs,bits=wavread(高山流水2.wav);sound(x,fs,bits);N=length(x); % 計算信號x的長度fn=5000; % 單頻噪聲頻率t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數(shù)除以采樣頻率x=x(:,1);x=x;y=x+ 0.1*sin(fn*2*pi*t); %加單頻噪聲sound(y,fs,bits); % 應該可以明顯聽出有尖銳的單頻嘯叫聲X=abs(fft(x); Y=abs(fft(y); % 對原始信號和加噪信號進行fft變換X=X(1:length(X)/2); Y=Y(1:length(Y)/2); %截取前半部分deltaf=fs/2/length(X); %計算頻譜的譜線間隔 f=0:deltaf:fs/2-deltaf; %計算頻譜頻率范圍figure (1);subplot(2,2,1);plot(t,x);grid on; xlabel(時間（單位：s）);ylabel(幅值);title( 原始語音信號);axis tight;subplot(2,2,2);plot(f,X);grid on; xlabel(頻率（單位：Hz）);ylabel(幅度譜);title(語音信號幅度譜圖);axis tight;subplot(2,2,3);plot(t,y);grid on; xlabel(時間（單位：s）);ylabel(幅值);title( 加入單頻干擾后語音信號);axis tight;subplot(2,2,4);plot(f,Y);grid on; xlabel(頻率（單位：Hz）);ylabel(幅度譜);title(加入干擾后的語音信號幅度譜圖);axis(-inf inf 0 250);fpd=4700;fsd=4850;fsu=5150;fpu=5300;Rp=1;As=60; %帶阻濾波器設計指標fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); %計算上下帶中心頻率和頻率間隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(10.8*pi/dw)+1; %計算GAUSSWIN窗設計該濾波器時需要的階數(shù)n=0:M-1; %定義時間范圍gauss=gausswin(M);hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 調(diào)用自編函數(shù)計算理想帶阻濾波器的脈沖響應h_bs=gauss.*hd_bs; % 用窗口法計算實際濾波器脈沖響應db,mag,pha,grd,w=freqz_m(h_bs,1); %調(diào)用自編函數(shù)計算濾波器的頻率特性figure (2);subplot(2,2,1);plot(w/pi,db,LineWidth,2);axis(-inf,inf,-80,10);a=0 3.15;b=-60 -60;line(a, b,Color,r,LineStyle,-, LineWidth,2);a=0 3.15;b=-1 -1;line(a, b,Color,r,LineStyle,-, LineWidth,2);a=wsd/pi wsd/pi;b=-150 10;line(a, b,Color,r,LineStyle,-, LineWidth,2);a=wsu/pi wsu/pi;b=-150 10;line(a, b,Color,r,LineStyle,-, LineWidth,2);xlabel(w/pi);ylabel(dB);title(濾波