數(shù)字信號課程設計--用nuttallw設計的fir濾波器

1、<p><b>　　吹管音樂濾波去噪</b></p><p>　　—用NUTTALLW設計的FIR濾波器</p><p>　　摘 要 本課程設計主要是通過使用NUTTALLW窗設計一個FIR濾波器以對音樂信號進行濾波去噪處理。本設計首先通過搜狗軟件采集一段音樂信號，依據(jù)對該信號的頻譜分析，給定相關指標。以MATLAB軟件為平臺，采用NUTTALLW窗設計

2、滿足指標的FIR濾波器，以該音樂信號進行濾波去噪處理。通過對比濾波前后的波形圖，深入了解濾波器的基本方法。通過程序調(diào)試及完善，該設計基本滿足設計要求。</p><p>　　關鍵詞 濾波去噪；FIR濾波器；NUTTALLWIN窗函數(shù)；MATLAB</p><p><b>　　1 引 言</b></p><p>　　數(shù)字濾波器是一種用來過濾時間

3、離散信號的數(shù)字系統(tǒng),它是通過對抽樣數(shù)據(jù)進行數(shù)學處理來達到頻域濾波的目的。隨著現(xiàn)代通信的數(shù)字化,數(shù)字濾波器變得更加重要。數(shù)字濾波器的種類很多，但總的來說可以分成兩大類，一類是經(jīng)典濾波器，另一類可稱為現(xiàn)代濾波器。從濾波特性方面考慮，數(shù)字濾波器可分成數(shù)字高通、數(shù)字低通、數(shù)字帶通和數(shù)字帶阻等濾波器。從實現(xiàn)方法上考慮，將濾波器分成兩種，一種稱為無限脈沖響應濾波器，簡稱IIR（Infinite Impulse Response）濾波器，另一種稱為F

4、IR（Finite Impulse Response）濾波器[1]。設計FIR數(shù)字濾波器的方法有窗函數(shù)法、頻率采樣法和等波紋最佳逼近法等。</p><p>　　1.1 課程設計目的</p><p>　　數(shù)字信號處理（Digital Signal Processing，簡稱DSP）是一門涉及許多學科而又廣泛應用于許多領域的新興學科。20世紀60年代以來，隨著計算機和信息技術的飛速發(fā)展，數(shù)字信

5、號處理技術應運而生并得到迅速的發(fā)展。在過去的二十多年時間里，數(shù)字信號處理已經(jīng)在通信等領域得到極為廣泛的應用。數(shù)字信號處理是利用計算機或?qū)Ｓ锰幚碓O備，以數(shù)字形式對信號進行采集、變換、濾波、估值、增強、壓縮、識別等處理，以得到符合人們需要的信號形式。</p><p>　　在本次課程設計中，最主要的設計是設計FIR濾波器，F(xiàn)IR濾波器的設計方法主要分為兩類：第一類是基于逼近理想濾波器器特性的方法包括窗函數(shù)法、頻率采樣法

6、、和等波紋最佳逼近法；第二類是最優(yōu)設計法。</p><p>　　本次的課程設計主要采用的是第一類設計方法，是利用漢寧窗函數(shù)法設計FIR濾波器對一段語音進行濾波去噪，通過這一過程，對濾波前后波形進行對比分析得到結論。此課程設計比較簡單，主要是將書本中的知識運用到現(xiàn)實中，并且根據(jù)自己對設計題目的理解，運用軟件編寫出程序?qū)崿F(xiàn)這一設計，也是我們對數(shù)字信號處理的原理進行驗證的一個過程。對此，也可以加深我們對所學知識的理解，

7、培養(yǎng)我們的動手能力。</p><p>　　1.2課程設計的要求</p><p>　?。?）通過利用各種不同的開發(fā)工具實現(xiàn)模擬信號數(shù)字化、信道編解碼、基帶數(shù)字信號編解碼、數(shù)字信號的調(diào)制解調(diào)和音樂信號的濾波去噪等課題，掌握數(shù)字信號的分析方法和處理方法。</p><p>　?。?）按要求編寫課程設計報告書，能正確闡述設計和實驗結果等等。</p><p&

