第一篇：數(shù)字信號處理實驗-FFT的實現(xiàn)

實

驗

報

告

學生姓名：

學 號：

指導教師：

一、實驗室名稱：數(shù)字信號處理實驗室

二、實驗項目名稱：FFT的實現(xiàn)

三、實驗原理：

一．FFT算法思想：

1．DFT的定義：

對于有限長離散數(shù)字信號{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，稱為旋轉因子。

2．直接計算DFT的問題及FFT的基本思想：

由DFT的定義可以看出，在x[n]為復數(shù)序列的情況下，完全直接運算N點DFT需要（N-1）2次復數(shù)乘法和N（N-1）次加法。因此，對于一些相當大的N值（如1024）來說，直接計算它的DFT所作的計算量是很大的。

FFT的基本思想在于，將原有的N點序列分成兩個較短的序列，這些序列的DFT可以很簡單的組合起來得到原序列的DFT。例如，若N為偶數(shù)，將原有的N

22點序列分成兩個（N/2）點序列，那么計算N點DFT將只需要約[(N/2)·2]=N/2次復數(shù)乘法。即比直接計算少作一半乘法。因子（N/2）2表示直接計算（N/2）點DFT所需要的乘法次數(shù)，而乘數(shù)2代表必須完成兩個DFT。上述處理方法可以反復使用，即（N/2）點的DFT計算也可以化成兩個（N/4）點的DFT（假定N/2為偶數(shù)），從而又少作一半的乘法。這樣一級一級的劃分下去一直到最后就劃分成兩點的FFT運算的情況。

3．基2按時間抽取（DIT）的FFT算法思想：

設序列長度為N?2L，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。

將長度為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上式中利用了旋轉因子的可約性，即：WNN/?21NrkN?/21rk?WN/2。又令

rkX1[k]??r?0x[1r]W,/X2[k]?2?r?0x[r]WN2，則上式可以寫成： /2X[k]?X1[k]?WNX2[k](k=0,1,…,N/2-1)

k可以看出，X1[k],X2[k]分別為從X[k]中取出的N/2點偶數(shù)點和奇數(shù)點序列的N/2點DFT值，所以，一個N點序列的DFT可以用兩個N/2點序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了X[k]前N/2點的DFT值，還需要繼續(xù)利用X1[k],X2[k]表示X[k]的后半段本算法推導才完整。利用旋轉因子的周期性，有：WN/2?WN/2X1[N2N/2?1rkr(k?N/2)，則后半段的DFT值表達式：

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值表達式獲得，中間還利用到WN(N2?k)N?WN2Wk得到后半段的X[k]值表達式??W，k為：X[k]?X1[k]?WNkX2[k]（k=0,1,…,N/2-1）。

這樣，通過計算兩個N/2點序列x1[n],x2[n]的N/2點DFTX1[k],X2[k]，可以組合得到N點序列的DFT值X[k]，其組合過程如下圖所示：

X1[k] X1[k]?WNkX2[k]

X2[k] WNnk-1 X1[k]?WNkX2[k]

比如，一個N = 8點的FFT運算按照這種方法來計算FFT可以用下面的流程圖來表示：

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按頻率抽取（DIF）的FFT算法思想：

設序列長度為N?2L，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。

在把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因為W2N??1，則有WX[k]???(?1)，所以：

kkN/2?1?n?0[x[n]?(?1)x[n?N2]]?WN,k?0,1,...,N?1

nk按k的奇偶來討論，k為偶數(shù)時：

N/2?1X[2r]??n?0[x[n]?x[n?N2]]?WN,k?0,1,...,N?1 N22rnN/2?1k為奇數(shù)時：X[2r?1]?前面已經(jīng)推導過WNN/2?1?n?0[x[n]?x[n?]]?WN(2r?1)n,k?0,1,...,N?1

2rk?WN/2，所以上面的兩個等式可以寫為：

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通過上面的推導，X[k]的偶數(shù)點值X[2r]和奇數(shù)點值X[2r?1]分別可以由組合而成的N/2點的序列來求得，其中偶數(shù)點值X[2r]為輸入x[n]的前半段和后半段之和序列的N/2點DFT值，奇數(shù)點值X[2r?1]為輸入x[n]的前半段和后半段之差再與WN相乘序列的N/2點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

這樣，也可以用兩個N/2點DFT來組合成一個N點DFT，組合過程如下圖所示：

x[n] x[n]?x[n?N2]

x[n?N2]-1 WNn [x[n]?x[n?N2]]WNn

二．在FFT計算中使用到的MATLAB命令：

函數(shù)fft(x)可以計算R點序列的R點DFT值；而fft(x,N)則計算R點序列的N點DFT，若R>N，則直接截取R點DFT的前N點，若R

四、實驗目的：

離散傅氏變換（DFT）的目的是把信號由時域變換到頻域，從而可以在頻域分析處理信息，得到的結果再由逆DFT變換到時域。FFT是DFT的一種快速算法。在數(shù)字信號處理系統(tǒng)中，F(xiàn)FT作為一個非常重要的工具經(jīng)常使用，甚至成為DSP運算能力的一個考核因素。

本實驗通過直接計算DFT，利用FFT算法思想計算DFT，以及使用MATLAB函數(shù)中的FFT命令計算離散時間信號的頻譜，以加深對離散信號的DFT變換及FFT算法的理解。

五、實驗內(nèi)容：

a)計算實數(shù)序列x(n)?cos5?16n,0?n?256的256點DFT。

b)計算周期為1kHz的方波序列（占空比為50％，幅度取為＋/-512，采樣頻率為25kHz，取256點長度）256點DFT。

