使用布萊克曼窗設計的FIR濾波器正文

1、音樂信號濾波去噪使用布萊克曼窗設計的 FIR 濾波器摘 要 本課程設計主要是用麥克風采集一段語音信號，繪制其波形并觀察其頻譜。然后 在該語言信號中加一個噪音，利用布萊克曼窗設計一個 FIR 濾波器，對該語音信號進行慮 噪處理，然后比較濾波前后的波形與頻譜。在本課程設計中，是用MATLAB的集成環(huán)境完成一系列的設計。首先對加噪的語音信號進行慮波去噪處理，再比較濾波前后的頻率響應 曲線，若一樣則滿足所設計指標，否則不滿足。也可以調用函數sound 聽濾波前后其語音信號是否帶有噪聲。若無噪聲也說明該濾波器的設置也是成功的。關鍵詞 音樂信號； MATLA；B FIR 濾波器；濾波去噪1 引 言人們在語

2、音通信的過程中將不可避免的會受到來自周圍環(huán)境的干擾，例如傳輸媒介引 入的噪聲，通信設備部的電噪聲，乃至其他講話者的話音等。正因為有這些干擾噪聲的存 在，接受者接受到的語音已不是原始的純凈語音信號，而是受噪聲干擾污染的帶噪聲語音 信號。而本課程設計就是利用 MATLAB!成環(huán)境用布萊克曼窗的方法設計一個 FIR濾波器， 對語音信號進行濾波去噪處理，并將慮噪前后的頻譜圖進行對比。1.1 課程設計目的數字信號處理課程設計是數字信號處理課程的重要實踐性環(huán)節(jié)，是學生在校期間一次 較全面的工程師能力訓練，在實現學生總體培養(yǎng)目標中占有重要地位。綜合運用本課程的 理論知識進行頻譜分析以及濾波器設計，通過理論推

3、導得出相應結論，并利用MATLAB乍為編程工具進行計算機實現，從而復習鞏固了課堂所學的理論知識，提高了對所學知識的 綜合應用能力，并從實踐上初步實現了對數字信號的處理。本課程設計能使學生對通信工 程領域各種技術的DSP實現的設計有較熟練的掌握。且通過自身的實踐，對 DSP的設計程 序、容和方法有更深入的掌握，提高實際運用的能力。并可綜合運用這些知識解決一定的 實際問題，使學生在所學知識的綜合運用能力上以及分析問題、解決問題能力上得到一定 的提高。1.2 課程設計的要求（1）濾波器指標必須符合工程實際。（2）設計完后應檢查其頻率響應曲線是否滿足指標。（3）處理結果和分析結論應該一致，而且應符合理

4、論。（4）獨立完成課程設計并按要求編寫課程設計報告書。1.3 設計平臺MATLA名稱是有兩個英文單詞 Matrix和Laboratory的前三個字母組成。MATLAB7.0 是美國MathWorks公司開發(fā)的優(yōu)秀計算軟件 MATLA的最新版本。MATLA自20世紀80年 代面世以來，以其強大的數值計算能力、優(yōu)秀的繪圖功能以及與其他軟件良好的交互功能 在眾多的數學計算軟件中獨領風騷，特別是它源代碼的開放性使用戶可以二次開發(fā)，受到 了廣大使用者的格外贊賞。MATLA是 一個為科學和工程計算機專門設計的交互式大型軟件，是一個可以完成各種 精確計算和數據處理的、可視化的、強大的計算工具。它集圖和精確計

5、算與一身，在應用 數學、物理、化工、機電工程、醫(yī)藥、金融和其他需要進行復雜數值計算的領域得到了廣 泛應用。它不僅是一個在各類工程設計中便于使用的計算工具，在世界各地的高等院校中 十分流行，在各類工業(yè)應用中更有不俗的表現。MATLA可以幾乎所有的PC機和大型計算機上運行，適用于WindowUNIX等多種系統平臺。本課程設計我們就可以直接誒使用 MATLAB 提供的模塊，實現模擬通信系統的仿真。MATLAB軟件有很強的開放性和適應性。在保持核不變的情況下，MATLAB可以針對不同的應用學科推出相應的工具箱，目前已經推出了圖像處理工具箱、信號處理工具箱、小 波工具箱、神經網絡工具箱以及通信工具箱等多