8、gt;　　（3）通過課程設計培養(yǎng)學生嚴謹?shù)目茖W態(tài)度，認真的工作作風和團隊協(xié)作精神。</p><p>　?。?）在老師的指導下，要求每個學生獨立完成課程設計的全部內(nèi)容。</p><p><b>　　1.3設計平臺</b></p><p>　　課程設計的主要設計平臺式MATLAB 7.0。如下圖1.1所示：MATLAB是矩陣實驗室（Matrix L

9、aboratory）的簡稱，是美國MathWorks公司出品的商業(yè)數(shù)學軟件，用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算的高級技術計算語言和交互式環(huán)境，主要包括MATLAB和Simulink兩大部分。MATLAB和Mathematica、Maple并稱為三大數(shù)學軟件。它在數(shù)學類科技應用軟件中在數(shù)值計算方面首屈一指。MATLAB可以進行矩陣運算、繪制函數(shù)和數(shù)據(jù)、實現(xiàn)算法、創(chuàng)建用戶界面、連接其他編程語言的程序等，主要應用于工程計算、控制設

10、計、信號處理與通訊、圖像處理、信號檢測、金融建模設計與分析等領域。</p><p>　　圖1.1 MATLAB 7.0的設計平臺</p><p>　　MATLAB的基本數(shù)據(jù)單位是矩陣，它的指令表達式與數(shù)學、工程中常用的形式十分相似，故用MATLAB來解算問題要比用C，F(xiàn)ORTRAN等語言完成相同的事情簡捷得多，并且MathWork也吸收了像Maple等軟件的優(yōu)點,使MATLAB成為一個強大

11、的數(shù)學軟件。在新的版本中也加入了對C，F(xiàn)ORTRAN，C++ ，JAVA的支持?？梢灾苯诱{(diào)用,用戶也可以將自己編寫的實用程序?qū)氲組ATLAB函數(shù)庫中方便自己以后調(diào)用，此外許多的MATLAB愛好者都編寫了一些經(jīng)典的程序，用戶可以直接進行下載就可以用[2]。</p><p><b>　　2 設計原理</b></p><p>　　2.1 FIR濾波器</p>

12、<p>　　FIR(Finite Impulse Response)濾波器：有限長單位沖激響應濾波器，是數(shù)字信號處理系統(tǒng)中最基本的元件，它可以在保證任意幅頻特性的同時具有嚴格的線性相頻特性，同時其單位抽樣響應是有限長的，因而濾波器是穩(wěn)定的系統(tǒng)。因此，F(xiàn)IR濾波器在通信、圖像處理、模式識別等領域都有著廣泛的應用。</p><p>　　有限長單位沖激響應（FIR）濾波器有以下特點：</p>

13、<p>　　(1) 系統(tǒng)的單位沖激響應h(n)在有限個n值處不為零；</p><p>　　(2) 系統(tǒng)函數(shù)H(z)在|z|>0處收斂，極點全部在z=0處（因果系統(tǒng)）；</p><p>　　(3) 結構上主要是非遞歸結構，沒有輸出到輸入的反饋，但有些結構中（例如頻率抽樣結構）也包含有反饋的遞歸部分。</p><p>　　FIR濾波器的系統(tǒng)函數(shù)用下式表示

14、：。H(n)就是FIR濾波器的單位脈沖響應。FIR濾波器最重要的優(yōu)點就是由于不存在系統(tǒng)極點，F(xiàn)IR濾波器是絕對穩(wěn)定的系統(tǒng)。</p><p>　　相較于IIR濾波器,FIR濾波器有以下的優(yōu)點：（1）可以很容易地設計線性相位的濾波器。線性相位濾波器延時輸入信號，卻并不扭曲其相位。（2）實現(xiàn)簡單。在大多數(shù)DSP處理器，只需要對一個指令積習循環(huán)就可以完成FIR計算。（3）適合于多采樣率轉(zhuǎn)換,它包括抽取(降低采樣率)