六、實驗器材（設備、元器件）：

安裝MATLAB軟件的PC機一臺，DSP實驗演示系統(tǒng)一套。

七、實驗步驟：

（1）先利用DFT定義式，編程直接計算2個要求序列的DFT值。

（2）利用MATLAB中提供的FFT函數(shù)，計算2個要求序列的DFT值。（3）(拓展要求)不改變序列的點數(shù)，僅改變DFT計算點數(shù)（如變?yōu)橛嬎?024點DFT值），觀察畫出來的頻譜與前面頻譜的差別，并解釋這種差別。通過這一步驟的分析，理解頻譜分辨力的概念，解釋如何提高頻譜分辨力。

（4）利用FFT的基本思想（基2－DIT或基2－DIF），自己編寫FFT計算函數(shù)，并用該函數(shù)計算要求序列的DFT值。并對前面3個結果進行對比。

（5）(拓展要求)嘗試對其他快速傅立葉變換算法（如Goertzel算法）進行MATLAB編程實現(xiàn)，并用它來計算要求的序列的DFT值。并與前面的結果進行對比。

（6）（拓展要求）在提供的DSP實驗板上演示要求的2種序列的FFT算法（基2－DIT），用示波器觀察實際計算出來的頻譜結果，并與理論結果對比。

八、實驗數(shù)據(jù)及結果分析：

程序：（1）對要求的2種序列直接進行DFT計算的程序

（2）對要求的2種序列進行基2－DIT和基2－DIF FFT算法程序（3）對要求的2種序列用MATLAB中提供的FFT函數(shù)進行計算的程序

結果：（1）對2種要求的序列直接進行DFT計算的頻域波形

（2）對2種要求的序列進行基2－DIT和基2－DIF FFT算法頻域波形（3）對2種要求的序列用MATLAB中提供的FFT函數(shù)計算的頻域波形。（4）（拓展要求）分析利用上面的方法畫出的信號頻譜與理論計算出來的頻譜之間的差異，并解釋這種差異。

（5）（拓展要求）保持序列點數(shù)不變，改變DFT計算點數(shù)（變?yōu)?024點），觀察頻譜的變化，并分析這種變化，由此討論如何提高頻譜分辨力的問題。

九、實驗結論：

十、總結及心得體會：

十一、對本實驗過程及方法、手段的改進建議：

第二篇：MATLAB實現(xiàn)數(shù)字信號處理

數(shù)字信號處理

說 明 書

目錄

一．摘要…………………………………3 二．課程設計目的………………………3 三．設計內(nèi)容……………………………3 四．設計原理……………………………4 4.1．語音信號的采集…………………………….4 4.2．濾波器……………………………………….4 4．21．IIR濾波器原理…………………………………….4 4.22．FIR濾波器原理………………………………………5 五．設計步驟……………………………6 5.1錄制女音………………………………………6 5.2采樣語音信號并畫出時域波形和頻譜圖……7 5.3采用雙線性變換法設計IIR濾波器…………10 5.4窗函數(shù)法設計FFR濾波器………………......12 5.5用IIR濾波器對信號進行濾波………………14 5.6用FIR濾波器對信號進行濾波………………16 5.7男女聲語音信號頻譜特點分析………………19 5.8有背景噪聲的信號分析………………………20 六．心得體會…………………………….22 七．參考文獻…………………………….23

一．摘要：