6、個學科的專用工具箱，極大的方便了不同 學科的研究工作。國已有越來越多的科研和技術人員認識到MABLAB的強大作用，并在不同領域使用MATLA來快速實現科研構想和提高工作效率。2 設計原理2.1 FIR 濾波器FIR（Finite Impulse Response） 濾波器：有限長單位沖激響應濾波器，是數字信號處 理系統中最基本的元件，它可以在保證任意幅頻特性的同時具有嚴格的線性相頻特性，同 時其單位抽樣響應是有限長的，因而濾波器是穩(wěn)定的系統。因此， FIR 濾波器在通信、圖 像處理、模式識別等領域都有著廣泛的應用。FIR 數字濾波器以其良好的線性特性被廣泛應用于現代電子通信系統中，是數字信號

7、處理的重要容之一。在實際信號處理中，往往要求系統兼具實時性和靈活性，而已有的一 些軟件或硬件實現方案（如DSP）則難以同時達到這兩方面的要求。使用具有并行處理特性 的FPGA來實現FIR濾波器，既有很強的實時性，又兼顧了靈活性，為數字信號處理提供 了一種很好的解決方案。 FIR 濾波器系數計算較為繁瑣，在設計時借助 Matlab 工具箱，選 擇合適的窗函數，可以方便地計算濾波器系數，并分析其幅頻、相頻特性。有限長單位沖激響應(FIR)濾波器有以下特點：(1) 、系統的單位沖激響應h (n)在有限個n值處不為零；、系統函數H在|z|>0處收斂，極點全部在z = 0處(因果系統)；(3) 、

8、 結構上主要是非遞歸結構，沒有輸出到輸入的反饋，但有些結構中(例如頻率 抽樣結構)也包含有反饋的遞歸部分。設FIR濾波器的單位沖激響應h (n)為一個N點序列，0 w n < N 1，則濾波器的 系統函數為：N1-nH(z) =h(n)z-n(2-1)n0就是說，它有(N-1)階極點在z = 0處，有(N-1)個零點位于有限z平面的任 何位置因此H(z)是永遠穩(wěn)定的。穩(wěn)定和相位特性是 FIR濾波器突出的優(yōu)點。FIR濾波器有 以下幾種基本結構：橫截型(卷積型、直接型) 、級聯型、頻率抽樣型、快速卷積結構。FIRDF 的設計方法主要分為兩類：第一類是基于逼近理想濾波器特性的方法，包括窗 函數

9、法、頻率采樣法和等波紋最佳逼近法；第二類是最優(yōu)設計法。2.2 窗口設計法數字信號處理的主要數學工具是博里葉變換 而傅里葉變換是研究整個時間域和 頻率域的關系。不過，當運用計算機實現工程測試信號處理時，不可能對無限長的信 號進行測量和運算， 而是取其有限的時間片段進行分析。 做法是從信號中截取一個時 間片段，然后用觀察的信號時間片段進行周期延拓處理，得到虛擬的無限長的信號， 然后就可以對信號進行傅里葉變換、 相關分析等數學處理。 無線長的信號被截斷以后， 其頻譜發(fā)生了畸變， 原來集中在 f(0) 處的能量被分散到兩個較寬的頻帶中去了(這種現象稱之為頻譜能量泄漏 )。為了減少頻譜能量泄漏， 可采用