15、，插值(增加采樣率)操作。無論是抽取或者插值，運用FIR濾波器可以省去一些計算，提高計算效率。相反，如果使用IIR濾波器，每個輸出都要逐一計算，不能省略,即使輸出要丟棄。（4）具有理想的數(shù)字特性。在實際中，所有的DSP濾波器必須用有限精度（有限bit數(shù)目）實現(xiàn)，而在IIR濾波器中使用有限精度會產(chǎn)生很大的問題，由于采用的是反饋電路，因此IIR通常用非常少的bit實現(xiàn)，設計者就能解決更少的與非理想算術有關的問題。（5）可以用小數(shù)實現(xiàn)。不

16、像IIR濾波器，F(xiàn)IR濾波器通?？赡苡眯∮?的系數(shù)來實現(xiàn)。（如果需要，F(xiàn)IR濾波器的總的增益可以在輸出調(diào)整）。當使用定點DSP的時候，這也是一個考慮因素，它能使得實現(xiàn)更加地簡單。</p><p>　　FIR濾波器的缺點在于它的性能不如同樣階數(shù)的IIR濾波器，不過由于數(shù)字計算硬件的飛速發(fā)展，這一點已經(jīng)不成為問題。再加上引入計算機輔助設計，F(xiàn)IR濾波器的設計也得到極大的簡化?；谏鲜鲈?，F(xiàn)IR濾波器比IIR濾波器的

17、應用更廣。</p><p><b>　　2.2窗口設計法</b></p><p>　　窗函數(shù)設計法的基本思想是用FIRDF逼近洗完的濾波特性。設希望逼近的濾波器的頻率響應函數(shù)為，其單位脈沖響應為表示。為了設計簡單方便，通常選擇為具有片段常數(shù)特性的理想濾波器。因此是無限長非因果序列，不能直接作為FIRDF的單位脈沖響應。窗函數(shù)設計法就是截取為有限長的一段因果序列，并用合

18、適的窗函數(shù)進行加權作為FIRDF的單位脈沖響應h(n)。</p><p>　　用窗函數(shù)法設計FIRDF的具體設計步驟如下：</p><p>　?。?）構造希望逼近的頻率響應函數(shù)。</p><p><b>　?。?）求出。</b></p><p>　　（3）加窗得到FIRDF的單位脈沖響應h(n)。h(n)=w(n)。式中

19、，w(n)稱為窗函數(shù)，其長度為。如果要求設計第一類線性相位FIRDF，則要求h(n)關于(N-1)/2點偶對稱。而關于n=t點偶對稱，所以要求t=(N-1)/2。同時要求w(n)關于(N-1)/2點偶對稱。</p><p>　　用窗函數(shù)法設計第一類線性相位FIRDF的步驟：</p><p>　?。?）選擇窗函數(shù)類型和長度，寫出才窗函數(shù)的表達式。</p><p>　　

20、根據(jù)阻帶最小衰減選擇窗函數(shù)的的類型，再根據(jù)過渡帶寬度確定所選窗函數(shù)的長度。用窗函數(shù)法設計的FIRDF通帶波紋幅度近似等于阻帶波紋幅度。一般阻帶最小衰減達到40dB以上，則通帶最大衰減就小于0.1dB。所以用窗函數(shù)法設計FIRDF時，通常只考慮阻帶最小衰減就可以了。</p><p>　　（2）構造希望逼近的頻率響應函數(shù)。</p><p>　　根據(jù)設計需要，一般選擇線性相位理想濾波器（理想低通

21、、理想高通、理想帶通、理想帶阻）。理想濾波器的截止頻率近似為最終設計的FIRDF的過渡帶中心頻率，幅度函數(shù)衰減一半（約-6dB）。所以一般取，和分別為通帶邊界頻率和阻帶邊界頻率。</p><p><b>　　（3）計算。</b></p><p>　　（4）加窗得到設計結果：h(n)=w(n)。</p><p><b>　　3設計步驟&

22、lt;/b></p><p><b>　　3.1 設計流程</b></p><p>　　3.2 采集音樂信號</p><p>　　在網(wǎng)絡中找到自己想要的歌曲文件，并用搜狗下載下來，即可對音樂文件進行文件格式轉(zhuǎn)換及聲音屬性的變化。</p><p>　　在轉(zhuǎn)換過程中，將音樂格式設置為wav格式，屬性設置為8kHz，8位