這次課程設計的主要目的是綜合運用本課程的理論知識進行頻譜分析以及濾波器設計，通過理論推導得出相應結論，并利用MATLAB或者ＤＳＰ開發(fā)系統(tǒng)作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。通過對聲音的采樣，將聲音采樣后的頻譜與濾波。

MATLAB全稱是Matrix Laboratory，是一種功能強大、效率高、交互性好的數(shù)值和可視化計算機高級語言，它將數(shù)值分析、矩陣運算、信號處理和圖形顯示有機地融合為一體，形成了一個極其方便、用戶界面友好的操作環(huán)境。經(jīng)過多年的發(fā)展，已經(jīng)發(fā)展成為一種功能全面的軟件，幾乎可以解決科學計算中所有問題。MATLAB軟件還提供了非常廣泛和靈活的用于處理數(shù)據(jù)集的數(shù)組運算功能。

在本次課程設計中，主要通過MATLAB來編程對語音信號處理與濾波，設計濾波器來處理數(shù)字信號并對其進行分析。

二．課程設計目的：

綜合運用本課程的理論知識進行頻譜分析以及濾波器設計，通過理論推導得出相應結論，并利用MATLAB作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。

三．設計內(nèi)容：

內(nèi)容：錄制一段個人自己的語音信號，并對錄制的信號進行采樣；畫出采樣后語音信號的時域波形和頻譜圖；給定濾波器的性能指標，采用窗函數(shù) 法和雙線性變換法設計濾波器，并畫出濾波器的頻率響應；然后用自己設計的濾波器對采集的信號進行濾波，畫出濾波后信號的時域波形和頻譜，并對濾波前后的信號進行對比，分析信號的變化；回放語音信號；換一個與你性別相異的人錄制同樣一段語音內(nèi)容，分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點；再錄制一段同樣長時間的背景噪聲疊加到你的語音信號中，分析疊加前后信號頻譜的變化，設計一個合適的濾波器，能夠把該噪聲濾除。

四．設計原理：

4.1．語音信號的采集

熟悉并掌握MATLAB中有關聲音（wave）錄制、播放、存儲和讀取的函數(shù)，在MATLAB環(huán)境中，有關聲音的函數(shù)有：

a：y=wavrecord(N,fs,Dtype);利用系統(tǒng)音頻輸入設備錄音，以fs為采樣頻率，默認值為11025，即以11025HZ進行采樣。Dtype為采樣數(shù)據(jù)的存儲格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度進行采樣，其它三種都是16位采樣結果轉換為指定的MATLAB數(shù)據(jù)；

b：wavplay(y,fs)；利用系統(tǒng)音頻輸出設備播放，以fs為播放頻率，播放語音信號y；

c：wavwrite((y,fs,wavfile);創(chuàng)建音頻文件； d：y=wavread(file);讀取音頻文件；

關于聲音的函數(shù)還有sound()；soundsc()；等。4.2濾波器： 4．21.IIR濾波器原理

沖激響應不變法是使數(shù)字濾波器在時域上模擬濾波器，但是它們的缺點是產(chǎn)生頻率響應的混疊失真，這是由于從s平面到z平面是多值的映射關系所造成的。

雙線性變換法是使數(shù)字濾波器的頻率響應與模擬濾波器的頻率響應相似的一種變換方法。為了克服多值映射這一缺點，我們首先把整個s平面壓縮變換到某一中介的s1平面的一條橫帶里，再通過變換關系將此橫帶變換到整個z平面上去，這樣就使得s平面與z平面是一一對應的關系，消除了多值變換性，也 就消除了頻譜混疊現(xiàn)象。

雙線性法設計IIR數(shù)字濾波器的步驟：

1）將數(shù)字濾波器的頻率指標{ ?k}由Wk＝（2/T）*tan(wk),轉換為模擬濾波器的頻率指標{?k}.2)由模擬濾波器的指標設計H(s).3)由H(s)轉換為H(z)21?z?1H(z)?H(s)s?T1?z?1

4.22.FIR濾波器原理

FIR濾波器與IIR濾波器特點不同，IIR濾波器的相位是非線性的，若需線性相位則要采用全通網(wǎng)絡進行相位校正。而有限長單位沖激響應（FIR）數(shù)字濾波器就可以做成具有嚴格的線性相位，同時又可以具有任意的幅度特性。

由于FIR系統(tǒng)的沖激響應就是其系統(tǒng)函數(shù)各次項的系數(shù)，所以設計FIR濾波器的方法之一可以從時域出發(fā)，截取有限長的一段沖激響應作為H（z）的系數(shù)，沖激響應長度N就是系統(tǒng)函數(shù)H（z）的階數(shù)。只要N足夠長，截取的方法合理，總能滿足頻域的要求。這種時域設計、頻域檢驗的方法一般要反復幾個回合，不像IIR DF設計靠解析公式一次計算成功。給出的理想濾波器頻率響應是，它是w的周期函數(shù)，周期

