吹管音樂濾波去噪基于漢寧創(chuàng)的FIR濾波器

1、 . 吹管音樂濾波去噪基于漢寧窗的FIR濾波器學(xué)生： 指導(dǎo)老師：摘 要 從上下載一段吹管樂器演奏音樂，利用CE軟件對音樂進行編輯。繪制波形并觀察其頻譜特點，加入一個帶外單頻噪聲，用漢寧窗設(shè)計一個滿足指標(biāo)的FIR濾波器，對該含噪音樂信號進行濾波去噪處理，比較濾波前后波形和頻譜并進行分析，根據(jù)結(jié)果和學(xué)過的理論得出合理結(jié)論。與不同信源一樣濾波方法的同學(xué)比較各種信源的特點，與一樣信源不同濾波方法的同學(xué)比較各種濾波方法性能優(yōu)劣。關(guān)鍵詞 濾波去噪；FIR濾波器；漢寧窗；MATLAB1 引 言本課程設(shè)計主要針對一段吹管音樂信號，在進行加噪后，利用窗函數(shù)設(shè)計法選擇漢寧窗設(shè)計的FIR濾波器，對加噪后的吹管音樂信

2、號進行濾波去噪處理，并對前后時域波形和頻域波形進行對比分析的程序設(shè)計。1.1 課程設(shè)計目的本次課設(shè)中的主要目的是讓學(xué)生在熟悉Matlab語言環(huán)境，掌握其語言編程規(guī)則的前提下，利用漢寧窗設(shè)計一個符合要求的FIR濾波器來實現(xiàn)音樂信號的濾波去噪，并繪制濾波前后的時域波形和頻譜圖。根據(jù)圖形分析判斷濾波器設(shè)計的正確性。通過本次課設(shè)，我們能夠?qū)W會如何綜合運用課堂上學(xué)會的理論知識，增強自己的動能力與聯(lián)系實際的能力，為以后的工作奠定基礎(chǔ)。1.2 課程設(shè)計的要求（1）濾波器指標(biāo)必須符合工程實際。（2）設(shè)計完后應(yīng)檢查其頻率響應(yīng)曲線是否滿足指標(biāo)。（3）處理結(jié)果和分析結(jié)論應(yīng)該一致，而且應(yīng)符合理論。（4）獨立完成課程設(shè)

3、計并按要求編寫課程設(shè)計報告書。1.3 課程設(shè)計平臺課程設(shè)計的主要設(shè)計平臺是MATLAB 7.0。MATLAB 的名稱源自 Matrix Laboratory ，它是美國MathWorks公司生產(chǎn)的一個為科學(xué)和工程計算專門設(shè)計的交互式大型軟件，是一個可以完成各種精確計算和數(shù)據(jù)處理的、可視化的、強大的計算工具。它集圖示和精確計算于一身，在應(yīng)用數(shù)學(xué)、物理、化工、機電工程、醫(yī)藥、金融和其他需要進行復(fù)雜數(shù)值計算的領(lǐng)域得到廣泛應(yīng)用。它不僅是一個在各類工程設(shè)計中便于使用的計算工具，而且也是一個在數(shù)學(xué)、數(shù)值分析和工程計算等課程教學(xué)中的優(yōu)秀的教學(xué)工具，在世界各地的高等院校中十分流行，在各類工業(yè)應(yīng)用中更有不俗的表

4、現(xiàn)。MATLAB可以在幾乎所有的PC機和大型計算機上運行，適用于Windows、UNIX等各種系統(tǒng)平臺。1MATLAB軟件包括五大通用功能：數(shù)值計算功能（Nemeric）；符號運算功能（Symbolic）；數(shù)據(jù)可視化功能（Graphic）；數(shù)據(jù)圖形文字統(tǒng)一處理功能（Notebook）和建模仿真可視化功能（Simulink）。其中，符號運算功能的實現(xiàn)是通過請求MAPLE核計算并將結(jié)果返回到MATLAB命令窗口。該軟件有三大特點：一是功能強大；二是界面友善、語言自然；三是開放性強。2MATLAB在信號與系統(tǒng)中的應(yīng)用主要包括符號運算和數(shù)值計算仿真分析。由于信號與系統(tǒng)課程的許多容都是基于公式演算，而M