23、，單聲道，128KB/秒，然后點擊確定。然后再點擊文件/另存為，即可將該語音文件存在MATLAB文件下的work文件夾里。</p><p>　　然后在Matlab軟件平臺下，利用函數(shù)wavread對音樂信號進行采樣，記住采樣頻率和采樣點數(shù)。通過wavread函數(shù)的使用，讓我們很快理解了采樣頻率、采樣位數(shù)等概念。采集完成后在信號中加入一個單頻噪聲。</p><p><b>　　對應

24、程序如下：</b></p><p>　　[x,fs,bits]=wavread('d:\tangtao.wav'); % 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)。</p><p>　　sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放<

25、/p><p>　　N=length(x); % 計算信號x的長度</p><p>　　fn=2100; % 單頻噪聲頻率，此參數(shù)可改</p><p>　　t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數(shù)除以采樣頻率</p><p>　　x=x'; y=x+sin(fn*2*pi*t); </p><

26、;p>　　sound(y,fs,bits); % 應該可以明顯聽出有尖銳的單頻嘯叫聲</p><p>　　通過該程序得到的圖如下圖3.2.1所示</p><p>　　圖3.2.1加入噪音前后對比圖</p><p><b>　　3.3濾波器設計</b></p><p>　　在該濾波器的設計中，老師給出濾波器的性能

27、指標如下：</p><p>　　fpd=2000; fsd=2050; fsu=2150; fpu=2200; Rp=1dB; As=自定義dB;</p><p>　　截止頻率也可以任意自選，在單頻噪聲干擾附近即可。在這里，很重要的是通帶截止db值的設置。這個值一定要根據(jù)我們使用的設計濾波器的方法來設定。因為我使用的是NUTTALLWIN窗法，NUTTALLW

28、IN窗函數(shù)中，濾波器的過渡帶寬為15.4，最小阻帶衰減為108db。所以，一定要將通帶截止db值設置的小于108，所以，我將其設置為105db。</p><p>　　在這里我是使用窗函數(shù)法設計上面要求的FIR濾波器。在Matlab中，利用NUTTALLWIN窗設計FIR濾波器，利用Matlab中的函數(shù)freqz畫出該濾波器的頻率響應。</p><p>　　首先，我們利用數(shù)字信號處理里面學過

29、的知識，根據(jù)自己選定的參數(shù)，用指定的方法設計FIR濾波器，得到FIR濾波器的階數(shù)M。隨后調(diào)用nuttallwin (M)函數(shù)產(chǎn)生M階的NUTTALLWIN窗。然后，調(diào)用自編ideal_lp函數(shù)計算理想帶阻濾波器的脈沖響應。最后，再調(diào)用自編freqz_m(h_bs,1)函數(shù)即可計算得到該濾波器的頻率特性。</p><p>　　我設計濾波器過程對應的程序如下：</p><p>　　fpd=24

30、00;fsd=2450;fsu=2550;fpu=2600;Rp=1;As=105; % 帶阻濾波器設計指標</p><p>　　fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min((fsd-fpd),(fpu-fsu));</p><p>　　% 計算上下邊帶中心頻率，和頻率間隔</p><p>　　wcd=fcd/fs*2*

31、pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;</p><p>　　% 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率</p><p>　　wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;</p><p>　　M=ceil(15.4*pi/dw)+1; % 計算NUTTALLWIN窗設計該濾波器時需要的階數(shù)</p

32、><p>　　n=0:M-1; % 定義時間范圍</p><p>　　w_nut= nuttallwin (M); % 產(chǎn)生M階的NUTTALLWIN窗 </p><p>　　hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); </p><p>　　% 調(diào)用自編函數(shù)計

33、算理想帶阻濾波器的脈沖響應</p><p>　　h_bs=w_nut'.*hd_bs; % 用窗口法計算實際濾波器脈沖響應</p><p>　　[db,mag,pha,grd,w]=freqz_m(h_bs,1); % 調(diào)用自編函數(shù)計算濾波器的頻率特性</p><p>　　在將濾波器設計好，頻率特性求出來之后，我編寫如下程序，來畫出該帶阻濾波