由傅立葉反變換導出，即

hd(n)?1Hd(ejw)ejwndw?2?，再將hd(n)與窗函數(shù)，因此可展開成傅氏級數(shù)w(n)相乘就可以得到h(n)。、的計算可采用傅氏變換的現(xiàn)成公式和程序，窗函數(shù)也是現(xiàn)成的。但整個設計過程不能一次完成，因為窗口類型和大小的選擇沒有解析公式可一次算，整個設計可用計算機編程來做。

窗函數(shù)的傅式變換W(ejω)的主瓣決定了H(ejω)過渡帶寬。W(ejω)的旁瓣大小和多少決定了H(ejω)在通帶和阻帶范圍內(nèi)波動幅度，常用的幾種窗函數(shù)有：

矩形窗

w(n)=RN(n)；

Hanning窗

；

Hamming窗

；

Blackmen窗

；

Kaiser窗。

式中Io(x)為零階貝塞爾函數(shù)。

五．設計步驟：

5.1錄制女音：

利用MATLAB中的函數(shù)錄制聲音。function nvyin()fs=11025;

%采樣頻率

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

開始錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double');

%錄制聲音3秒

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

錄音結束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs);

%播放錄音

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放錄音結束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音');

%聲音的存儲

5.2采樣語音信號并畫出時域波形和頻譜圖

讀取語音信號，畫出其時域波形和頻譜圖，與截取后的語音信號的時域波形和頻譜圖比較，觀察其變化。程序如下：

[x,fs,bits]=wavread('女音.wav');

%讀取聲音

N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半 figure(1);subplot(2,2,1);

%把畫圖區(qū)域劃分為2行2列，指定第一個圖 plot(t, x);

%畫出聲音采樣后的時域波形 title('原女音信號的時域波形');

%給圖形加注標簽說明 xlabel('時間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,2,3);

%把畫圖區(qū)域劃分為2行2列，指定第三個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('FFT變換后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%添加網(wǎng)格

subplot(2,2,4);

%把畫圖區(qū)域劃分為2行2列，指定第四個圖 if y~=0

%判斷指數(shù)是否為0

plot(k,20*log10(abs(y(n))));

%畫信號頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid;

%添加網(wǎng)格

%實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象的解決 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]);

%截取語音信號 N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

figure(2);subplot(2,2,1);

%把畫圖區(qū)域劃分為2行2列，指定第一個圖 plot(t,x1);

%畫出聲音采樣后的時域波形 title('截取后女音信號的時域波形');

%給圖形加注標簽說明 xlabel('時間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y1=fft(x1);

%對信號做N點FFT變換

subplot(2,2,3);

%把畫圖區(qū)域劃分為2行2列，指定第三個圖 k=(n-1)*f0;

%頻域采樣點

plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('FFT變換后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%添加網(wǎng)格

subplot(2,2,4);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 if y1~=0

%判斷指數(shù)是否為0

plot(k,20*log10(abs(y1(n))));

%畫信號頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid;

%添加網(wǎng)格

原女音信號的時域波形10.5A/幅0振-0.5-10123時間/tFFT變換后聲音的頻譜特性FFT變換后聲音的頻譜特性30050A200貝/值分/幅0幅100振00200040006000-***頻率/HzHz 截取后女音信號的時域波形10.5振幅/A0-0.5-10123FFT變換后聲音的頻譜特性50時間/tFFT變換后聲音的頻譜特性300200振幅/分貝幅值/A01000020004000頻率/Hz6000-5002000Hz40006000

結果分析：

由原女音信號的時域波形可知錄取開始時實際發(fā)出聲音大概落后3500個采樣點，我們把前3500點去除即可解決實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象。由原女音的的頻譜圖和截取后聲音的頻譜圖可看出，對聲音的截取并不會影響它們頻譜分布。

5.3采用雙線性變換法設計IIR濾波器：

人的聲音頻率一般在（1~~4）kHZ之間，則我們只需要設計一個帶通濾波器即可濾去聲音頻帶以外的無用噪聲，得到比較清晰的聲音。根據(jù)聲音的頻譜圖分析，設計一個帶通濾波器性能指標如下：

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);

%估計契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化

截 %止頻率，Ap通帶波動

[h,w]=freqz(b,a);

%求數(shù)字濾波器的復頻率響應 figure(1);subplot(2,1,1);plot(w/pi,abs(h));

%繪制數(shù)字濾波器的頻譜圖 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應');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('契比雪夫Ⅰ型帶通濾波器的幅頻響應');契比雪夫Ⅰ型帶通濾波器的幅頻響應1振幅（幅值）0.5000.10.20.50.60.70.8?/?契比雪夫Ⅰ型帶通濾波器的幅頻響應0.30.40.910振幅(分貝)-200-400-60000.10.20.30.40.5?/?0.60.70.80.91

5.4窗函數(shù)法設計FFR濾波器

線性相位FIR濾波器通常采用窗函數(shù)法設計。窗函數(shù)法設計FIR濾波器的基本思想是：根據(jù)給定的濾波器技術指標，選擇濾波器長度N和窗函數(shù)ω（n），使其具有最窄寬度的主瓣和最小的旁瓣。其核心是從給定的頻率特性，通過加窗確定有限長單位脈沖響應序列h(n)。工程中常用的窗函數(shù)共有6種，即矩形窗、巴特利特（Bartlett）窗、漢寧（Hanning）窗、漢明（Hamming）窗、布萊克曼（Blackman）窗和凱澤（Kaiser）窗。

這次設計我采用的是布萊克曼來設計給定數(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));