5、ATLAB借助符號數(shù)學(xué)工具箱提供的符號運算功能能基本滿足信號與系統(tǒng)課程的需求。例如，解微分方程、傅里葉正反變換、拉普拉斯正反變換、z正反變換等。MATLAB在信號與系統(tǒng)中的另一主要應(yīng)用是數(shù)值計算與仿真分析，主要包括函數(shù)波形繪制、函數(shù)運算、沖激響應(yīng)與階躍響應(yīng)仿真分析、信號的時域分析、信號的頻譜分析、系統(tǒng)的S域分析、零極點圖繪制等容。數(shù)值計算仿真分析可以幫助學(xué)生更深入理解信號與系統(tǒng)的理論知識，并為將來使用MATLAB進行信號處理領(lǐng)域的各種分析和實際應(yīng)用打下基礎(chǔ)3。2 設(shè)計原理2.1 數(shù)字信號處理數(shù)字信號處理（Digital Signal Processing，簡稱DSP）是一門涉與許多學(xué)科而又廣泛

6、應(yīng)用于許多領(lǐng)域的新興學(xué)科4。20世紀(jì)60年代以來，隨著計算機和信息技術(shù)的飛速發(fā)展，數(shù)字信號處理技術(shù)應(yīng)運而生并得到迅速的發(fā)展。在過去的二十多年時間里，數(shù)字信號處理已經(jīng)在通信等領(lǐng)域得到極為廣泛的應(yīng)用。數(shù)字信號處理是利用計算機或?qū)Ｓ锰幚碓O(shè)備，以數(shù)字形式對信號進行采集、變換、濾波、估值、增強、壓縮、識別等處理，以得到符合人們需要的信號形式。52.2 FIR濾波器有限長單位脈沖響應(yīng)數(shù)字濾波器（Finite Impulse Response Digital Filter,縮寫FIRDF）簡稱FIR濾波器，是數(shù)字信號處理系統(tǒng)中最基本的原件，其最大優(yōu)點是可以實現(xiàn)線性相位濾波，可以在保證任意幅頻特性的同時具有嚴(yán)

7、格的線性相頻特性，滿足了在數(shù)字通信和圖像傳輸與處理等應(yīng)用場合對線性相位的要求。FIR濾波器是全零點濾波器，硬件和軟件實現(xiàn)結(jié)構(gòu)簡單，因而是十分穩(wěn)定的系統(tǒng)。6FIR濾波器的設(shè)計方法主要分為兩類：第一類是基于逼近理想濾波器器特性的方法包括窗函數(shù)法、頻率采樣法、和等波紋最佳逼近法；第二類是最優(yōu)設(shè)計法。本次課設(shè)采用的是第一類設(shè)計法中的窗函數(shù)法。設(shè)FIR濾波器的單位脈沖響應(yīng)的長度為，則其頻率響應(yīng)函數(shù)為（2-1）一般將表示成如下形式： （2-2）式中，是的實函數(shù)（可以去負值）。與前面的表示形式，即相比，與不同。與不同。為了區(qū)別于幅頻響應(yīng)函數(shù)和相頻響應(yīng)函數(shù)，稱為幅頻特性函數(shù)，稱為相頻特性函數(shù)。第一類線性相位F

8、IR濾波器的相位特性函數(shù)是的嚴(yán)格線性函數(shù)：（2-3）2.3 窗口設(shè)計法窗口設(shè)計法是一種通過截斷和計權(quán)的方法使無限長非因果序列成為有限長脈沖響應(yīng)序列的設(shè)計方法。通常在設(shè)計濾波器之前，應(yīng)該先根據(jù)具體的工程應(yīng)用確定濾波器的技術(shù)指標(biāo)。在大多數(shù)實際應(yīng)用中，數(shù)字濾波器常常被用來實現(xiàn)選頻操作，所以指標(biāo)的形式一般為在頻域中以分貝值給出的相對幅度響應(yīng)和相位響應(yīng)。6窗口設(shè)計法基本步驟如下：（1）根據(jù)過渡帶寬與阻帶衰減要求，選擇窗函數(shù)的類型并估計窗口長度N。窗函數(shù)的類型可根據(jù)最小阻帶衰減AS獨立選擇。（2）根據(jù)待求濾波器的理想頻率響應(yīng)求出理想單位脈沖響應(yīng)hd(n)。（3）由性能指標(biāo)確定窗函數(shù)W(n)和長度N。（4）