34、器的特性圖。</p><p>　　通過程序得到的圖如下圖3.3.1所示。</p><p>　　圖 3.3.1 濾波器的特性圖</p><p>　　3.4 信號濾波處理</p><p>　　在將濾波器設計好后，我們用自己設計的帶阻濾波器對采集的語音信號進行濾波。在Matlab中，F(xiàn)IR濾波器利用函數(shù)fftfilt對信號進行濾波。</p&g

35、t;<p>　　我對信號進行濾波處理對應的程序如下：</p><p>　　y_fil=filter(h_bs,1,y);% 用設計好的濾波器對y進行濾波</p><p>　　Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半</p><p><b>　　3.5 結果分析</b></

36、p><p>　　在將加噪信號濾波之后，我們將濾波前后語音信號的波形及頻譜圖相互比較。在同一張大圖里分別繪制原始信號x，加噪信號y，濾波去噪信號y_fil的時域波形和頻譜，以便比較和分析。</p><p>　　經(jīng)過程序畫出來的三個信號的時域波形和頻譜圖如下圖3.5.1所示。</p><p>　　圖 3.5.1 三個信號的時域波形和頻譜圖比較</p><

37、;p>　　從圖3.5.1中我們可以看出，原信號與濾波去噪信號的時域圖基本相似，只有邊緣部分有點差異；原信號與濾波去噪信號的頻譜圖波形也大致相似。通過觀察可以看到，加噪信號的時域圖中大部分都被加入的噪聲給遮蓋了，加噪信號的頻譜圖中，我們可以很明顯地看到與原信號頻譜圖相比，它在2400Hz左右處有一個尖脈沖，而濾波去噪信號的頻譜圖中該尖脈沖已經(jīng)消失，波形大致與原圖相似，可見濾波去噪效果基本不錯。</p><p>

38、;　　在將三個信號的時域波形和頻譜圖比較之后，我們還要通過回放去濾波去噪音樂信號，來跟原信號相比，以檢驗濾波器的效果。在Matlab中，函數(shù)sound可以對聲音進行回放。其調(diào)用格式為：sound (x，fs，bits)。我用sound(y_fil,fs,bits)語句回放該濾波去噪信號，便可以感覺到濾波后的語音信號與原信號差不多，但仍有一點點變化。</p><p>　　4出現(xiàn)的問題及解決方法</p>

39、<p>　　在這前面一周的時間里，我們每天都是待在學校的機房上機，同時，老師跟我們進行講解，當我們有不懂的問題隨時向老師請教。在老師沒有進行講解時，我總以為這次的課程設計會有不小的難度，可是當老師跟我們講解了之后，我卻發(fā)現(xiàn)，我們的課程設計其實也不是很難，只要我們認真看書、查資料，再細心一點，并有耐心，就能解決我們所遇到的問題。</p><p>　　首先，我按照課程設計的要求在網(wǎng)絡上找到一段音樂，并對其

40、進行格式轉(zhuǎn)換及屬性變換，但是當我真正制作完成之后并存盤后，運行MATLAB后，才感覺音樂文件過于長，造成出現(xiàn)的頻譜圖很模糊，根本看不清楚，后面看了一下任務指導書及問了同學之后，截取其中的一小段，效果就好了許多。</p><p>　　在此次課程設計中，雖然沒有預計的那么難，可還是遇到了一些小小的問題。慶幸的是，能夠在老師和同學的幫助下解決了。我覺得最重要的就是要有耐心，要細致，要不懂就問，這樣，只要自己堅持，一切問

41、題都會解決的。</p><p><b>　　5 結束語</b></p><p>　　在本次課程設計中，讓我知道了以前所學《信號與系統(tǒng)》與《數(shù)字信號處理》的結合，并在實際運用中設計濾波器的過程。課程設計不僅要求對濾波器理論的研究，更重要的是實際設計中遇到的問題。</p><p>　　因為有了這次課程設計，我不得不對其設計原理進行更深一層次的理解，

