第一篇：數(shù)字信號(hào)處理實(shí)驗(yàn)-FFT的實(shí)現(xiàn)
實(shí)
驗(yàn)
報(bào)
告
學(xué)生姓名：
學(xué) 號(hào)：
指導(dǎo)教師：
一、實(shí)驗(yàn)室名稱：數(shù)字信號(hào)處理實(shí)驗(yàn)室
二、實(shí)驗(yàn)項(xiàng)目名稱：FFT的實(shí)現(xiàn)
三、實(shí)驗(yàn)原理：
一．FFT算法思想：
1．DFT的定義：
對(duì)于有限長(zhǎng)離散數(shù)字信號(hào){x[n]}，0 ? n ? N-1，其離散譜{x[k]}可以由離散付氏變換（DFT）求得。DFT的定義為：
N?1X[k]?通常令e?j2?N?x[n]en?0?j2?Nnk，k=0,1,…N-1 ?WN，稱為旋轉(zhuǎn)因子。
2．直接計(jì)算DFT的問(wèn)題及FFT的基本思想：
由DFT的定義可以看出，在x[n]為復(fù)數(shù)序列的情況下，完全直接運(yùn)算N點(diǎn)DFT需要（N-1）2次復(fù)數(shù)乘法和N（N-1）次加法。因此，對(duì)于一些相當(dāng)大的N值（如1024）來(lái)說(shuō)，直接計(jì)算它的DFT所作的計(jì)算量是很大的。
FFT的基本思想在于，將原有的N點(diǎn)序列分成兩個(gè)較短的序列，這些序列的DFT可以很簡(jiǎn)單的組合起來(lái)得到原序列的DFT。例如，若N為偶數(shù)，將原有的N
22點(diǎn)序列分成兩個(gè)（N/2）點(diǎn)序列，那么計(jì)算N點(diǎn)DFT將只需要約[(N/2)·2]=N/2次復(fù)數(shù)乘法。即比直接計(jì)算少作一半乘法。因子（N/2）2表示直接計(jì)算（N/2）點(diǎn)DFT所需要的乘法次數(shù)，而乘數(shù)2代表必須完成兩個(gè)DFT。上述處理方法可以反復(fù)使用，即（N/2）點(diǎn)的DFT計(jì)算也可以化成兩個(gè)（N/4）點(diǎn)的DFT（假定N/2為偶數(shù)），從而又少作一半的乘法。這樣一級(jí)一級(jí)的劃分下去一直到最后就劃分成兩點(diǎn)的FFT運(yùn)算的情況。
3．基2按時(shí)間抽取（DIT）的FFT算法思想：
設(shè)序列長(zhǎng)度為N?2L，L為整數(shù)（如果序列長(zhǎng)度不滿足此條件，通過(guò)在后面補(bǔ)零讓其滿足）。
將長(zhǎng)度為N?2L的序列x[n](n?0,1,...,N?1)，先按n的奇偶分成兩組：
x[2r]?x1[r]x[2r?1]?x2[r]，r=0,1,…,N/2-1 DFT化為：
N?1N/2?1N/2?1X[k]?DFT{x[n]}?N/2?1?n?0x[n]WnkN?2rk?r?0x[2r]W2rkN??r?0x[2r?1]WN(2r?1)kN/2?1???r?0N/2?1x1[r]Wx1[r]W2rkN?W?WkN?r?0N/2?1x2[r]WN
?r?0rkN/2kN?r?0x2[r]WN/22rkrk上式中利用了旋轉(zhuǎn)因子的可約性，即：WNN/?21NrkN?/21rk?WN/2。又令
rkX1[k]??r?0x[1r]W,/X2[k]?2?r?0x[r]WN2，則上式可以寫(xiě)成： /2X[k]?X1[k]?WNX2[k](k=0,1,…,N/2-1)
k可以看出，X1[k],X2[k]分別為從X[k]中取出的N/2點(diǎn)偶數(shù)點(diǎn)和奇數(shù)點(diǎn)序列的N/2點(diǎn)DFT值，所以，一個(gè)N點(diǎn)序列的DFT可以用兩個(gè)N/2點(diǎn)序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了X[k]前N/2點(diǎn)的DFT值，還需要繼續(xù)利用X1[k],X2[k]表示X[k]的后半段本算法推導(dǎo)才完整。利用旋轉(zhuǎn)因子的周期性，有：WN/2?WN/2X1[N2N/2?1rkr(k?N/2)，則后半段的DFT值表達(dá)式：
rk?k]??r?0x1[r]W2N/2r(N?k)N/2?1??r?0x1[r]WN/2?X1[k]，同樣，X2[N2?k]?X2[k]
（k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表達(dá)式獲得，中間還利用到WN(N2?k)N?WN2Wk得到后半段的X[k]值表達(dá)式??W，k為：X[k]?X1[k]?WNkX2[k]（k=0,1,…,N/2-1）。
這樣，通過(guò)計(jì)算兩個(gè)N/2點(diǎn)序列x1[n],x2[n]的N/2點(diǎn)DFTX1[k],X2[k]，可以組合得到N點(diǎn)序列的DFT值X[k]，其組合過(guò)程如下圖所示：
X1[k] X1[k]?WNkX2[k]
X2[k] WNnk-1 X1[k]?WNkX2[k]
比如，一個(gè)N = 8點(diǎn)的FFT運(yùn)算按照這種方法來(lái)計(jì)算FFT可以用下面的流程圖來(lái)表示：
x(0)W0x(1)W0x(2)W0x(3)W2W0W1W0x(5)W0x(6)W0x(7)W2X(7)W3X(6)W2X(5)X(3)X(2)X(1)X(0)x(4)X(4)
4．基2按頻率抽?。―IF）的FFT算法思想：
設(shè)序列長(zhǎng)度為N?2L，L為整數(shù)（如果序列長(zhǎng)度不滿足此條件，通過(guò)在后面補(bǔ)零讓其滿足）。
在把X[k]按k的奇偶分組之前，把輸入按n的順序分成前后兩半：
N?1N/2?1nkNN?1X[k]?DFT{x[n]}?N/2?1N/2?1?x[n]Wn?0?(n??n?0N2)kx[n]WnkN??n?N/2x[n]WNnk??n?0N/2?1x[n]WnkN??n?0x[n?NkN2]WNnk
?N?n?0[x[n]?x[n?N2NkN2]W2N]?WN,k?0,1,...,N?1因?yàn)閃2N??1，則有WX[k]???(?1)，所以：
kkN/2?1?n?0[x[n]?(?1)x[n?N2]]?WN,k?0,1,...,N?1
nk按k的奇偶來(lái)討論，k為偶數(shù)時(shí)：
N/2?1X[2r]??n?0[x[n]?x[n?N2]]?WN,k?0,1,...,N?1 N22rnN/2?1k為奇數(shù)時(shí)：X[2r?1]?前面已經(jīng)推導(dǎo)過(guò)WNN/2?1?n?0[x[n]?x[n?]]?WN(2r?1)n,k?0,1,...,N?1
2rk?WN/2，所以上面的兩個(gè)等式可以寫(xiě)為：
N2]]?WN/2,r?0,1,...,N/2?1 N2rnrkX[2r]??n?0[x[n]?x[n?N/2?1X[2r?1]??n?0{[x[n]?x[n?]]?WN}WN/2,r?0,1,...,N/2?1
nnr通過(guò)上面的推導(dǎo)，X[k]的偶數(shù)點(diǎn)值X[2r]和奇數(shù)點(diǎn)值X[2r?1]分別可以由組合而成的N/2點(diǎn)的序列來(lái)求得，其中偶數(shù)點(diǎn)值X[2r]為輸入x[n]的前半段和后半段之和序列的N/2點(diǎn)DFT值，奇數(shù)點(diǎn)值X[2r?1]為輸入x[n]的前半段和后半段之差再與WN相乘序列的N/2點(diǎn)DFT值。
令x1[n]?x[n]?x[n?N/2?1nN2]，x2[n]?[x[n]?x[n?N/2?1N2]]?WN，則有：
nX[2r]??n?0x1[n]?WrnN/2,X[2r?1]??n?0x2[n]?WrnN/2,r?0,1,...,N2?1
這樣，也可以用兩個(gè)N/2點(diǎn)DFT來(lái)組合成一個(gè)N點(diǎn)DFT，組合過(guò)程如下圖所示：
x[n] x[n]?x[n?N2]
x[n?N2]-1 WNn [x[n]?x[n?N2]]WNn
二．在FFT計(jì)算中使用到的MATLAB命令：
函數(shù)fft(x)可以計(jì)算R點(diǎn)序列的R點(diǎn)DFT值；而fft(x,N)則計(jì)算R點(diǎn)序列的N點(diǎn)DFT，若R>N，則直接截取R點(diǎn)DFT的前N點(diǎn)，若R 四、實(shí)驗(yàn)?zāi)康模?/p> 離散傅氏變換（DFT）的目的是把信號(hào)由時(shí)域變換到頻域，從而可以在頻域分析處理信息，得到的結(jié)果再由逆DFT變換到時(shí)域。FFT是DFT的一種快速算法。在數(shù)字信號(hào)處理系統(tǒng)中，F(xiàn)FT作為一個(gè)非常重要的工具經(jīng)常使用，甚至成為DSP運(yùn)算能力的一個(gè)考核因素。 本實(shí)驗(yàn)通過(guò)直接計(jì)算DFT，利用FFT算法思想計(jì)算DFT，以及使用MATLAB函數(shù)中的FFT命令計(jì)算離散時(shí)間信號(hào)的頻譜，以加深對(duì)離散信號(hào)的DFT變換及FFT算法的理解。 五、實(shí)驗(yàn)內(nèi)容： a)計(jì)算實(shí)數(shù)序列x(n)?cos5?16n,0?n?256的256點(diǎn)DFT。 b)計(jì)算周期為1kHz的方波序列（占空比為50％，幅度取為＋/-512，采樣頻率為25kHz，取256點(diǎn)長(zhǎng)度）256點(diǎn)DFT。 六、實(shí)驗(yàn)器材（設(shè)備、元器件）： 安裝MATLAB軟件的PC機(jī)一臺(tái)，DSP實(shí)驗(yàn)演示系統(tǒng)一套。 七、實(shí)驗(yàn)步驟： （1）先利用DFT定義式，編程直接計(jì)算2個(gè)要求序列的DFT值。 （2）利用MATLAB中提供的FFT函數(shù)，計(jì)算2個(gè)要求序列的DFT值。（3）(拓展要求)不改變序列的點(diǎn)數(shù)，僅改變DFT計(jì)算點(diǎn)數(shù)（如變?yōu)橛?jì)算1024點(diǎn)DFT值），觀察畫(huà)出來(lái)的頻譜與前面頻譜的差別，并解釋這種差別。通過(guò)這一步驟的分析，理解頻譜分辨力的概念，解釋如何提高頻譜分辨力。 （4）利用FFT的基本思想（基2－DIT或基2－DIF），自己編寫(xiě)FFT計(jì)算函數(shù)，并用該函數(shù)計(jì)算要求序列的DFT值。并對(duì)前面3個(gè)結(jié)果進(jìn)行對(duì)比。 （5）(拓展要求)嘗試對(duì)其他快速傅立葉變換算法（如Goertzel算法）進(jìn)行MATLAB編程實(shí)現(xiàn)，并用它來(lái)計(jì)算要求的序列的DFT值。并與前面的結(jié)果進(jìn)行對(duì)比。 （6）（拓展要求）在提供的DSP實(shí)驗(yàn)板上演示要求的2種序列的FFT算法（基2－DIT），用示波器觀察實(shí)際計(jì)算出來(lái)的頻譜結(jié)果，并與理論結(jié)果對(duì)比。 八、實(shí)驗(yàn)數(shù)據(jù)及結(jié)果分析： 程序：（1）對(duì)要求的2種序列直接進(jìn)行DFT計(jì)算的程序 （2）對(duì)要求的2種序列進(jìn)行基2－DIT和基2－DIF FFT算法程序（3）對(duì)要求的2種序列用MATLAB中提供的FFT函數(shù)進(jìn)行計(jì)算的程序 結(jié)果：（1）對(duì)2種要求的序列直接進(jìn)行DFT計(jì)算的頻域波形 （2）對(duì)2種要求的序列進(jìn)行基2－DIT和基2－DIF FFT算法頻域波形（3）對(duì)2種要求的序列用MATLAB中提供的FFT函數(shù)計(jì)算的頻域波形。（4）（拓展要求）分析利用上面的方法畫(huà)出的信號(hào)頻譜與理論計(jì)算出來(lái)的頻譜之間的差異，并解釋這種差異。 （5）（拓展要求）保持序列點(diǎn)數(shù)不變，改變DFT計(jì)算點(diǎn)數(shù)（變?yōu)?024點(diǎn)），觀察頻譜的變化，并分析這種變化，由此討論如何提高頻譜分辨力的問(wèn)題。 九、實(shí)驗(yàn)結(jié)論： 十、總結(jié)及心得體會(huì)： 十一、對(duì)本實(shí)驗(yàn)過(guò)程及方法、手段的改進(jìn)建議： 數(shù)字信號(hào)處理 說(shuō) 明 書(shū) 目錄 一．摘要…………………………………3 二．課程設(shè)計(jì)目的………………………3 三．設(shè)計(jì)內(nèi)容……………………………3 四．設(shè)計(jì)原理……………………………4 4.1．語(yǔ)音信號(hào)的采集…………………………….4 4.2．濾波器……………………………………….4 4．21．IIR濾波器原理…………………………………….4 4.22．FIR濾波器原理………………………………………5 五．設(shè)計(jì)步驟……………………………6 5.1錄制女音………………………………………6 5.2采樣語(yǔ)音信號(hào)并畫(huà)出時(shí)域波形和頻譜圖……7 5.3采用雙線性變換法設(shè)計(jì)IIR濾波器…………10 5.4窗函數(shù)法設(shè)計(jì)FFR濾波器………………......12 5.5用IIR濾波器對(duì)信號(hào)進(jìn)行濾波………………14 5.6用FIR濾波器對(duì)信號(hào)進(jìn)行濾波………………16 5.7男女聲語(yǔ)音信號(hào)頻譜特點(diǎn)分析………………19 5.8有背景噪聲的信號(hào)分析………………………20 六．心得體會(huì)…………………………….22 七．參考文獻(xiàn)…………………………….23 一．摘要： 這次課程設(shè)計(jì)的主要目的是綜合運(yùn)用本課程的理論知識(shí)進(jìn)行頻譜分析以及濾波器設(shè)計(jì)，通過(guò)理論推導(dǎo)得出相應(yīng)結(jié)論，并利用MATLAB或者ＤＳＰ開(kāi)發(fā)系統(tǒng)作為工具進(jìn)行實(shí)現(xiàn)，從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識(shí)，提高對(duì)所學(xué)知識(shí)的綜合應(yīng)用能力，并從實(shí)踐上初步實(shí)現(xiàn)對(duì)數(shù)字信號(hào)的處理。通過(guò)對(duì)聲音的采樣，將聲音采樣后的頻譜與濾波。 MATLAB全稱是Matrix Laboratory，是一種功能強(qiáng)大、效率高、交互性好的數(shù)值和可視化計(jì)算機(jī)高級(jí)語(yǔ)言，它將數(shù)值分析、矩陣運(yùn)算、信號(hào)處理和圖形顯示有機(jī)地融合為一體，形成了一個(gè)極其方便、用戶界面友好的操作環(huán)境。經(jīng)過(guò)多年的發(fā)展，已經(jīng)發(fā)展成為一種功能全面的軟件，幾乎可以解決科學(xué)計(jì)算中所有問(wèn)題。MATLAB軟件還提供了非常廣泛和靈活的用于處理數(shù)據(jù)集的數(shù)組運(yùn)算功能。 在本次課程設(shè)計(jì)中，主要通過(guò)MATLAB來(lái)編程對(duì)語(yǔ)音信號(hào)處理與濾波，設(shè)計(jì)濾波器來(lái)處理數(shù)字信號(hào)并對(duì)其進(jìn)行分析。 二．課程設(shè)計(jì)目的： 綜合運(yùn)用本課程的理論知識(shí)進(jìn)行頻譜分析以及濾波器設(shè)計(jì)，通過(guò)理論推導(dǎo)得出相應(yīng)結(jié)論，并利用MATLAB作為工具進(jìn)行實(shí)現(xiàn)，從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識(shí)，提高對(duì)所學(xué)知識(shí)的綜合應(yīng)用能力，并從實(shí)踐上初步實(shí)現(xiàn)對(duì)數(shù)字信號(hào)的處理。 三．設(shè)計(jì)內(nèi)容： 內(nèi)容：錄制一段個(gè)人自己的語(yǔ)音信號(hào)，并對(duì)錄制的信號(hào)進(jìn)行采樣；畫(huà)出采樣后語(yǔ)音信號(hào)的時(shí)域波形和頻譜圖；給定濾波器的性能指標(biāo)，采用窗函數(shù) 法和雙線性變換法設(shè)計(jì)濾波器，并畫(huà)出濾波器的頻率響應(yīng)；然后用自己設(shè)計(jì)的濾波器對(duì)采集的信號(hào)進(jìn)行濾波，畫(huà)出濾波后信號(hào)的時(shí)域波形和頻譜，并對(duì)濾波前后的信號(hào)進(jìn)行對(duì)比，分析信號(hào)的變化；回放語(yǔ)音信號(hào)；換一個(gè)與你性別相異的人錄制同樣一段語(yǔ)音內(nèi)容，分析兩段內(nèi)容相同的語(yǔ)音信號(hào)頻譜之間有什么特點(diǎn)；再錄制一段同樣長(zhǎng)時(shí)間的背景噪聲疊加到你的語(yǔ)音信號(hào)中，分析疊加前后信號(hào)頻譜的變化，設(shè)計(jì)一個(gè)合適的濾波器，能夠把該噪聲濾除。 四．設(shè)計(jì)原理： 4.1．語(yǔ)音信號(hào)的采集 熟悉并掌握MATLAB中有關(guān)聲音（wave）錄制、播放、存儲(chǔ)和讀取的函數(shù)，在MATLAB環(huán)境中，有關(guān)聲音的函數(shù)有： a：y=wavrecord(N,fs,Dtype);利用系統(tǒng)音頻輸入設(shè)備錄音，以fs為采樣頻率，默認(rèn)值為11025，即以11025HZ進(jìn)行采樣。Dtype為采樣數(shù)據(jù)的存儲(chǔ)格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度進(jìn)行采樣，其它三種都是16位采樣結(jié)果轉(zhuǎn)換為指定的MATLAB數(shù)據(jù)； b：wavplay(y,fs)；利用系統(tǒng)音頻輸出設(shè)備播放，以fs為播放頻率，播放語(yǔ)音信號(hào)y； c：wavwrite((y,fs,wavfile);創(chuàng)建音頻文件； d：y=wavread(file);讀取音頻文件； 關(guān)于聲音的函數(shù)還有sound()；soundsc()；等。4.2濾波器： 4．21.IIR濾波器原理 沖激響應(yīng)不變法是使數(shù)字濾波器在時(shí)域上模擬濾波器，但是它們的缺點(diǎn)是產(chǎn)生頻率響應(yīng)的混疊失真，這是由于從s平面到z平面是多值的映射關(guān)系所造成的。 雙線性變換法是使數(shù)字濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)相似的一種變換方法。為了克服多值映射這一缺點(diǎn)，我們首先把整個(gè)s平面壓縮變換到某一中介的s1平面的一條橫帶里，再通過(guò)變換關(guān)系將此橫帶變換到整個(gè)z平面上去，這樣就使得s平面與z平面是一一對(duì)應(yīng)的關(guān)系，消除了多值變換性，也 就消除了頻譜混疊現(xiàn)象。 雙線性法設(shè)計(jì)IIR數(shù)字濾波器的步驟： 1）將數(shù)字濾波器的頻率指標(biāo){ ?k}由Wk＝（2/T）*tan(wk),轉(zhuǎn)換為模擬濾波器的頻率指標(biāo){?k}.2)由模擬濾波器的指標(biāo)設(shè)計(jì)H(s).3)由H(s)轉(zhuǎn)換為H(z)21?z?1H(z)?H(s)s?T1?z?1 4.22.FIR濾波器原理 FIR濾波器與IIR濾波器特點(diǎn)不同，IIR濾波器的相位是非線性的，若需線性相位則要采用全通網(wǎng)絡(luò)進(jìn)行相位校正。而有限長(zhǎng)單位沖激響應(yīng)（FIR）數(shù)字濾波器就可以做成具有嚴(yán)格的線性相位，同時(shí)又可以具有任意的幅度特性。 由于FIR系統(tǒng)的沖激響應(yīng)就是其系統(tǒng)函數(shù)各次項(xiàng)的系數(shù)，所以設(shè)計(jì)FIR濾波器的方法之一可以從時(shí)域出發(fā)，截取有限長(zhǎng)的一段沖激響應(yīng)作為H（z）的系數(shù)，沖激響應(yīng)長(zhǎng)度N就是系統(tǒng)函數(shù)H（z）的階數(shù)。只要N足夠長(zhǎng)，截取的方法合理，總能滿足頻域的要求。這種時(shí)域設(shè)計(jì)、頻域檢驗(yàn)的方法一般要反復(fù)幾個(gè)回合，不像IIR DF設(shè)計(jì)靠解析公式一次計(jì)算成功。給出的理想濾波器頻率響應(yīng)是，它是w的周期函數(shù)，周期 由傅立葉反變換導(dǎo)出，即 hd(n)?1Hd(ejw)ejwndw?2?，再將hd(n)與窗函數(shù)，因此可展開(kāi)成傅氏級(jí)數(shù)w(n)相乘就可以得到h(n)。、的計(jì)算可采用傅氏變換的現(xiàn)成公式和程序，窗函數(shù)也是現(xiàn)成的。但整個(gè)設(shè)計(jì)過(guò)程不能一次完成，因?yàn)榇翱陬愋秃痛笮〉倪x擇沒(méi)有解析公式可一次算，整個(gè)設(shè)計(jì)可用計(jì)算機(jī)編程來(lái)做。 窗函數(shù)的傅式變換W(ejω)的主瓣決定了H(ejω)過(guò)渡帶寬。W(ejω)的旁瓣大小和多少?zèng)Q定了H(ejω)在通帶和阻帶范圍內(nèi)波動(dòng)幅度，常用的幾種窗函數(shù)有： 矩形窗 w(n)=RN(n)； Hanning窗 ； Hamming窗 ； Blackmen窗 ； Kaiser窗。 式中Io(x)為零階貝塞爾函數(shù)。 五．設(shè)計(jì)步驟： 5.1錄制女音： 利用MATLAB中的函數(shù)錄制聲音。function nvyin()fs=11025; %采樣頻率 str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 開(kāi)始錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double'); %錄制聲音3秒 str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 錄音結(jié)束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 播放錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs); %播放錄音 str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 播放錄音結(jié)束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音'); %聲音的存儲(chǔ) 5.2采樣語(yǔ)音信號(hào)并畫(huà)出時(shí)域波形和頻譜圖 讀取語(yǔ)音信號(hào)，畫(huà)出其時(shí)域波形和頻譜圖，與截取后的語(yǔ)音信號(hào)的時(shí)域波形和頻譜圖比較，觀察其變化。程序如下： [x,fs,bits]=wavread('女音.wav'); %讀取聲音 N=length(x); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; %信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N; %采樣間隔 n=1:N/2; %取信號(hào)的一半 figure(1);subplot(2,2,1); %把畫(huà)圖區(qū)域劃分為2行2列，指定第一個(gè)圖 plot(t, x); %畫(huà)出聲音采樣后的時(shí)域波形 title('原女音信號(hào)的時(shí)域波形'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('時(shí)間/t');ylabel('振幅/A');grid; %添加網(wǎng)格 y=fft(x); %對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0; %頻域采樣點(diǎn) subplot(2,2,3); %把畫(huà)圖區(qū)域劃分為2行2列，指定第三個(gè)圖 plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('FFT變換后聲音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid; %添加網(wǎng)格 subplot(2,2,4); %把畫(huà)圖區(qū)域劃分為2行2列，指定第四個(gè)圖 if y~=0 %判斷指數(shù)是否為0 plot(k,20*log10(abs(y(n)))); %畫(huà)信號(hào)頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid; %添加網(wǎng)格 %實(shí)際發(fā)出聲音落后錄制動(dòng)作半拍的現(xiàn)象的解決 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]); %截取語(yǔ)音信號(hào) N=length(x1); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; %信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N; %采樣間隔 n=1:N/2; %取信號(hào)的一半 figure(2);subplot(2,2,1); %把畫(huà)圖區(qū)域劃分為2行2列，指定第一個(gè)圖 plot(t,x1); %畫(huà)出聲音采樣后的時(shí)域波形 title('截取后女音信號(hào)的時(shí)域波形'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('時(shí)間/t');ylabel('振幅/A');grid; %添加網(wǎng)格 y1=fft(x1); %對(duì)信號(hào)做N點(diǎn)FFT變換 subplot(2,2,3); %把畫(huà)圖區(qū)域劃分為2行2列，指定第三個(gè)圖 k=(n-1)*f0; %頻域采樣點(diǎn) plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('FFT變換后聲音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid; %添加網(wǎng)格 subplot(2,2,4); %把畫(huà)圖區(qū)域劃分為1行2列，指定第二個(gè)圖 if y1~=0 %判斷指數(shù)是否為0 plot(k,20*log10(abs(y1(n)))); %畫(huà)信號(hào)頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid; %添加網(wǎng)格 原女音信號(hào)的時(shí)域波形10.5A/幅0振-0.5-10123時(shí)間/tFFT變換后聲音的頻譜特性FFT變換后聲音的頻譜特性30050A200貝/值分/幅0幅100振00200040006000-***頻率/HzHz 截取后女音信號(hào)的時(shí)域波形10.5振幅/A0-0.5-10123FFT變換后聲音的頻譜特性50時(shí)間/tFFT變換后聲音的頻譜特性300200振幅/分貝幅值/A01000020004000頻率/Hz6000-5002000Hz40006000 結(jié)果分析： 由原女音信號(hào)的時(shí)域波形可知錄取開(kāi)始時(shí)實(shí)際發(fā)出聲音大概落后3500個(gè)采樣點(diǎn)，我們把前3500點(diǎn)去除即可解決實(shí)際發(fā)出聲音落后錄制動(dòng)作半拍的現(xiàn)象。由原女音的的頻譜圖和截取后聲音的頻譜圖可看出，對(duì)聲音的截取并不會(huì)影響它們頻譜分布。 5.3采用雙線性變換法設(shè)計(jì)IIR濾波器： 人的聲音頻率一般在（1~~4）kHZ之間，則我們只需要設(shè)計(jì)一個(gè)帶通濾波器即可濾去聲音頻帶以外的無(wú)用噪聲，得到比較清晰的聲音。根據(jù)聲音的頻譜圖分析，設(shè)計(jì)一個(gè)帶通濾波器性能指標(biāo)如下： fp1＝1000 Hz，fp2＝3000 Hz，fsc1＝500 Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000 程序如下： %iir帶通的代碼： %w=2*pi*f/fs Ap=1; %通帶波紋系數(shù) Az=100; %最小阻帶衰減 wp=[0.2 0.6]; %歸一化通帶數(shù)字截止頻率 wz=[0.1 0.7]; %歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az); %估計(jì)契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn); %N指定濾波器階數(shù)，wn歸一化 截 %止頻率，Ap通帶波動(dòng) [h,w]=freqz(b,a); %求數(shù)字濾波器的復(fù)頻率響應(yīng) figure(1);subplot(2,1,1);plot(w/pi,abs(h)); %繪制數(shù)字濾波器的頻譜圖 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)');subplot(2,1,2);if abs(h)~=0 %判斷指數(shù)是否為0 plot(w/pi,20*log10(abs(h))); %繪制數(shù)字濾波器頻譜的分貝圖 end grid;xlabel('omega/pi');ylabel('振幅(分貝)');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)');契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)1振幅（幅值）0.5000.10.20.50.60.70.8?/?契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)0.30.40.910振幅(分貝)-200-400-60000.10.20.30.40.5?/?0.60.70.80.91 5.4窗函數(shù)法設(shè)計(jì)FFR濾波器 線性相位FIR濾波器通常采用窗函數(shù)法設(shè)計(jì)。窗函數(shù)法設(shè)計(jì)FIR濾波器的基本思想是：根據(jù)給定的濾波器技術(shù)指標(biāo)，選擇濾波器長(zhǎng)度N和窗函數(shù)ω（n），使其具有最窄寬度的主瓣和最小的旁瓣。其核心是從給定的頻率特性，通過(guò)加窗確定有限長(zhǎng)單位脈沖響應(yīng)序列h(n)。工程中常用的窗函數(shù)共有6種，即矩形窗、巴特利特（Bartlett）窗、漢寧（Hanning）窗、漢明（Hamming）窗、布萊克曼（Blackman）窗和凱澤（Kaiser）窗。 這次設(shè)計(jì)我采用的是布萊克曼來(lái)設(shè)計(jì)給定數(shù)字帶通濾波器的參數(shù)如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下： Ap=1; %通帶波紋系數(shù) Az=100; %最小阻帶衰減 fs=10000; %采樣頻率 wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2)); %---取兩個(gè)過(guò)渡帶中的小者 N0=ceil(2*5.5*pi/deltaW); %---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2); %---確保N為奇數(shù) hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想帶通濾波器 wdWindow=blackman(N); %布拉克曼窗 hr=wdWindow.*hdWindow'; %點(diǎn)乘 n=0:N-1; %階數(shù) subplot(2,2,1);stem(n,wdWindow); %繪制布拉克曼窗時(shí)域波形 xlabel('時(shí)間');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1); %求濾波器頻率響應(yīng) subplot(2,2,3);plot(W/pi,abs(H)) %繪制濾波器頻域波形 xlabel('omega/pi');ylabel('振幅');title('FIR帶通濾波器幅頻特性');subplot(2,2,4); if abs(H)~=0 %判斷指數(shù)是否為0 plot(W/pi,20*log10(abs(H))); %畫(huà)濾波器頻譜的分貝圖 end xlabel('omega/');ylabel('振幅/分貝');title('FIR帶通濾波器幅頻特性');grid; %添加網(wǎng)格 %---ideallp（）函數(shù)（非系統(tǒng)自有函數(shù)）在系統(tǒng)安裝目錄的WORK子目錄ideallp.m function hd = ideallp(wc,N);% 理想低通濾波器的脈沖響應(yīng)子程序 % hd = 點(diǎn)0 到 N-1之間的理想脈沖響應(yīng) % wc = 截止頻率(弧度)% N = 理想濾波器的長(zhǎng)度 tao =(N-1)/2; % 理想脈沖響應(yīng)的對(duì)稱中心位置 n = [0:(N-1)]; % 設(shè)定脈沖響應(yīng)長(zhǎng)度 m = n-tao + eps; % 加一個(gè)小數(shù)以避免零作除數(shù) hd = sin(wc*m)./(pi*m); % 理想脈沖響應(yīng) 布拉克曼窗1振幅0.500406080時(shí)間FIR帶通濾波器幅頻特性500振幅/分貝20FIR帶通濾波器幅頻特性1.51振幅-50-100-15000.5?/10.5000.5?/?1 5.5用IIR濾波器對(duì)信號(hào)進(jìn)行濾波 用自己設(shè)計(jì)的IIR濾波器分別對(duì)采集的信號(hào)進(jìn)行濾波，在Matlab中，IIR濾波器利用函數(shù)filter對(duì)信號(hào)進(jìn)行濾波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; %信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N; %采樣間隔 n=1:N/2; %取信號(hào)的一半 y=fft(x); %對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0; %頻域采樣點(diǎn) subplot(2,1,1); %把畫(huà)圖區(qū)域劃分為2行1列，指定第一個(gè)圖 plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid; %iir帶通的代碼： Ap=1; %通帶波紋系數(shù) Az=100; %最小阻帶衰減 wp=[0.2 0.6]; %歸一化通帶數(shù)字截止頻率 wz=[0.1 0.7]; %歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az); %估計(jì)契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn); %N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動(dòng) x1=filter(b,a,x); %對(duì)聲音濾波 wavplay(x1)wavwrite(x1,'IIR濾波后女音.wav');N=length(x1); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; %信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N; %采樣間隔 n=1:N/2; %取信號(hào)的一半 y=fft(x1); %對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0; %頻域采樣點(diǎn) subplot(2,1,2); %把畫(huà)圖區(qū)域劃分為2行1列，指定第一個(gè)圖 plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('l濾波后女音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid; 濾波前女音的頻譜特性300幅值/A***030004000頻率/Hz濾波后女音的頻譜特性500060006040幅值/A***0頻率/Hz400050006000 結(jié)果分析： 由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3400Hz的頻譜成分?；胤耪Z(yǔ)音信號(hào)，由于低頻和高頻成分被濾除，聲音變得較低沉。 5.6用FIR濾波器對(duì)信號(hào)進(jìn)行濾波 用自己設(shè)計(jì)的FIR濾波器分別對(duì)采集的信號(hào)進(jìn)行濾波，在Matlab中,FIR濾波器利用函數(shù)fftfilt對(duì)信號(hào)進(jìn)行濾波 程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; %信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N; %采樣間隔 n=1:N/2; %取信號(hào)的一半 y=fft(x); %對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0; %頻域采樣點(diǎn) subplot(2,1,1); %把畫(huà)圖區(qū)域劃分為2行1列，指定第一個(gè)圖 plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/omega');ylabel('幅值/A');grid; %FIR帶通濾波器代碼 fs=10000;wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2)); %---取兩個(gè)過(guò)渡帶中的小者 N0=ceil(2*5.5*pi/deltaW); %---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2); %---確保N為奇數(shù) hdWindow=ideallp(wc2,N)-ideallp(wc1,N);wdWindow=blackman(N);hr=wdWindow.*hdWindow';x1=fftfilt(hr,x); %對(duì)聲音濾波 wavplay(x1)wavwrite(x1,'FIR濾波后女音.wav');N=length(x1); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; %信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N; %采樣間隔 n=1:N/2; %取信號(hào)的一半 y=fft(x1); %對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0; %頻域采樣點(diǎn) subplot(2,1,2); %把畫(huà)圖區(qū)域劃分為2行1列，指定第一個(gè)圖 plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('l濾波后女音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid; 濾波前女音的頻譜特性300200幅值/A***004000頻率/?l濾波后女音的頻譜特性500060006040幅值/A20005001000***03000頻率/Hz***0 結(jié)果分析： 由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3500Hz的頻譜成分。和用IIR濾波器濾波一樣，回放語(yǔ)音信號(hào)，由于低頻和高頻成分被濾除，聲音變得較低沉。5.7男女聲語(yǔ)音信號(hào)頻譜特點(diǎn)分析 換一個(gè)男音錄制同樣一段語(yǔ)音內(nèi)容，分析兩段內(nèi)容相同的語(yǔ)音信號(hào)頻譜之間有什么特點(diǎn)。程序如下： [x,fs,bits]=wavread('女音.wav'); %讀取聲音 N=length(x); %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs; f0=fs/N; n=1:N/2; y=fft(x); k=(n-1)*f0; subplot(2,1,1); plot(k,abs(y(n))); title('FFT變換后女音的頻譜特性');xlabel('頻率/omega');ylabel('幅值/A');grid; [x,fs,bits]=wavread('明明.wav'); N=length(x); t=(1:N)/fs; f0=fs/N; n=1:N/2; y=fft(x); k=(n-1)*f0; subplot(2,1,2); plot(k,abs(y(n))); title('FFT變換后男音的頻譜特性');xlabel('頻率/omega');ylabel('幅值/A');grid; %信號(hào)的時(shí)域采樣點(diǎn) %采樣間隔 %取信號(hào)的一半 %對(duì)信號(hào)做N點(diǎn)FFT變換 %頻域采樣點(diǎn) %把畫(huà)圖區(qū)域劃分為2行1列，指定第一個(gè)圖%繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 %給圖形加注標(biāo)簽說(shuō)明 %添加網(wǎng)格 %讀取聲音 %計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) %信號(hào)的時(shí)域采樣點(diǎn) %采樣間隔 %取信號(hào)的一半 %對(duì)信號(hào)做N點(diǎn)FFT變換 %頻域采樣點(diǎn) %把畫(huà)圖區(qū)域劃分為2行1列，指定第二個(gè)圖%繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 %給圖形加注標(biāo)簽說(shuō)明 %添加網(wǎng)格 axis([0 6000 0 300])； %改變橫縱坐標(biāo)便于比較頻譜圖 FFT變換后女音的頻譜特性300200幅值/A***00頻率/?FFT變換后男音的頻譜特性***200幅值/A***00頻率/?400050006000 結(jié)果分析： 通過(guò)比較上面女音頻譜圖和男音頻譜圖可知，男音的頻譜集中在低頻部分，高頻成分底，譜線較平滑，聲音聽(tīng)起來(lái)低沉。5.8有背景噪聲的信號(hào)分析 從硬盤(pán)中把一段噪聲（頻譜能量集中在某個(gè)小范圍內(nèi)）疊加到語(yǔ)音信號(hào)中，分析疊加前后信號(hào)頻譜的變化，設(shè)計(jì)一個(gè)合適的濾波器，能夠把該噪聲濾除； 程序如下： z=wavread('女音.wav',[1 24000]); %讀取聲音在1-24000之間 f=wavread('noise.wav',[1 24000]);x=z+f;wavplay(x);fs=11025;N=length(x);f0=fs/N; %采樣間隔 n=1:N; %取信號(hào)的一半 y=fft(x,N);%對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0; %頻域采樣點(diǎn) subplot(2,1,1); %把畫(huà)圖區(qū)域劃分為1行2列，指定第二個(gè)圖 plot(k,abs(y(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('加噪聲后聲音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格 %iir帶通濾波器的代碼： Ap=1; %通帶波紋系數(shù) Az=70; %最小阻帶衰減 wp=[0.2 0.7]; %歸一化通帶數(shù)字截止頻率 wz=[0.1 0.8]; %歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az); %估計(jì)契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn); %N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動(dòng) x1=filter(b,a,x); %對(duì)聲音濾波 wavplay(x1)； wavwrite(x1,'濾除噪音后女音.wav');N=length(x1);f0=fs/N; %采樣間隔 n=1:N; %取信號(hào)的一半 y1=fft(x1,N); %對(duì)信號(hào)做fs點(diǎn)FFT變換 subplot(2,1,2); %把畫(huà)圖區(qū)域劃分為1行2列，指定第二個(gè)圖 k=(n-1)*f0; %頻域采樣點(diǎn) plot(k,abs(y1(n))); %繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('濾除噪聲后聲音的頻譜特性'); %給圖形加注標(biāo)簽說(shuō)明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格 加噪聲后聲音的頻譜特性3000幅值/A***0008000頻率/Hz濾除噪聲后聲音的頻譜特性***030幅值/A***000頻率/Hz80001000012000 結(jié)果分析 觀察加噪聲后聲音的頻譜圖可知，噪音頻率主要在4000Hz處，只要我們?cè)O(shè)計(jì)一個(gè)，濾波器濾除大概在4000Hz的頻譜即可，回放濾波后的語(yǔ)音信號(hào)，可證噪音基本濾除。 六．心得體會(huì)： 通過(guò)這次課程設(shè)計(jì)，讓我對(duì)MATLAB的基本應(yīng)用有了更深的了解，還有數(shù)字信號(hào)處理在MATLAB中的一些函數(shù)的用法。通過(guò)理論推導(dǎo)得出相應(yīng)結(jié)論，并利用MATLAB作為工具進(jìn)行實(shí)現(xiàn)，從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識(shí)，提高對(duì)所學(xué)知識(shí)的綜合應(yīng)用能力，并從實(shí)踐上初步實(shí)現(xiàn)對(duì)數(shù)字信號(hào)的處理。 在這次實(shí)驗(yàn)中，也遇到了很多問(wèn)題，比如畫(huà)信號(hào)頻譜的分貝圖時(shí)（20*log10(abs(y))）指數(shù)為零時(shí)的處理。濾波器的設(shè)計(jì)也花了好大的功夫，剛開(kāi)始不會(huì)設(shè)計(jì)參數(shù)，一頭霧水，通過(guò)同學(xué)的指導(dǎo)和討論，得知通過(guò)觀察信號(hào)的頻譜圖，看噪音頻率集中在那一部分，設(shè)計(jì)濾波器把其濾除即可。可反復(fù)設(shè)置參數(shù)直到濾波后語(yǔ)音信號(hào)的效果好為止。 七．參考文獻(xiàn)： （1）《MATLAB LabVIEW SystemView》翁劍楓 葉志前 編著, 機(jī)械工業(yè)出版社； （2）《MATLAB及在電子信息課程中的應(yīng)用》陳懷琛 吳大正 高西全編著，電子工業(yè)出版社； （3）《MATLAB在數(shù)字信號(hào)處理中的應(yīng)用》（弟2版）薛年喜 編著，清華大學(xué)出版社； （4）《MATLAB擴(kuò)展編程》何強(qiáng) 何英 編著，清華大學(xué)出版社；（5）《MATLAB7簡(jiǎn)明教程》吳清 曹輝林 編著，清華大學(xué)出版社；（6）MATLAB5.3精要.編程及高級(jí)應(yīng)用》程衛(wèi)國(guó) 馮峰 王雪梅 劉藝 編著，機(jī)械工程出版社。 邯 鄲 學(xué) 院 講 稿 2010 ～2011 學(xué)年 第 一 學(xué)期 分院(系、部)： 信息工程學(xué)院 教 研 室： 電子信息工程 課 程 名 稱： 數(shù)字信號(hào)處理 授 課 班 級(jí)： 07級(jí)電子信息工程 主 講 教 師： 王苗苗 職 稱： 助教（研究生） 使 用 教 材： 《數(shù)字信號(hào)處理》 制 作 系 統(tǒng)： Word2003 邯鄲學(xué)院制 實(shí)驗(yàn)一..Matlab仿真軟件介紹 一、實(shí)驗(yàn)?zāi)康?/p> 熟悉Matlab仿真軟件 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、學(xué)習(xí)Matlab仿真軟件的安裝 2、熟悉Matlab仿真軟件的操作環(huán)境 3、直接在Matlab仿真軟件的命令窗口實(shí)現(xiàn)數(shù)值計(jì)算 4、編寫(xiě)M文件 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 1、熟悉Matlab仿真軟件 2、參閱Matlab及在電子信息類課程中的應(yīng)用（第2版）唐向宏 電子工業(yè)出版社 實(shí)驗(yàn)二 離散信號(hào)和系統(tǒng)分析的Matlab實(shí)現(xiàn) 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)離散信號(hào)和系統(tǒng)分析 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、利用Matlab產(chǎn)生離散信號(hào) 2、利用Matlab計(jì)算離散卷積 3、利用Matlab求解離散LTI系統(tǒng)響應(yīng) 4、利用Matlab計(jì)算DTFT 5、利用Matlab實(shí)現(xiàn)部分分式法 6、利用Matlab計(jì)算系統(tǒng)的零極點(diǎn) 7、利用Matlab進(jìn)行簡(jiǎn)單數(shù)字濾波器設(shè)計(jì) 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)三 利用Matlab實(shí)現(xiàn)信號(hào)DFT的計(jì)算 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)信號(hào)DFT的計(jì)算 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、利用Matlab計(jì)算信號(hào)的DFT 2、利用Matlab實(shí)現(xiàn)由DFT計(jì)算線性卷積 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)四 利用Matlab實(shí)現(xiàn)濾波器設(shè)計(jì) 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)實(shí)現(xiàn)濾波器設(shè)計(jì) 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、利用Matlab實(shí)現(xiàn)模擬低通濾波器的設(shè)計(jì) 2、利用Matlab實(shí)現(xiàn)模擬域頻率變換 3、利用Matlab實(shí)現(xiàn)脈沖響應(yīng)不變法 4、利用Matlab實(shí)現(xiàn)雙線性變換法 5、利用Matlab實(shí)現(xiàn)數(shù)字濾波器設(shè)計(jì) 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)五 利用Matlab實(shí)現(xiàn)FIR濾波器設(shè)計(jì) 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)實(shí)現(xiàn)濾波器設(shè)計(jì) 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、利用Matlab實(shí)現(xiàn)窗函數(shù)法 2、利用Matlab實(shí)現(xiàn)頻率取樣法 3、利用Matlab實(shí)現(xiàn)優(yōu)化設(shè)計(jì) 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)六..隨機(jī)信號(hào)功率譜估計(jì)的Matlab實(shí)現(xiàn) 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)實(shí)現(xiàn)濾波器設(shè)計(jì) 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、利用Matlab實(shí)現(xiàn)隨機(jī)序列 2、利用Matlab計(jì)算相關(guān)函數(shù)的估計(jì) 3、利用Matlab進(jìn)行非參數(shù)功率譜估計(jì) 4、利用Matlab進(jìn)行AR模型功率譜估計(jì) 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)七..數(shù)字濾波器結(jié)構(gòu)的Matlab實(shí)現(xiàn) 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)實(shí)現(xiàn)濾波器設(shè)計(jì) 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、利用Matlab實(shí)現(xiàn)數(shù)字濾波器直接型設(shè)計(jì) 2、利用Matlab實(shí)現(xiàn)數(shù)字濾波器級(jí)聯(lián)設(shè)計(jì) 3、利用Matlab實(shí)現(xiàn)數(shù)字濾波器并聯(lián)型設(shè)計(jì) 4、利用Matlab實(shí)現(xiàn)數(shù)字濾波器格型設(shè)計(jì) 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)八....利用Matlab實(shí)現(xiàn)信號(hào)小波分析 一、實(shí)驗(yàn)?zāi)康?/p> 1、Matlab實(shí)現(xiàn)實(shí)現(xiàn)濾波器設(shè)計(jì) 2、進(jìn)一步熟悉Matlab軟件操作 二、實(shí)驗(yàn)設(shè)備和元器件 含Matlab仿真軟件的計(jì)算機(jī) 三、實(shí)驗(yàn)內(nèi)容和步驟 1、小波測(cè)試信號(hào) 2、分解與重構(gòu)濾波器組 3、離散小波變換 4、離散小波反變換 5、基于小波的信號(hào)去噪 6、基于小波的信號(hào)壓縮 四、實(shí)驗(yàn)報(bào)告要求 按照《Matlab程序設(shè)計(jì)》模板提交實(shí)驗(yàn)報(bào)告 五、預(yù)習(xí)要求 預(yù)習(xí)課本上的相關(guān)內(nèi)容 實(shí)驗(yàn)五 FFT 算法的應(yīng)用 1.進(jìn)一步加深對(duì)離散信號(hào) DFT 的理解。 2.運(yùn)用其 FFT 算法解決實(shí)際問(wèn)題。 1.微機(jī)。 2.Matlab 編程環(huán)境。 1.熟悉 Matlab 的編程環(huán)境和編程語(yǔ)言。 2.學(xué)習(xí)教材 P213-227，P249-263，掌握快速傅里葉變換（FFT）的原理。 1.實(shí)驗(yàn)重點(diǎn)、難點(diǎn)、特點(diǎn) 快速傅里葉變換（FFT）的原理及應(yīng)用。難點(diǎn)在 FFT 的應(yīng)用以及 Matlab 編程中矩陣乘和數(shù)乘的應(yīng)用。 2.實(shí)驗(yàn)原理 一、實(shí)驗(yàn)?zāi)康?/p> 二、實(shí)驗(yàn)儀器設(shè)備 三、實(shí)驗(yàn)學(xué)時(shí) 學(xué)時(shí) 四、預(yù)習(xí)要求 五、實(shí)驗(yàn)特點(diǎn)及實(shí)驗(yàn)原理簡(jiǎn)介 1.已知 2N 點(diǎn)實(shí)數(shù)序列 N=64。理論計(jì)算 X(k)=DFT[x(n)]2N。 六、實(shí)驗(yàn)內(nèi)容及步驟 N=64;n=0:1:N-1;n2=0:1:2*N-1;k=0:1:2*N-1;xn1=cos(2*pi/N*7*2*n)+1/2*cos(2*pi/N*19*2*n);xn2=cos(2*pi/N*7*(2*n+1))+1/2*cos(2*pi/N*19*(2*n+1));wn=xn1+i*xn2;[wk,kk]=dft_my(wn,N);xk1=1/2*(wk+conj(wk(mod(-kk,N)+1)));xk2=1/2*(wk-conj(wk(mod(-kk,N)+1)))/i;XK1=fft(xn1,N);XK2=fft(xn2,N);%subplot(2,2,1);%stem(kk,abs(xk1));grid on;%subplot(2,2,3);%stem(kk,abs(XK1));grid on;%subplot(2,2,2);%stem(kk,abs(xk2));grid on;%subplot(2,2,4);%stem(kk,abs(XK2));grid on; xk=[xk1+(exp(-1i*pi/N).^kk).*xk2,xk1-(exp(-1i*pi/N).^kk).*xk2];subplot(2,1,1);stem(k,abs(xk),'r');grid on;xlabel('k');ylabel('2*NFFT');title('n點(diǎn)DFT完成2N點(diǎn)FFT');xn=cos(2*pi/N*7*n2)+1/2*cos(2*pi/N*19*n2);xkk=fft(xn,2*N);subplot(2,1,2);stem(k,abs(xkk),'g');grid on;xlabel('k');ylabel('2*NFFT');title('Matlab2N點(diǎn)FFT'); 七、問(wèn)題思考 1.兩次編程計(jì)算與理論計(jì)算相比較，結(jié)果一致嗎？說(shuō)明原因。 2.兩次編程計(jì)算在編程方面、運(yùn)算量上的比較。 八、心得體會(huì) 北京信息科技大學(xué) 數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告 題 目： 數(shù)字信號(hào)處理課程設(shè)計(jì)實(shí)驗(yàn) 學(xué) 院： 信息與通信工程學(xué)院 專 業(yè)： 通信工程專業(yè) 姓 名： 班 級(jí)： 學(xué) 號(hào)： 指導(dǎo)老師： 實(shí)驗(yàn)?zāi)康?/p> 1、熟悉IIR數(shù)字濾波器的設(shè)計(jì)原理與方法。 2、掌握數(shù)字濾波器的計(jì)算機(jī)軟件實(shí)現(xiàn)方法。 3、通過(guò)觀察對(duì)實(shí)際心電圖信號(hào)的濾波作用，學(xué)習(xí)數(shù)字濾波器在實(shí)際中的應(yīng)用。 實(shí)驗(yàn)儀器及材料 計(jì)算機(jī)，MATLAB軟件 實(shí)驗(yàn)內(nèi)容及要求 1．設(shè)計(jì)巴特沃斯低通數(shù)字濾波器對(duì)人體心電信號(hào)進(jìn)行濾波 (1)人體心電圖信號(hào)在測(cè)量過(guò)程中會(huì)受到工業(yè)高頻干擾，所以必須經(jīng)過(guò)低通濾波處理，才能作為判斷心臟功能的有用信息。以下為一個(gè)實(shí)際心電圖信號(hào)采樣序列x(n)，其中存在高頻干擾，抽樣周期Ts=1秒。在實(shí)驗(yàn)中，以x(n)作為輸入序列，濾除其中干擾成分。 x(n)=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0] 對(duì)序列x(n)用FFT做頻譜分析，生成x(n)的頻譜圖。 (2)設(shè)計(jì)一個(gè)巴特沃斯低通IIR數(shù)字濾波器H(z)。 設(shè)計(jì)指標(biāo)參數(shù)為：在通帶內(nèi)頻率低于0.2π時(shí)，最大衰減小于1dB； 在阻帶內(nèi) [0.3π, π]頻率區(qū)間上，最小衰減大于15dB。 j?|H(e)|。寫(xiě)出數(shù)字濾波器H(z)的表達(dá)式，畫(huà)出濾波器的幅頻響應(yīng)曲線 (3)用所設(shè)計(jì)的濾波器對(duì)實(shí)際心電圖信號(hào)采樣序列x(n)進(jìn)行濾波處理，編寫(xiě)程序，求濾波后的序列y(n)，并分別畫(huà)出濾波前后的心電圖信號(hào)波形圖和頻譜圖。 y(n)= [0,0,0,0, 0,0,0,0,-0.14025,0.40279,-0.56085 ,0.33328,0.023981,-0.18809,0.11843,-0.1038,0.11576,-0.1225,0.099815 ,-0.13769 ,0.095249,-0.0070273,0.018867,0.090543,-0.11257,-0.070884 ,0.17676,-0.55407,0.24813,-0.34732,-0.30428,0.59426,-0.29574,-0.063869,0.34018,-0.73334,1.0293,-0.57107,-0.2461,0.83605,-0.83026,0.45459,0.011551,-0.25667,0.23896,-0.17361,0.20829,-0.28417,0.28765 ,-0.2035,0.02865,0.066164,0.077916,-0.36052, 0.53517,-0.5571] 源程序 clear all，clc x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];%未經(jīng)濾波的心電圖信號(hào) L=length(x);l=0:L-1; y=fft(x,L);Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=15; [N,Wn] = buttord(Wp,Ws,Rp,Rs,'s');[b,a] = butter(N,Wn,'s');[numa,dena]=impinvar(b,a,1);w=linspace(0,pi,1024);h=freqz(numa,dena,w);norm=max(abs(h));numa=numa/norm;[z,p]=tf2zp(b,a);figure(1) plot(w,20*log10(abs(h)/norm));grid; xlabel('數(shù)字頻率');ylabel('幅度響應(yīng)dB');figure(2)plot(w,abs(h));grid; xlabel('數(shù)字頻率'); ylabel('幅度響應(yīng)|H(e^(jw))|');figure(3)zplane(z,p);xx=filter(b,a,x);yy=fft(xx,L);figure(4)subplot(2,1,1)stem(l,x); title('未經(jīng)濾波的心電圖信號(hào)');xlabel('n');subplot(2,1,2)stem(l,xx); title('經(jīng)濾波之后的心電圖信號(hào)');xlabel('n');figure(5)subplot(2,1,1)plot(l,abs(y)); title('未經(jīng)濾波的心電圖信號(hào)的頻譜');subplot(2,1,2)plot(l,abs(yy)); title('經(jīng)濾波處理的心電圖信號(hào)的頻譜'); 2．用help查看內(nèi)部函數(shù)cheb1ord.m及cheby1.m，了解調(diào)用格式。 編程設(shè)計(jì)教材習(xí)題6-2，求模擬濾波器Ha(s)的表達(dá)式。 源程序 close all clear all clc Wp=2*pi*3000;Rp=2;Ws=2*pi*12000;Rs=50;Fs=24000; w=linspace(0,pi,1024); [N,Wn]=cheb1ord(Wp,Ws,Rp,Rs,'s');e=sqrt(10^(Rp/10)-1);[b,a]=cheby1(N,e,Wn,'s')[numa,dena]=impinvar(b,a,Fs);h=freqz(numa,dena,w);norm=max(abs(h)); plot(w*Fs/pi,20*log10(abs(h)/norm))title('幅度響應(yīng)')xlabel('頻率(Hz)')ylabel('幅度(dB)')grid 3．模擬濾波器的數(shù)字化 用內(nèi)部函數(shù)impinvar及bilinear實(shí)現(xiàn)教材習(xí)題6-5，求數(shù)字濾波器H(z)的表達(dá)式。 源程序 close all clear all clc b1=[0 0 1];a1=[1 1 1];b2=[0 0 1];a2=[2 3 1]; [numa1,dena1]=impinvar(b1,a1,0.5)[numa2,dena2]=bilinear(b1,a1,0.5)[numa3,dena3]=impinvar(b2,a2,0.5) [numa4,dena4]=bilinear(b2,a2,0.5) 本實(shí)驗(yàn)所用的部分MATLAB函數(shù) ? L=length(x)：求序列x長(zhǎng)度。 ? y=fft(x,L)：將序列x(n)做L點(diǎn)快速傅立葉變換，結(jié)果賦給序列y(n)。? [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s')：計(jì)算模擬Butterworth濾波器的最小階次n和截止頻率為Wn。 ? [b,a] = butter(n,Wn,'s')：設(shè)計(jì)模擬截止頻率為Wn（rad/s）的n階 Butterworth低通濾波器，返回值為模擬濾波器的系數(shù)。 ? y=filter(b,a,x): 將序列x(n)通過(guò)濾波器濾波后生成序列y(n)，濾波器的分母多項(xiàng)式系數(shù)構(gòu)成a向量，分子多項(xiàng)式系數(shù)構(gòu)成b向量。 ? [BZ,AZ] = impinvar(B,A,Fs)：沖激響應(yīng)不變法，返回值為數(shù)字濾波器的系數(shù)。? [BZ,AZ] = bilinear(B,A,fs)：雙線性變換，返回值為數(shù)字濾波器的系數(shù)。? [H w]=freqz(b,a)：由濾波器分母多項(xiàng)式系數(shù)構(gòu)成的a向量和分子多項(xiàng)式系數(shù)構(gòu)成的b向量求系統(tǒng)頻響。 截圖 實(shí)驗(yàn)體會(huì) 這次的實(shí)驗(yàn)讓這讓我看到了理論與實(shí)踐相結(jié)合的優(yōu)勢(shì)與用處，讓我受益匪淺。我認(rèn)識(shí)到了自己理論知識(shí)的不足，也認(rèn)識(shí)到了我們學(xué)習(xí)的基礎(chǔ)知識(shí)究竟能運(yùn)用于什么領(lǐng)域，如何運(yùn)用。我們?cè)诶蠋煹哪托闹笇?dǎo)下調(diào)試電路，直到得到要求的效果。讓我們?cè)趯W(xué)習(xí)電路、信號(hào)等理論知識(shí)的同時(shí)，明白如何把這些應(yīng)用于實(shí)際。第二篇：MATLAB實(shí)現(xiàn)數(shù)字信號(hào)處理
第三篇：數(shù)字信號(hào)處理實(shí)驗(yàn)講稿
第四篇：數(shù)字信號(hào)處理實(shí)驗(yàn)5
第五篇：數(shù)字信號(hào)處理實(shí)驗(yàn)二
點(diǎn)擊咨詢 