9、求得實際濾波器的單位脈沖響應(yīng)h(n)，h(n)即為所設(shè)計FIR濾波器系數(shù)向量b(n)。常見的窗函數(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/M113d

10、BNUTTALLWIN15.4/M108dBTUKEYWIN2.4/M22dB2.4 漢寧窗（Hanning window）漢寧窗函數(shù)是余弦平方函數(shù)，又稱之為升余弦函數(shù)，它的時域形式可以表為：（2-8）其中k=1,2,，k。它的頻域幅度特性函數(shù)為：（2-9）其中為矩形窗函數(shù)的幅度頻率特性函數(shù)。漢寧窗函數(shù)的最大旁瓣值比主瓣值低31dB，但是主瓣寬度比矩形窗函數(shù)的主瓣寬度增加了一倍，為。漢寧窗函數(shù)的時域幅度與頻域幅度特性曲線的MATLAB實現(xiàn)的曲線圖如圖2-1所示。圖2-1 漢寧窗函數(shù)的時域幅度與頻域幅度特性曲線3設(shè)計步驟3.1 設(shè)計流程圖本課程設(shè)計主要是從上下載一段吹管樂器演奏音樂，利用CE軟件

11、對音樂進行編輯。繪制波形并觀察其頻譜特點，加入一個帶外單頻噪聲，用漢寧窗設(shè)計一個滿足指標(biāo)的FIR濾波器，對該含噪音樂信號進行濾波去噪處理，比較濾波前后波形和頻譜并進行分析，根據(jù)結(jié)果和學(xué)過的理論得出合理結(jié)論。程序的設(shè)計流程圖如下圖3-1所示。 開始下載一段吹管音樂信號，用CE軟件編輯格式為wav.加入單頻干擾噪聲對吹管音樂信號進行頻譜分析，畫出干擾前后的時域和頻域波形圖利用漢寧窗設(shè)計FIR濾波器對吹管音樂信號進行濾波比較濾波前后的時域波形和頻譜圖，并回放音樂信號，驗證是否達到去噪效果是否達到去噪效果結(jié)束NOYES圖3-1 程序設(shè)計流程圖3.2 編輯語音信號在網(wǎng)上下載一段音樂，再利用CE軟件將其轉(zhuǎn)

12、換成單聲道的.格式文件，再將此.格式音樂控制在10秒，以減少設(shè)計中的誤差。然后在Matlab軟件平臺下，利用函數(shù)wavread對語音信號進行采樣，記住采樣頻率和采樣點數(shù)。CE軟件操作界面如圖3-2所示。圖3-2 CE軟件操作界面3.3語音加噪處理采集完成后在信號中加入一個單頻噪聲，繪制原音樂信號和加噪后的音樂信號的時域和頻域的波形圖。首先，輸入原始音樂信號并播放一次。調(diào)用程序如下：x,fs,bits=wavread('h:2013DSPPurpleBambooTune.wav'); % 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣

13、率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放計算信號長度并加入噪聲。調(diào)用程序如下：N=length(x); % 計算信號x的長度fn=2100; % 單頻噪聲頻率t=0:1/fs:(N-1)/fs; % 計算時間圍，樣本數(shù)除以采樣頻率x=x(:,1)'y=x+sin(fn*2*pi*t); % 加入一個單頻噪聲sound(y,fs,bits); % 可以明顯聽出有尖銳的單頻嘯叫聲繪制原始音樂信號和加入噪聲后的音樂信號的時域和頻譜波形圖。調(diào)用程序如下：X=abs(fft(x); Y=abs(fft(y); % 對原始信

14、號和加噪信號進行fft變換X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/ N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率圍figure(1);subplot(2,2,1);plot(t,x); % 布局為2*2的四個小圖title('原始音樂信號');xlabel('時間(t)');ylabel('幅度'); %改變橫縱坐標(biāo)的圍axis(0,2,-1.5,1.5); %加上標(biāo)題和橫坐標(biāo)名稱grid on; % 加上網(wǎng)格subplot(2,2,2);plot(f,X

15、);title('原始音樂信號頻譜');xlabel('頻率(f)');ylabel('幅度譜');axis(0,3000,0,3000);grid on; subplot(2,2,3);plot(t,y);title('加入干擾后的音樂信號');xlabel('時間(t)');ylabel('幅度');axis(0,2,-1.5,1.5);grid on;subplot(2,2,4);plot(f,Y);title('加入單頻干擾后的音樂信號頻譜');xlabel('頻率(

16、f)');ylabel('幅度譜');axis(0,3000,0,3000);grid on;用繪圖命令分別畫出加噪前后的時域和頻域波形，如下圖3-3所示。圖3-3 吹管音樂信號加入單頻噪聲前后的時域與頻譜波形圖由上圖可以看到，語音信號加入單頻噪聲后的時域波形比未加之前在幅度圍有了明顯的增加，在頻譜方面可以看到除了在加了噪聲之后的頻譜圖上的2100Hz出現(xiàn)一個明顯的沖激信號外，其它地方均與未加時的原始吹管音樂信號頻譜一樣，這一現(xiàn)象表現(xiàn)在音樂播放時，可以聽見一聲尖銳的嘯叫聲。3.4濾波器設(shè)計本次課程設(shè)計中主要應(yīng)用漢寧窗設(shè)計出FIR濾波器。利用Matlab中的函數(shù)freqz

17、畫出各濾波器的頻率響應(yīng)，首先利用數(shù)字信號處理里面學(xué)過的知識，根據(jù)選定的參數(shù)，用漢寧窗函數(shù)法設(shè)計FIR數(shù)字濾波器，得到數(shù)字濾波器的參數(shù)b，a。其中b為系統(tǒng)函數(shù)的分子系數(shù)，a為系統(tǒng)函數(shù)分母系數(shù)。再調(diào)用freqz(b,a,512,fs)即可得到該濾波器的頻率響應(yīng)。主程序如下：fpd=1450;fsd=1650;fsu=2250;fpu=2350; % FIR濾波器的上下截止頻率Rp=1;As=37; % 帶阻濾波器設(shè)計指標(biāo)fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 計算上下邊帶中心頻率，和頻率間隔wcd=fcd/fs*2*

18、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(6.2*pi/dw)+1; % 計算漢寧窗設(shè)計該濾波器時需要的階數(shù)n=0:M-1; % 定義時間圍w_ham=hanning(M); % 產(chǎn)生M階的漢寧窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 調(diào)用自編函數(shù)計算理想帶阻濾波器的脈沖響應(yīng)h_bs=w_ham'.*hd_bs; % 用窗口法計算實際濾波器脈沖響應(yīng)db

19、,mag,pha,grd,w=freqz_m(h_bs,1); % 調(diào)用自編函數(shù)計算濾波器的頻率特性通過繪圖工具可得出濾波器的波形圖，如圖3-4所示。圖3-4 FIR濾波器的頻率響應(yīng)上圖為用漢寧窗設(shè)計的FIR濾波器圖，可以看出，阻帶最大衰減為-100dB，F(xiàn)IR濾波器的主瓣寬度很小，這樣可以使過渡寬度很陡，旁瓣相對于主瓣也比較小。3.5 信號濾波處理用自己設(shè)計的各濾波器分別對采集的信號進行濾波，在Matlab中，F(xiàn)IR濾波器利用函數(shù)fftfilt對信號進行濾波，IIR濾波器利用函數(shù)filter對信號進行濾波，對語音信號進行濾波后，仔細對比濾波前和濾波后的語音信號圖，得出結(jié)論。主程序如下：y_f

20、il=filter(h_bs,1,y); % 用設(shè)計好的濾波器對y進行濾波Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半由繪圖工具可以得出濾波前后的吹管音樂信號波形圖、原始吹管音樂信號波形圖和加入噪聲后的吹管音樂信號波形圖，如圖3-5所示。圖3-5濾波前后的波形圖由上圖可以看出，加噪后的吹管音樂信號經(jīng)過FIR濾波器的濾噪處理，時域和頻域圖幾乎一樣，這說明噪聲被完全濾掉，同時也說明FIR濾波器設(shè)計很理想，能滿足課設(shè)要求。3.6 結(jié)果分析語音信號經(jīng)過FIR濾波器的濾除噪聲的處理，在Matlab中，函數(shù)sound可以對聲音進行回放。其調(diào)用格式：soun

21、d(x，fs，bits)我們可以明顯感覺濾波前后的聲音有變化。聲音中刺耳的聲音沒有了，幾乎恢復(fù)成原始的聲音，但較原始的聲音更平滑一些。這說明用漢寧窗設(shè)計FIR濾波器濾掉了語音中的噪聲同時，也把原始語音的很小的一部分也濾掉了，所以回放語音的時候聽起來比以前的更加平滑，說明這段程序設(shè)計是成功的。4出現(xiàn)的問題與解決方法 在設(shè)計課程設(shè)計時，出現(xiàn)了以下幾個問題：1.、在編輯.wav音樂時，由于沒有控制好所需音樂信號的時間，導(dǎo)致結(jié)果不理想。2、對各時域、頻譜圖的圍沒有好的預(yù)計，導(dǎo)致出圖時波形效果不理想。3、因為多次改動單頻噪聲頻率值，所設(shè)計的濾波器的性能指標(biāo)沒有隨之變化，出現(xiàn)了濾波上的錯誤。4、在確定單頻

22、噪聲頻率值后，由于不能很好的掌握其它參數(shù)的調(diào)試指標(biāo)，導(dǎo)致多次調(diào)試都無法得到理想的濾波。5 結(jié)束語在此次課程設(shè)計中，我的任務(wù)是利用漢寧窗函數(shù)設(shè)計FIR濾波器，對吹管音樂信號進行濾波去噪處理。課設(shè)開始之前，我認真復(fù)習(xí)了有關(guān)窗函數(shù)特別是漢寧窗的相關(guān)知識以與濾波器的設(shè)計方法，了解課設(shè)流程。在這兩周里，我利用老師給出的模板，結(jié)合相關(guān)的專業(yè)知識，比較輕松的完成了課設(shè)的任務(wù)。不同于之前的理論課，雖然這次的課設(shè)容并不是很難，但是仍然很考驗我的動手實踐能力。此間，我得到了以下收獲：首先，在學(xué)習(xí)方面，雖然已經(jīng)學(xué)習(xí)過了DSP課程，但是如何融洽的把理論和實際結(jié)合仍然是我需要面對的問題。其次，在matlab編程方面，由

23、于老師給出了相關(guān)的程序資料，所以設(shè)計過程中并不算困難，可以說是順利過關(guān)。最后，在團結(jié)合作方面，雖然每個人都單獨分配了課題，但是總能找到與自己課題相關(guān)或者類似的同學(xué)，跟同學(xué)交流經(jīng)驗成為本次課設(shè)中的一個重要環(huán)節(jié)。在這次設(shè)計過程中，既體現(xiàn)出了我自己單獨設(shè)計的能力以與綜合運用知識的能力，又讓我體會到了學(xué)以致用、理論與實際貫通的喜悅，提高了我的團隊協(xié)作能力。同時，我也從中也能發(fā)現(xiàn)自己平時學(xué)習(xí)的不足。在此，特別感指導(dǎo)老師高明，在您的指導(dǎo)下我成功完成了本次課設(shè)的任務(wù),！參考文獻1威. MATLAB基礎(chǔ)與編程入門. ：電子科技大學(xué)，2008.12 圣勤.MATLAB7.0實用教程.：機械工業(yè)，20083 志涌

24、.精通MATLAB 6.5版M：航空航天大學(xué)，20034 程佩青.數(shù)字信號處理教程.：清華大學(xué)，20025維納·K·英格爾，約翰·G普羅克斯. 數(shù)字信號處理M. ：交通大學(xué)，2008.16小虹.信號系統(tǒng)與數(shù)字信號處理M.第版.：電子科技，2002附錄1：語音信號濾波去噪設(shè)計源程序清單% 程序名稱：DSPYFL.m% 程序功能：采用基于漢寧的窗口設(shè)計法，設(shè)計FIR濾波器對含噪語音進行濾波去噪處理。% 程序 余霏霖/% 最后修改日期：2013-3-8x,fs,bits=wavread('h:2013DSPPurpleBambooTune.wav');

25、% 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放N=length(x); % 計算信號x的長度fn=2100; % 單頻噪聲頻率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); % 對原始信號和加噪信號進