42、對書中原來學到的理論，僅知道了其表面，而不知其原因。在設計中也使我對一些概念有了更深刻的認識。例如：在指標方面，我混淆了模擬指標和數(shù)字指標的概念，經(jīng)過老師的點撥，自己更加明確，而且記憶深刻。還有在課程設計中每一次的數(shù)據(jù)輸入都有其重要意義，用MATLAB編譯程序時，可以根據(jù)濾波器指標的要求實時知道對濾波器的影響。例如，編好程序后，調(diào)試成功，任意改變輸入阻帶或者通帶衰減，可以看到輸出波形的變化，改變截止頻率wc，同樣可以看到輸出波形的變化。

43、由此，對理論的理解就更加簡單方便，而且記憶力深刻。</p><p>　　除此之外，對程序的編譯不是一蹴而就的，而是經(jīng)過多次的編譯與調(diào)試。以前用MATLAB就是簡單的輸入，可以說都不是自己的勞動成果，但這次不一樣，課程設計沒有別人給你編好，而是自己寫，出錯率就大大提高了，但這是過程，學習就是在過程中進行的，經(jīng)過自己幾天的腦力勞動，再加上同學們的幫助，不僅對讀程序有了很大提高，而且自己的編譯水平也上了一個新臺階，更加

44、熟系了MATLAB的應用，也對其中的函數(shù)有了大概的了解，對其中一些函數(shù)也相當熟練，濾波器設計中用到了ideal_lp（）、freqz_m（）等很重要的函數(shù)，雖然可以用fir1（）、fir2（）可以直接調(diào)用，但那樣就達不到對真正理論設計過程的 理解和運用。</p><p>　　完成整個設計過程后，學到的東西已經(jīng)不僅僅上面的那些東西，還有就是同學之間的共同努力和探討和設計過程中的每一個細節(jié)，也許每一個細節(jié)的錯誤就可能

45、導致結局的失敗，所以我認為這次收獲最大的莫過于靜心，學習不能急，一定要冷靜，心無旁騖，不放過任何一個細節(jié)，就能帶來凱旋的消息。</p><p>　　在此，感謝老師一學年以來對我們不僅學習上的關心，還有生活中的照顧，我也不會辜負老師的期望，繼續(xù)努力深造。</p><p><b>　　參考文獻</b></p><p>　　[1] 程佩青．數(shù)字信號

46、處理教程〔M〕．北京：清華大學出版社，2002．</p><p>　　[2] 劉敏，魏玲．Matlab通信仿真與應用〔M〕．北京：國防工業(yè)出版社，2001．</p><p>　　［3］ Sanjit K.Mitra 著，孫洪，余翔宇 譯，《數(shù)字信號處理實驗指導書》，電子工業(yè)</p><p>　　出版社，2005 年。</p><p>　?。?/p>

47、4］ Edward W. Kamen, Bonnie S.Heck 編，《信號與系統(tǒng)基礎——應用Web 和MATLAB</p><p>　?。ǖ诙妫罚茖W出版社，2002 年。</p><p>　　［5］ 張賢達 編, 《現(xiàn)代信號處理（第二版）》，清華大學出版社，2002 年。</p><p>　　附錄1：語音信號濾波去噪——用NUTTALLWIN窗設計的濾波器

48、源程序清單</p><p>　　// 程序名稱：1.m</p><p>　　// 程序功能：在Matlab中，用NUTTALLWIN窗法設計濾波器</p><p>　　// 程序作者：湯濤</p><p>　　// 最后修改日期：2012-3-15</p><p>　　[x,fs,bits]=wavread('

49、d:\tangtao.wav');</p><p>　　sound(x,fs,bits);</p><p>　　N=length(x);</p><p><b>　　fn=2500;</b></p><p>　　t=0:1/fs:(N-1)/fs;</p><p>　　x=x';y

50、=x+sin(fn*2*pi*t);</p><p>　　sound(x,fs,bits);</p><p>　　X=abs(fft(x)); Y=abs(fft(y)); % 對原始信號和加噪信號進行fft變換,取幅度譜</p><p>　　X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分</p><p>　

51、　deltaf=fs/N; % 計算頻譜的譜線間隔</p><p>　　f=0:deltaf:fs/2-deltaf;</p><p>　　subplot(2,2,1);plot(t,x); grid on;xlabel('時間s');ylabel('幅度');title('原始音樂信號'); axis([0,3,-2,2])