10、不同的截取函數對信號進行截斷， 截斷函數稱為 窗函數，簡稱為窗。窗函數設計法的基本思路是用 FIRDF逼近希望的濾波特性。設希望逼近的濾波器的頻 率響應函數為hd(ej°)，其單位脈沖響應用hd(n)表示。為了設計簡單方便，通常選擇hd(ej°) 為具有片段常數特性的理想濾波器。 因此 hd(n) 是無限長非因果序列， 不能直接作為 FIRDF 的單位脈沖響應。窗函數設計法就是截取 hd(n)為有限長的一段因果序列，并用合適的窗 函數進行加權做為FIRDF的單位脈沖響應h(n)。下面介紹窗函數設計法的基本設計過程。窗口設計法的主要工作是計算hd(n)和w(n),但當FL(e

11、ju)較為復雜時，hd(n)就不容易 由反付里葉變換求得。這時一般可用離散付里葉變換代替連續(xù)付里葉變換，求得近似值。 窗口法的設計步驟如下：(1) 、通過傅里葉變換忽的理想濾波器的單位脈沖響應hd(n)。(2) 、根據指標選擇窗口形狀、大小和位置。確定窗口類型的主要依據是過渡帶寬和阻帶 最小衰耗的指標。(3) 、給定理想頻響由hd(ejJ和hd(n)，加窗得h(n)=w(n)h d(n)。(4) 、檢驗濾波器的性能。由h (n)求H(ej J是否在誤差容限之。如果不滿足，則返回 第(2)步。以上步驟中hd(n)、H(ej w)的計算可采用傅氏變換的現成公式和程序，窗函數w(n)也是現成的。但整

12、個設計過程不能一次完成，因為窗口類型和大小的選擇沒有解析公式可一次算 出。整個設計可用計算機編程來做。窗口法的優(yōu)點是簡單，有閉合的公式可用，性能及參數都有表格資料可查，計算程序 簡單，較為實用。缺點是當H(ej J較為復雜時，hd(n)就不容易由反付里葉變換求得。邊界頻率因為加窗的影響而不易控制。窗口函數對理想特性的影響：改變了理想頻響的邊沿特性，形成過渡帶，寬為4 N， 等于 WR©)的主瓣寬度；過渡帶兩旁產生肩峰和余振(帶、帶外起伏)，取決于 WR© )的旁瓣，旁瓣多，余振多；旁瓣相對值大，肩峰強，與 N無關；N增加,過渡帶寬減小，肩 峰值不變。因主瓣附近W sin(

13、Nd)“弘 n 沁(2-2)sin( /2)N /2x其中x=Nw/2,所以N的改變不能改變主瓣與旁瓣的比例關系，只能改變 W(3)的絕 對值大小和起伏的密度，當N增加時，幅值變大，頻率軸變密，而最大肩峰永遠為 8.95%， 這種現象稱為吉布斯(Gibbs)效應。肩峰值的大小決定了濾波器通帶的平穩(wěn)程度和阻帶的衰減，所以對濾波器的性能有很 大的影響。改變窗函數的形狀，可改善濾波器的特性，窗函數有許多種，但要滿足以下兩 點要求：窗譜主瓣寬度要窄，以獲得較陡的過渡帶；相對于主瓣幅度，旁瓣要盡可能小， 使能量盡量集中在主瓣中，這樣就可以減小肩峰和余振，以提高阻帶衰減和通帶平穩(wěn)性。但實際上對同樣長度的窗

14、這兩點不能兼得，一般總是通過增加主瓣寬度來換取對旁瓣的抑2.3布萊克曼窗布萊克曼窗的時域形式可表示為:2 n4 n,、c、w(n) 0.42 0.5(cos) 0.08cos()Rn( n)(2-3)N 1N 1它的頻域特性為：W )0.42W r( )0.25W R(-)W R(10.04W r()Wr(N 11)(2-4)其中Wr()為矩形窗函數的幅度頻率特性。增加一個二次諧波余弦分量，可進一步降低旁瓣，但主瓣寬度進一步增加，為 12 N 加N可減少過渡帶。布萊克曼窗函數的最大旁瓣之比主瓣值低 57db,但是主瓣寬度是矩形 窗函數的主瓣寬度的三倍。布萊克曼窗主瓣寬，旁瓣小，頻率識別精度最低

15、，但幅值 識別精度最高3設計步驟3.1設計流程圖本課程設計是對錄制的語音信號進行加噪處理并分析加噪前后語音信號的時域圖與 頻域圖，再用布萊克曼窗設計一個 FIR濾波器，而濾波器的設計必須符合其相應的指標, 否則將不能濾掉加入的噪聲。最后將濾波前后的波形圖進行比較看是否相同。下面是整個 課程設計的流程圖如圖3.1所示：錄制語音信號在語音信號中加入噪聲分析語音信號的頻譜用布萊克曼窗設計FIR濾波器用濾波器對信號進行濾波比較濾波前后的波型及頻回放語音信號圖3.1語音信號的整個流程圖3.2錄制語音信號點擊桌面上的“開始”菜單，再選擇“程序”中的“附件”，在“附件“的菜單欄中選擇“娛樂”，最后點擊“錄音

16、機”。就可以得到如圖3.2的圖。然后點擊開始錄制語音信號，時間大約在23秒之間。圖3.2錄制語音的錄音機錄制好語音信號后，打開MATLA軟件平臺，利用函數waveread對語音信號進行采樣， 記住采樣頻率的和采樣點數。再調用函數sou nd此時可以聽見錄制的語音。采樣完后再語音信號中加入一個單頻噪聲，單頻的噪聲的頻率可以自己設置。按照加入噪聲后的采樣頻 率調用sou nd函數，這時可以明顯的聽見播放的語音信號中有尖銳的單頻嘯叫聲。下面是調用該語言信號以及加入噪聲的程序：x,fs,bits=wavread('e:yuyin.wav');%輸入參數為文件的全路徑和文件名，輸出的第一

17、個參數是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數sou nd(x,fs,bits); %按指定的采樣率和每樣本編碼位數回放N=length(x); %計算信號x的長度fn=2100;% 單頻噪聲頻率，此參數可改t=O:1/fs:(N-1)/fs;%計算時間圍，樣本數除以采樣頻率X=x(:,1)：y=x+0.1*si n(fn*2*pi*t);sound (y,fs,bits); % 明顯聽出有尖銳的單頻嘯叫聲現在是對加入噪聲前后的語音信號進行頻譜分析，先對原始和加噪后的語音信號進行 傅里葉變換，再計算頻譜的頻率圍和譜線間隔。最后就可以畫出未加入噪聲和加入

18、噪聲后 的時域圖和頻域圖。將所有未加和加入噪聲的時域圖和頻域圖畫在同一個圖中，便于比較 和分析。F面是對未加和加入噪聲的頻譜分析的程序: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;%計算頻譜頻率圍原始語音信號語音信號幅度譜頻率(單位：血)加入單頻干擾后的語音信號幅度譜圖3.3語音信號的時域圖和頻率圖3.3濾波器設計而設本課程設計就是要設計一個濾波器慮掉加入的噪聲， 使其恢復原始的語音信號計濾波器的方法

19、有很多，例如：窗函數法、頻率采樣法、脈沖響應不變法和雙線性變換法。 而本課程設計采用的是窗函數法設計FIR濾波器。而FIR濾波器的設計也有很多方法。在Matlab中，可以利用矩形窗、三角窗、漢寧窗、漢明窗、布萊克曼窗、凱塞窗等設計FIR濾波器。而本次采用的是布萊克曼窗來設計濾波器。在用布萊克曼窗設計濾波器的時候，首先要確定濾波器的性能指標。從六種窗函數的 基本參數中我們可以得到旁瓣峰值n=-57,過度帶寬B=11 M ,最小阻帶衰減s=74db,這就表明在設置這些值時其參數必須不大于這些值。而其它帶阻濾波器的設計指標則要根 據加入噪聲的頻率來確定。若不能按照這些來設計濾波器則不可能慮掉噪聲。當

20、所有的指 標都設置完后，可以用這些數字來計算上下邊帶的中心頻率和頻率間隔，并計算布萊克曼 窗設計該濾波器所需要的階數和產生幾階的布萊克曼窗。當所有的準備工作完成后就可以 調用自編的函數計算理想帶阻濾波器的脈沖響應和用窗函數法計算實際的濾波器的脈沖 響應。最后調用freqz函數得到濾波器的頻率特性。從畫出的圖中可以清楚的看見濾波器 的幅頻和相頻特性。下面是用布萊克曼窗設計濾波器的整個程序： fpd=1800;fsd=2050;fsu=1950;fpu=2000;Rp=1;As=70;% 帶阻濾波器設計指標fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=mi n( (fsd-

21、fpd),(fpu-fsu);%計算上下邊帶中心頻率，和頻率間隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; %將 Hz 為單位的模擬頻率換算為rad為單位的數字頻率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(10*pi/dw)+1; %計算布萊克曼窗設計該濾波器時需要的階數n=0:M-1;%定義時間圍w_black=blackman(M); % 產生M階的布萊克曼窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %調用自編函數計算理想帶阻濾波器的脈沖

22、響應h_bs=w_black'.*hd_bs; %用窗口法計算實際濾波器脈沖響應db,mag,pha,grd,w=freqz_m(h_bs,1); %調用自編函數計算濾波器的頻率特性濾波器的北1C|'3000 05 0 10 15 C 2bE濾波器的脈沖響應d訓9濾滾器的相頻特性環(huán) 005115trt0403020 100200040006000圖3.4濾波器的頻率特性3.4信號濾波處理對語音信號信號進行濾波處理主要是濾掉加入的噪聲。不同的濾波器利用不同的函數對語音信號進行濾波，FIR濾波器利用函數fftfilt對信號進行濾波。對信號進行濾波處理要先利用函數filter 對y進

23、行濾波，然后對y進行傅里葉變換' 而畫頻譜時只取前面一半。最后在同一個圖中畫出原始信號的、加入噪聲的語音信號以及 濾波后語音信號的頻域圖和時域圖。這樣便于將所有的圖進行對比和分析，而且還可以直 觀的觀察該課程設計是否成功。當將設計好的濾波器濾掉噪聲后我們也可以再一次調用函 數sou nd,聽此時的聲音是否與原始語音信號基本一樣，若沒有單頻嘯叫聲則說明此次設 計是成功的，否則應重新設置指標。下面是對語音信號進行濾波的程序：y_fil=filter(h_bs,1,y); %用設計好的濾波器對y進行濾波Y_fil=fft(y_fil); %對y進行傅里葉變化圖3.5濾波前后語音信號對比sou

24、nd (x , fs , bits);可以感覺濾波前后的聲音有變化。3.5結果分析要確定本課程設計是否成功就得看原始信號的頻域圖和時域圖與經過濾波器后的語音 信號的頻域圖和時域圖是否一樣，若一樣則表示該設計是成功的，否則是不成功的。在第一個圖中：第一幅圖和第二幅圖是原始語音信號的時域圖和頻域圖，第三幅圖和第四幅圖是加入頻率為 2100的噪聲。從圖中可以看出，第一圖和第三圖相比因為加入噪 聲的緣故所以第三圖y軸的幅度要比第一圖要大，但其形狀還是基本沒有改變。而第二圖 與第四圖相比較在頻率f=2100時多了一個尖銳脈沖。說明原始語音信號加入噪聲是成功 的。在濾波器頻率特性的圖中可以看到：第一個圖是

25、以db為單位的幅頻特性，第二圖是幅頻特性，第三個圖是濾波器的相頻特性，最后一個圖是濾波器的脈沖響應。從圖中可以 清楚的了解濾波器的幅頻和相頻特性。在濾波前后信號比較的圖中我們可以得到：原始的語音信號與濾波后的信號的圖基本 一樣，只是濾波后的圖在原始信號的基礎上有所延遲。所以用布萊克曼窗設計的濾波器是 符合要求的，也就是說該課程設計是成功的。4 出現的問題及解決方法 雖然課程設計已經完成了，但是在設計的過程中還是遇到了許多的問題?？偨Y起來大 概有以下幾個方面：首先，最主要的是要把設計濾波器的參數設置正確才能濾除語音信號中的噪聲。有幾 次因為前面的噪聲頻率設置為 f=2100 ，而后面的帶阻濾波器

26、設計指標沒有與前面的頻率相 匹配導致畫出來的圖怎么也不能濾掉噪音。后來慢慢的改變指標使其在噪聲頻率的左右， 這時才能濾掉噪聲得到原始的語音信號。還有用布萊克曼窗設計的濾波器時期過渡帶寬度 和阻帶最小衰減必須符合布萊克曼窗的基本參數否則也不可能得到所期望的結果。其次，在對加入噪聲后的信號進行濾波時沒有使用正確的形式也不能濾去噪聲。后來 在同學的幫助下解決了這些問題。再次，在利用自編函數計算理想帶阻濾波器的脈沖響應，如果在MATLAB件中的Work下沒有定義，那么不能調用自編的函數，否則將會報錯。最后，在MATLAB件下編程時最好新建一個File文檔。因為在編程的過程過有可能 出現錯誤，如果建一個

27、文檔有助于程序出現錯誤時可以在文檔中直接修改，這樣可以省很 多的時間，又這個課程設計的程序多而繁雜，一不小心就有可能寫錯，如果在工作環(huán)境下 修改這樣利于將所有的程序復制在課程設計中，而且還要對復制后的程序進行刪除。5 結束語為期兩周的數字信號處理課程設計已經結束了， 但在這次設計中我學到了許多的東西。 通過這次的設計，不僅加深了我對課本基礎理論知識的理解，而且增強了我的實踐能力， 同時更加認識到理論知識和實踐結合的重要性。首先，更加深入理解了濾波器設計的各個 關鍵環(huán)節(jié)，包括在什么情況下使用哪種方法設計 FIR 濾波器最好以及在選擇特定的窗函數 進行濾波器的設計時我們應該怎樣確定其性能指標；其次

28、，更加深刻的認識了語音原始信 號與加噪后語音信號的波形及頻譜；再次，較提高了綜合運用專業(yè)基礎知識及軟件設計能 力，在一定程度上對自己的動手能力有很大的幫助。雖然這次課程設計已經完成了，但是遇到的困難也是很多的。其中最主要的問題要屬 怎樣設置濾波器的指標問題，如果指標的設置有問題那么后續(xù)的工作就不可能得到原始的 語音信號。在設置過程中有很多次因為設置的參數不合適而導致設計的濾波器不能慮出單 頻噪聲信號。所以在設計指標問題時一定要結合布萊克曼本身的特點還要考慮加入噪聲的頻率。其次就是一些函數的細節(jié)問題。雖然在這次課程設計中遇到很多的困難，但通過自 己查找有關資料以及老師和同學的幫助下都一一解決了，

29、而且在與同學交流的過程中使同 學之間的感情更進一步。這次設計不僅讓我學會如何獨立完成一項工作，而且提高了獨立 解決問題的能力，為以后的課程設計打下良好的基礎。在此向幫助我的老師及熱心同學表示忠心的感！希望今后還能參加更多的課程設計， 以鍛煉自己在各個方面的能力，尤其是綜合運用專業(yè)基礎知識和實踐結合的能力。設計的 過程中，我通過查閱大量有關資料，與同學交流經驗和自學，并向老師請教等方式，使我 學到了不少的東西，雖然有許多的辛酸，但是看到自己課程設計完成后心中的那份激動是 無法用言語來形容的。參考文獻1 丁玉美，高西全，闊永紅數字信號處理第一版：電子科技,2001年2 后金.數字信號處理.第三版.

30、：高等教育業(yè)，2004年3 程佩青.數字信號處理教程.第四版.:清華大學，2002年4 敏，玲.Matlab通信仿真與應用.第二版：國防工業(yè)，2001年 圣勤.MATLAB7.0實用教程.第三版.：機械工業(yè)，2006附錄：用布萊克曼窗設計FIR濾波器的整個源程序x,fs,bits=wavread('e:yuyin.wav');%輸入參數為文件的全路徑和文件名，輸出的第一個參數是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼 位數sou nd(x,fs,bits); %按指定的采樣率和每樣本編碼位數回放N=length(x); % 計算信號x的長度f

31、n=2100;% 單頻噪聲頻率，此參數可改t=0:1/fs:(N-1)/fs;%計算時間圍，樣本數除以采樣頻率x=x(:,1)'y=x+0.1*si n(fn*2*pi*t);%對原始信號和加噪信號進行fft變換sound (y,fs,bits); %明顯聽出有尖銳的單頻嘯叫聲X=abs(fft(x);Y=abs(fft(y);deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % subplot(2,2,1);plot(t,x);xlabel(' 語音信號 ');計算頻譜頻率圍s)');ylabel('幅度

32、');title('原始時間(單位：subplot(2,2,2);plot(f,X);xlabel('頻率(單位：Hz)');ylabel('幅度譜 ');title('語音信號幅度譜 ');subplot(2,2,3);plot(t,y);xlabel('時間(單位：s)');ylabel('幅度');title('加入單頻干擾后的語音信號 ');subplot(2,2,4);plot(f,Y);xlabel('頻率(單位：Hz)');ylabel('幅度譜

33、');title('加入單頻干擾后的語音信號幅度譜 '); fpd=1800;fsd=2050;fsu=1950;fpu=2000;Rp=1;As=70; % 帶阻濾波器設計指標 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 為單位的數字頻率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(10*pi/dw)+1; % 計算布萊克曼窗設計該濾波器時需要的階數 n=0:M-1; % 定義時間圍w_black=blackman(M); % 產生M階的布萊克曼窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %調用自編函數計算理想帶阻濾波器的脈沖響應h_bs=w_black'.*hd_bs; % 用窗口法計算實際濾波器脈沖響應db,mag,pha,grd,w=freqz_m(h_bs,1); % subplot(2,2,1);plot(w/pi,db);title('