26、行fft變換X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/ N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率圍figure(1)subplot(2,2,1);plot(t,x); % 布局為2*2的四個小圖title('原始音樂信號');xlabel('時間(t)');ylabel('幅度'); %改變橫縱坐標(biāo)的圍axis(0,2,-1.5,1.5); %加上標(biāo)題和橫坐標(biāo)名稱grid on; % 加上網(wǎng)格subplot(2,2,2);plot(f,X);title(

27、'原始音樂信號頻譜');xlabel('頻率(f)');ylabel('幅度譜');axis(0,3000,0,3000);grid on; subplot(2,2,3);plot(t,y);title('加入干擾后的音樂信號');xlabel('時間(t)');ylabel('幅度');axis(0,2,-1.5,1.5);grid on;subplot(2,2,4);plot(f,Y);title('加入單頻干擾后的音樂信號頻譜');xlabel('頻率(f)')

28、;ylabel('幅度譜');axis(0,3000,0,3000);grid on;fpd=1450;fsd=1650;fpu=2350;fsu=2250; % FIR濾波器的上下截止頻率Rp=1;As=37; % 帶阻濾波器設(shè)計指標(biāo)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/

29、fs*2*pi;M=ceil(6.2*pi/dw)+1; % 計算漢寧窗設(shè)計該濾波器時需要的階數(shù)n=0:M-1; % 定義時間圍w_ham=hanning(M); % 產(chǎn)生M階的漢寧窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 調(diào)用自編函數(shù)計算理想帶阻濾波器的脈沖響應(yīng)h_bs=w_ham'.*hd_bs; % 用窗口法計算實際濾波器脈沖響應(yīng)db,mag,pha,grd,w=freqz_m(h_bs,1); % 調(diào)用自編函數(shù)計算濾波器的頻率特性figure(2)subplot(2,2,1);plot(w/pi,db);