52、;</p><p>　　subplot(2,2,2);plot(f,X);grid on; xlabel('頻率Hz');ylabel('幅度譜'); title('原始音樂信號頻譜圖'); axis([0,3000,0,800]);subplot(2,2,3);plot(t,y);grid on; xlabel('時間s');ylabel(

53、9;幅度'); title('加入噪聲后信號'); axis([0,3,-2,2]);subplot(2,2,4);plot(f,Y);grid on; xlabel('頻率Hz');ylabel('幅度譜'); title('加入噪聲后信號頻譜圖'); axis([0,3000,0,800]);</p><p>　　fpd=2400;f

54、sd=2450;fsu=2550;fpu=2600;Rp=1;As=105; % 帶阻濾波器設計指標</p><p>　　fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min((fsd-fpd),(fpu-fsu));</p><p>　　% 計算上下邊帶中心頻率，和頻率間隔</p><p>　　wcd=fcd/fs*2*pi;w

55、cu=fcu/fs*2*pi;dw=df/fs*2*pi;</p><p>　　% 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率</p><p>　　wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;</p><p>　　M=ceil(15.4*pi/dw)+1; % 計算NUTTALLWIN窗設計該濾波器時需要的階數(shù)</p>

56、<p>　　n=0:M-1; % 定義時間范圍</p><p>　　w_nut= nuttallwin (M); % 產(chǎn)生M階的NUTTALLWIN窗 </p><p>　　hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); </p><p>　　% 調(diào)用自編函數(shù)計算理想帶

57、阻濾波器的脈沖響應</p><p>　　h_bs=w_nut'.*hd_bs; % 用窗口法計算實際濾波器脈沖響應</p><p>　　[db,mag,pha,grd,w]=freqz_m(h_bs,1); % 調(diào)用自編函數(shù)計算濾波器的頻率特性subplot(2,2,1);plot(w,db); grid on;xlabel('頻率');ylabel

58、('db');title('濾波器幅度響應圖');</p><p>　　subplot(2,2,2);plot(w,mag); grid on;xlabel('頻率');ylabel('幅度'); title('濾波器幅度響應圖');</p><p>　　subplot(2,2,3);plot(w,pha);

59、grid on;xlabel('頻率');ylabel('相位'); title('濾波器相位響應圖');</p><p>　　subplot(2,2,4);stem(h_bs);grid on;xlabel('頻率');ylabel('脈沖響應'); title('濾波器脈沖響應圖');</p><

60、;p>　　y_fil=filter(h_bs,1,y);% 用設計好的濾波器對y進行濾波</p><p>　　Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半</p><p>　　subplot(3,2,1);plot(t,x);xlabel('時間(t)');ylabel('幅度'); title(

61、9;原始音樂信號');axis([0,3,-2,2]);subplot(3,2,2);plot(f,X);xlabel('頻率(f)');ylabel('幅度譜');title('原始音樂信號幅度譜');axis([0,3000,0,800]);subplot(3,2,3);plot(t,y);xlabel('時間(t)');ylabel('幅度'

62、;);title('加干擾后的音樂信號'); axis([0,3,-2,2]);subplot(3,2,4);plot(f,Y);xlabel('頻率(f)');ylabel('幅度譜');title('加干擾后的音樂信號幅度譜');axis([0,3000,0,800]);subplot(3,2,5);plot(t,y_fil);xlabel('時間(t)&#

63、39;);ylabel('幅度'); title('濾波后音樂信號'); axis([0,3,-2,2]);subplot(3,2,6);plot(f,Y_fil);</p><p>　　subplot(211);plot(f,20*log10(abs(Y)/max(abs(X))));axis tight;grid on; xlabel('頻率');ylabel

64、('幅度（對數(shù)坐標）');title('加干擾后音樂信號幅度譜');</p><p>　　subplot(212);plot(f,20*log10(abs(Y_fil)/max(abs(X))));axis tight;grid on; xlabel('頻率');ylabel('幅度（對數(shù)坐標）');title('濾波后音樂信號幅度譜'