%---取兩個過渡帶中的小者 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';

%點乘

n=0:N-1;

%階數(shù) subplot(2,2,1);stem(n,wdWindow);

%繪制布拉克曼窗時域波形 xlabel('時間');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1);

%求濾波器頻率響應 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)));

%畫濾波器頻譜的分貝圖 end xlabel('omega/');ylabel('振幅/分貝');title('FIR帶通濾波器幅頻特性');grid;

%添加網(wǎng)格 %---ideallp（）函數(shù)（非系統(tǒng)自有函數(shù)）在系統(tǒng)安裝目錄的WORK子目錄ideallp.m function hd = ideallp(wc,N);% 理想低通濾波器的脈沖響應子程序 % hd = 點0 到 N-1之間的理想脈沖響應 % wc = 截止頻率(弧度)% N = 理想濾波器的長度

tao =(N-1)/2;

% 理想脈沖響應的對稱中心位置 n = [0:(N-1)];

% 設定脈沖響應長度 m = n-tao + eps;

% 加一個小數(shù)以避免零作除數(shù)

hd = sin(wc*m)./(pi*m);

% 理想脈沖響應

布拉克曼窗1振幅0.500406080時間FIR帶通濾波器幅頻特性500振幅/分貝20FIR帶通濾波器幅頻特性1.51振幅-50-100-15000.5?/10.5000.5?/?1

5.5用IIR濾波器對信號進行濾波

用自己設計的IIR濾波器分別對采集的信號進行濾波，在Matlab中，IIR濾波器利用函數(shù)filter對信號進行濾波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半 y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('濾波前女音的頻譜特性');

%給圖形加注標簽說明 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);

%估計契比雪夫I型濾波器階數(shù)

[b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動 x1=filter(b,a,x);

%對聲音濾波 wavplay(x1)wavwrite(x1,'IIR濾波后女音.wav');N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x1);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,2);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('l濾波后女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

濾波前女音的頻譜特性300幅值/A***030004000頻率/Hz濾波后女音的頻譜特性500060006040幅值/A***0頻率/Hz400050006000

結果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3400Hz的頻譜成分。回放語音信號，由于低頻和高頻成分被濾除，聲音變得較低沉。

5.6用FIR濾波器對信號進行濾波

用自己設計的FIR濾波器分別對采集的信號進行濾波，在Matlab中,FIR濾波器利用函數(shù)fftfilt對信號進行濾波 程序如下：

[x,fs,bits]=wavread('女音.wav');N=length(x);

%計數(shù)讀取信號的點數(shù)

t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('濾波前女音的頻譜特性');

%給圖形加注標簽說明 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));

%---取兩個過渡帶中的小者 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);