30、axis(0,0.4,-100,20);title('以db為單位的幅度特性');xlabel('w/pi');ylabel('db'); grid on;subplot(2,2,2);plot(w/pi,mag);axis(0,0.4,-0.5,1.25);title('以線性為單位的幅度特性');xlabel('w/pi');ylabel('mag');grid on;subplot(2,2,3);plot(w,pha);title('濾波器相位響應(yīng)圖');xlabel('

31、;w/pi');ylabel('相位（pha）');axis(0,3,-4,4);grid on;subplot(2,2,4);plot(h_bs);axis(0,800,-0.2,1);title('濾波器脈沖響應(yīng)圖');xlabel('n');ylabel('h(n)');grid on; y_fil=filter(h_bs,1,y); % 用設(shè)計好的濾波器對y進行濾波Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半figure(3)subplot(3,2,1);plot(

32、t,x); axis(0,2,-1.5,1.5);title('原始音樂信號時間x');xlabel('時間（t）');ylabel('幅度');grid on;subplot(3,2,2);plot(f,X);axis(0,10000,0,1000);title('原始音樂信號幅度譜X');xlabel('頻率（f）');ylabel('幅度');grid on;subplot(3,2,3);plot(t,y);axis(0,2,-1.5,1.5);title('加干擾音樂信號時間x1');xlabel('時間（t）');ylabel('幅度');grid on;subplot(3,2,4);plot(f,Y);axis(0,10000,0,1000);title('加干擾音樂信號幅度譜X1');xlabel('頻率（f）'