%對聲音濾波 wavplay(x1)wavwrite(x1,'FIR濾波后女音.wav');N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x1);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,2);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('l濾波后女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

濾波前女音的頻譜特性300200幅值/A***004000頻率/?l濾波后女音的頻譜特性500060006040幅值/A20005001000***03000頻率/Hz***0

結果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3500Hz的頻譜成分。和用IIR濾波器濾波一樣，回放語音信號，由于低頻和高頻成分被濾除，聲音變得較低沉。5.7男女聲語音信號頻譜特點分析

換一個男音錄制同樣一段語音內(nèi)容，分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點。程序如下：

[x,fs,bits]=wavread('女音.wav');

%讀取聲音

N=length(x);

%計數(shù)讀取信號的點數(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;

%信號的時域采樣點

%采樣間隔

%取信號的一半

%對信號做N點FFT變換

%頻域采樣點

%把畫圖區(qū)域劃分為2行1列，指定第一個圖%繪制原始語音信號的幅頻響應圖

%給圖形加注標簽說明

%添加網(wǎng)格

%讀取聲音

%計數(shù)讀取信號的點數(shù)

%信號的時域采樣點

%采樣間隔

%取信號的一半

%對信號做N點FFT變換

%頻域采樣點

%把畫圖區(qū)域劃分為2行1列，指定第二個圖%繪制原始語音信號的幅頻響應圖

%給圖形加注標簽說明

%添加網(wǎng)格

axis([0 6000 0 300])；

%改變橫縱坐標便于比較頻譜圖

FFT變換后女音的頻譜特性300200幅值/A***00頻率/?FFT變換后男音的頻譜特性***200幅值/A***00頻率/?400050006000

結果分析：

通過比較上面女音頻譜圖和男音頻譜圖可知，男音的頻譜集中在低頻部分，高頻成分底，譜線較平滑，聲音聽起來低沉。5.8有背景噪聲的信號分析

從硬盤中把一段噪聲（頻譜能量集中在某個小范圍內(nèi)）疊加到語音信號中，分析疊加前后信號頻譜的變化，設計一個合適的濾波器，能夠把該噪聲濾除； 程序如下：

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;

%取信號的一半 y=fft(x,N);%對信號做N點FFT變換

k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('加噪聲后聲音的頻譜特性');

%給圖形加注標簽說明 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);

%估計契比雪夫I型濾波器階數(shù)

[b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動 x1=filter(b,a,x);

%對聲音濾波 wavplay(x1)；

wavwrite(x1,'濾除噪音后女音.wav');N=length(x1);f0=fs/N;

%采樣間隔 n=1:N;

%取信號的一半

y1=fft(x1,N);

%對信號做fs點FFT變換

subplot(2,1,2);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 k=(n-1)*f0;

%頻域采樣點

plot(k,abs(y1(n)));

%繪制原始語音信號的幅頻響應圖 title('濾除噪聲后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格

加噪聲后聲音的頻譜特性3000幅值/A***0008000頻率/Hz濾除噪聲后聲音的頻譜特性***030幅值/A***000頻率/Hz80001000012000

結果分析

觀察加噪聲后聲音的頻譜圖可知，噪音頻率主要在4000Hz處，只要我們設計一個，濾波器濾除大概在4000Hz的頻譜即可，回放濾波后的語音信號，可證噪音基本濾除。

六．心得體會：

通過這次課程設計，讓我對MATLAB的基本應用有了更深的了解，還有數(shù)字信號處理在MATLAB中的一些函數(shù)的用法。通過理論推導得出相應結論，并利用MATLAB作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。

在這次實驗中，也遇到了很多問題，比如畫信號頻譜的分貝圖時（20*log10(abs(y))）指數(shù)為零時的處理。濾波器的設計也花了好大的功夫，剛開始不會設計參數(shù)，一頭霧水，通過同學的指導和討論，得知通過觀察信號的頻譜圖，看噪音頻率集中在那一部分，設計濾波器把其濾除即可。可反復設置參數(shù)直到濾波后語音信號的效果好為止。

七．參考文獻：

（1）《MATLAB LabVIEW SystemView》翁劍楓 葉志前 編著, 機械工業(yè)出版社；

（2）《MATLAB及在電子信息課程中的應用》陳懷琛 吳大正 高西全編著，電子工業(yè)出版社；

（3）《MATLAB在數(shù)字信號處理中的應用》（弟2版）薛年喜 編著，清華大學出版社；

（4）《MATLAB擴展編程》何強 何英

編著，清華大學出版社；（5）《MATLAB7簡明教程》吳清 曹輝林 編著，清華大學出版社；（6）MATLAB5.3精要.編程及高級應用》程衛(wèi)國 馮峰 王雪梅 劉藝 編著，機械工程出版社。

第三篇：數(shù)字信號處理實驗講稿

邯 鄲 學 院

講 稿

2010 ～2011 學年 第 一 學期

分院(系、部)： 信息工程學院 教 研 室： 電子信息工程 課 程 名 稱： 數(shù)字信號處理

授 課 班 級： 07級電子信息工程

主 講 教 師： 王苗苗 職

稱：

助教（研究生）

使 用 教 材： 《數(shù)字信號處理》

制 作 系 統(tǒng)：

Word2003

邯鄲學院制

實驗一..Matlab仿真軟件介紹

一、實驗目的

熟悉Matlab仿真軟件

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、學習Matlab仿真軟件的安裝

2、熟悉Matlab仿真軟件的操作環(huán)境

3、直接在Matlab仿真軟件的命令窗口實現(xiàn)數(shù)值計算

4、編寫M文件

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

1、熟悉Matlab仿真軟件

2、參閱Matlab及在電子信息類課程中的應用（第2版）唐向宏 電子工業(yè)出版社

實驗二 離散信號和系統(tǒng)分析的Matlab實現(xiàn)

一、實驗目的

1、Matlab實現(xiàn)離散信號和系統(tǒng)分析

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、利用Matlab產(chǎn)生離散信號

2、利用Matlab計算離散卷積

3、利用Matlab求解離散LTI系統(tǒng)響應

4、利用Matlab計算DTFT

5、利用Matlab實現(xiàn)部分分式法

6、利用Matlab計算系統(tǒng)的零極點

7、利用Matlab進行簡單數(shù)字濾波器設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

實驗三 利用Matlab實現(xiàn)信號DFT的計算

一、實驗目的

1、Matlab實現(xiàn)信號DFT的計算

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、利用Matlab計算信號的DFT

2、利用Matlab實現(xiàn)由DFT計算線性卷積

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

實驗四 利用Matlab實現(xiàn)濾波器設計

一、實驗目的

1、Matlab實現(xiàn)實現(xiàn)濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、利用Matlab實現(xiàn)模擬低通濾波器的設計

2、利用Matlab實現(xiàn)模擬域頻率變換

3、利用Matlab實現(xiàn)脈沖響應不變法

4、利用Matlab實現(xiàn)雙線性變換法

5、利用Matlab實現(xiàn)數(shù)字濾波器設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

實驗五 利用Matlab實現(xiàn)FIR濾波器設計

一、實驗目的

1、Matlab實現(xiàn)實現(xiàn)濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、利用Matlab實現(xiàn)窗函數(shù)法

2、利用Matlab實現(xiàn)頻率取樣法

3、利用Matlab實現(xiàn)優(yōu)化設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

實驗六..隨機信號功率譜估計的Matlab實現(xiàn)

一、實驗目的

1、Matlab實現(xiàn)實現(xiàn)濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、利用Matlab實現(xiàn)隨機序列

2、利用Matlab計算相關函數(shù)的估計

3、利用Matlab進行非參數(shù)功率譜估計

4、利用Matlab進行AR模型功率譜估計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

實驗七..數(shù)字濾波器結構的Matlab實現(xiàn)

一、實驗目的

1、Matlab實現(xiàn)實現(xiàn)濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、利用Matlab實現(xiàn)數(shù)字濾波器直接型設計

2、利用Matlab實現(xiàn)數(shù)字濾波器級聯(lián)設計

3、利用Matlab實現(xiàn)數(shù)字濾波器并聯(lián)型設計

4、利用Matlab實現(xiàn)數(shù)字濾波器格型設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

實驗八....利用Matlab實現(xiàn)信號小波分析

一、實驗目的

1、Matlab實現(xiàn)實現(xiàn)濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內(nèi)容和步驟

1、小波測試信號

2、分解與重構濾波器組

3、離散小波變換

4、離散小波反變換

5、基于小波的信號去噪

6、基于小波的信號壓縮

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內(nèi)容

第四篇：數(shù)字信號處理實驗5

實驗五 FFT 算法的應用

1.進一步加深對離散信號 DFT 的理解。

2.運用其 FFT 算法解決實際問題。

1.微機。

2.Matlab 編程環(huán)境。

1.熟悉 Matlab 的編程環(huán)境和編程語言。

2.學習教材 P213-227，P249-263，掌握快速傅里葉變換（FFT）的原理。

1.實驗重點、難點、特點

快速傅里葉變換（FFT）的原理及應用。難點在 FFT 的應用以及 Matlab 編程中矩陣乘和數(shù)乘的應用。

2.實驗原理

一、實驗目的

二、實驗儀器設備

三、實驗學時 學時

四、預習要求

五、實驗特點及實驗原理簡介

1.已知 2N 點實數(shù)序列

N=64。理論計算 X(k)=DFT[x(n)]2N。

六、實驗內(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點DFT完成2N點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點FFT');

七、問題思考

1.兩次編程計算與理論計算相比較，結果一致嗎？說明原因。

2.兩次編程計算在編程方面、運算量上的比較。

八、心得體會

第五篇：數(shù)字信號處理實驗二

北京信息科技大學

數(shù)字信號處理實驗報告

題 目：

數(shù)字信號處理課程設計實驗

學 院： 信息與通信工程學院 專 業(yè)： 通信工程專業(yè) 姓 名： 班 級：

學 號： 指導老師：

實驗目的

1、熟悉IIR數(shù)字濾波器的設計原理與方法。

2、掌握數(shù)字濾波器的計算機軟件實現(xiàn)方法。

3、通過觀察對實際心電圖信號的濾波作用，學習數(shù)字濾波器在實際中的應用。

實驗儀器及材料

計算機，MATLAB軟件

實驗內(nèi)容及要求

1．設計巴特沃斯低通數(shù)字濾波器對人體心電信號進行濾波

(1)人體心電圖信號在測量過程中會受到工業(yè)高頻干擾，所以必須經(jīng)過低通濾波處理，才能作為判斷心臟功能的有用信息。以下為一個實際心電圖信號采樣序列x(n)，其中存在高頻干擾，抽樣周期Ts=1秒。在實驗中，以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] 對序列x(n)用FFT做頻譜分析，生成x(n)的頻譜圖。

(2)設計一個巴特沃斯低通IIR數(shù)字濾波器H(z)。

設計指標參數(shù)為：在通帶內(nèi)頻率低于0.2π時，最大衰減小于1dB； 在阻帶內(nèi) [0.3π, π]頻率區(qū)間上，最小衰減大于15dB。

j?|H(e)|。寫出數(shù)字濾波器H(z)的表達式，畫出濾波器的幅頻響應曲線

(3)用所設計的濾波器對實際心電圖信號采樣序列x(n)進行濾波處理，編寫程序，求濾波后的序列y(n)，并分別畫出濾波前后的心電圖信號波形圖和頻譜圖。

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)濾波的心電圖信號 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('幅度響應dB');figure(2)plot(w,abs(h));grid;

xlabel('數(shù)字頻率');

ylabel('幅度響應|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)濾波的心電圖信號');xlabel('n');subplot(2,1,2)stem(l,xx);

title('經(jīng)濾波之后的心電圖信號');xlabel('n');figure(5)subplot(2,1,1)plot(l,abs(y));

title('未經(jīng)濾波的心電圖信號的頻譜');subplot(2,1,2)plot(l,abs(yy));

title('經(jīng)濾波處理的心電圖信號的頻譜');

2．用help查看內(nèi)部函數(shù)cheb1ord.m及cheby1.m，了解調(diào)用格式。

編程設計教材習題6-2，求模擬濾波器Ha(s)的表達式。

源程序

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('幅度響應')xlabel('頻率(Hz)')ylabel('幅度(dB)')grid

3．模擬濾波器的數(shù)字化

用內(nèi)部函數(shù)impinvar及bilinear實現(xiàn)教材習題6-5，求數(shù)字濾波器H(z)的表達式。

源程序

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)

本實驗所用的部分MATLAB函數(shù)

? L=length(x)：求序列x長度。

? y=fft(x,L)：將序列x(n)做L點快速傅立葉變換，結果賦給序列y(n)。? [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s')：計算模擬Butterworth濾波器的最小階次n和截止頻率為Wn。

? [b,a] = butter(n,Wn,'s')：設計模擬截止頻率為Wn（rad/s）的n階 Butterworth低通濾波器，返回值為模擬濾波器的系數(shù)。

? y=filter(b,a,x): 將序列x(n)通過濾波器濾波后生成序列y(n)，濾波器的分母多項式系數(shù)構成a向量，分子多項式系數(shù)構成b向量。

? [BZ,AZ] = impinvar(B,A,Fs)：沖激響應不變法，返回值為數(shù)字濾波器的系數(shù)。? [BZ,AZ] = bilinear(B,A,fs)：雙線性變換，返回值為數(shù)字濾波器的系數(shù)。? [H w]=freqz(b,a)：由濾波器分母多項式系數(shù)構成的a向量和分子多項式系數(shù)構成的b向量求系統(tǒng)頻響。

截圖

實驗體會

這次的實驗讓這讓我看到了理論與實踐相結合的優(yōu)勢與用處，讓我受益匪淺。我認識到了自己理論知識的不足，也認識到了我們學習的基礎知識究竟能運用于什么領域，如何運用。我們在老師的耐心指導下調(diào)試電路，直到得到要求的效果。讓我們在學習電路、信號等理論知識的同時，明白如何把這些應用于實際。