北郵數(shù)字信號(hào)處理MATLAB實(shí)驗(yàn)報(bào)告

1、數(shù)字信號(hào)處理軟件實(shí)驗(yàn)MatLab仿真實(shí)驗(yàn)報(bào)告學(xué)院:電子工程學(xué)院班級(jí)：2013211202姓名：學(xué)號(hào)：實(shí)驗(yàn)一：數(shù)字信號(hào)的 FFT 分析1、實(shí)驗(yàn)內(nèi)容及要求(1) 離散信號(hào)的頻譜分析： 設(shè)信號(hào) 此信號(hào)的0.3pi 和 0.302pi兩根譜線相距很近，譜線 0.45pi 的幅度很小，請選擇合適的序列長度 N 和窗函數(shù)，用 DFT 分析其頻譜，要求得到清楚的三根譜線。(2) DTMF 信號(hào)頻譜分析 用計(jì)算機(jī)聲卡采用一段通信系統(tǒng)中電話雙音多頻（DTMF）撥號(hào)數(shù)字 09的數(shù)據(jù)，采用快速傅立葉變換（FFT）分析這10個(gè)號(hào)碼DTMF撥號(hào)時(shí)的頻譜。 2、實(shí)驗(yàn)?zāi)康?通過本次實(shí)驗(yàn)，應(yīng)該掌握：(a) 用傅立葉變換進(jìn)行信

2、號(hào)分析時(shí)基本參數(shù)的選擇。 (b) 經(jīng)過離散時(shí)間傅立葉變換（DTFT）和有限長度離散傅立葉變換（DFT） 后信號(hào)頻譜上的區(qū)別，前者 DTFT 時(shí)間域是離散信號(hào)，頻率域還是連續(xù)的，而 DFT 在兩個(gè)域中都是離散的。(c) 離散傅立葉變換的基本原理、特性，以及經(jīng)典的快速算法（基2時(shí)間抽選法），體會(huì)快速算法的效率。(d) 獲得一個(gè)高密度頻譜和高分辨率頻譜的概念和方法，建立頻率分辨率和時(shí)間分辨率的概念，為將來進(jìn)一步進(jìn)行時(shí)頻分析（例如小波）的學(xué)習(xí)和研究打下基礎(chǔ)。(e) 建立 DFT 從整體上可看成是由窄帶相鄰濾波器組成的濾波器組的概念，此概念的一個(gè)典型應(yīng)用是數(shù)字音頻壓縮中的分析濾波器，例如 DVD AC3

3、 和MPEG Audio。3.設(shè)計(jì)思路及實(shí)驗(yàn)步驟1）離散信號(hào)的頻譜分析：該信號(hào)中要求能夠清楚的觀察到三根譜線。由于頻率0.3pi和0.302pi間隔非常小，要清楚的顯示，必須采取足夠大小的N，使得分辨率足夠好，至少到0.001單位級(jí)，而頻率0.45pi的幅度很小，要清楚的觀察到它的譜線，必須采取幅度夠大的窗函數(shù)，使得它的頻譜幅度變大一些。同時(shí)還要注意頻譜泄漏的問題，三個(gè)正弦函數(shù)的周期（2pi/w）分別為20,40,1000，所以為了避免產(chǎn)生頻譜泄漏(k=w/w0為整數(shù))，采樣點(diǎn)數(shù)N必須為1000的整數(shù)倍。2）DTMF 信號(hào)頻譜分析雙音多頻信號(hào)中4*4撥號(hào)盤中的每一個(gè)按鍵均對(duì)應(yīng)兩個(gè)頻率，一個(gè)高頻

4、，一個(gè)低頻，每個(gè)數(shù)字由兩個(gè)不同頻率的正弦波組成，高頻為：1209HZ,1336HZ,1633HZ,1477HZ,低頻為697HZ,770hz，852hz,941hz。分別用兩個(gè)數(shù)組來裝載高頻頻率和低頻頻率，然后再分別產(chǎn)生兩路頻率混疊的信號(hào)。再利用fft頻譜分析做出頻譜圖。4.實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果1)離散信號(hào)的頻譜分析代碼clf;close all;%關(guān)閉所有圖形窗口N=1000;%DFT點(diǎn)數(shù)n=1:1:N;x=0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-pi/4);y=fft(x,N);mag=abs(y);%對(duì)FFT結(jié)果求模w=2*pi/

5、N*0:1:N-1; %數(shù)字角頻率wsubplot(2,1,1); %將圖形窗分為2行1列stem(n,x,'.'); %畫脈沖圖title('時(shí)域波形');xlabel('n');ylabel('x(n)');subplot(2,1,2);stem(w/pi,mag); %歸一化角頻率axis(0.2 0.5 0 2); %控制坐標(biāo)范圍以使譜線幅度合適title('1000點(diǎn)DFT頻譜分析');xlabel('數(shù)字頻率');ylabel('X(k)');grid on;實(shí)驗(yàn)結(jié)果：2

6、) DTMF 信號(hào)頻譜分析代碼clearclose all%關(guān)閉所有窗口fh=1209,1336,1477,1633;%一個(gè)數(shù)組存放高頻fl=697,770,852,941;%一個(gè)數(shù)組存放低頻fs=8000;%采樣頻率N=1024;%采樣點(diǎn)數(shù)ts=1/fs;%采樣周期n=0:N-1;f=0:fs/N:fs/N*(N-1);key=zeros(10,N);%生成0矩陣key(1,:)=cos(2*pi*fh(1)*ts*n)+cos(2*pi*fl(1)*ts*n);%數(shù)字鍵1key(2,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(1)*ts*n); %數(shù)字鍵2key

7、(3,:)=cos(2*pi*fh(3)*ts*n)+cos(2*pi*fl(1)*ts*n); %數(shù)字鍵3key(4,:)=cos(2*pi*fh(1)*ts*n)+cos(2*pi*fl(2)*ts*n); %數(shù)字鍵4key(5,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(2)*ts*n); %數(shù)字鍵5key(6,:)=cos(2*pi*fh(3)*ts*n)+cos(2*pi*fl(2)*ts*n); %數(shù)字鍵6key(7,:)=cos(2*pi*fh(1)*ts*n)+cos(2*pi*fl(3)*ts*n); %數(shù)字鍵7key(8,:)=cos(2*pi*

8、fh(2)*ts*n)+cos(2*pi*fl(3)*ts*n); %數(shù)字鍵8key(9,:)=cos(2*pi*fh(3)*ts*n)+cos(2*pi*fl(3)*ts*n); %數(shù)字鍵9key(10,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(4)*ts*n); %數(shù)字鍵0figure;for i=1:9%畫出數(shù)字鍵1到9的頻譜subplot(3,4,i)plot(f,abs(fft(key(i,:);%利用fft求各個(gè)鍵的頻譜axis(500 1700 0 1000);%限定坐標(biāo)軸范圍xlabel('頻率（hz）');ylabel('

9、幅度');%橫縱坐標(biāo)標(biāo)識(shí)title(i);grid;subplot(3,4,10)%畫出數(shù)字鍵0的頻譜圖plot(f,abs(fft(key(10,:);axis(500 1700 0 1000);xlabel('頻率（hz）');ylabel(幅度');title(0);grid;endDTMF 信號(hào)頻譜分析結(jié)果實(shí)驗(yàn)二： DTMF 信號(hào)的編碼1、實(shí)驗(yàn)內(nèi)容及要求1）把您的聯(lián)系電話號(hào)碼 通過DTMF 編碼生成為一個(gè) .wav 文件。³ 技術(shù)指標(biāo)：± 根據(jù) ITU Q.23 建議，DTMF 信號(hào)的技術(shù)指標(biāo)是：傳送/接收率為每秒 10 個(gè)號(hào)碼，或每

10、個(gè)號(hào)碼 100ms。± 每個(gè)號(hào)碼傳送過程中，信號(hào)存在時(shí)間至少 45ms，且不多于 55ms，100ms 的其余時(shí)間是靜音。± 在每個(gè)頻率點(diǎn)上允許有不超過 ±1.5% 的頻率誤差。任何超過給定頻率 ±3.5% 的信號(hào)，均被認(rèn)為是無效的，拒絕接收。（其中關(guān)鍵是不同頻率的正弦波的產(chǎn)生。可以使用查表方式模擬產(chǎn)生兩個(gè)不同頻率的正弦波。正弦表的制定要保證合成信號(hào)的頻率誤差在±1.5%以內(nèi)，同時(shí)使取樣點(diǎn)數(shù)盡量少） 2）對(duì)所生成的DTMF文件進(jìn)行解碼。³ DTMF 信號(hào)解碼可以采用 FFT 計(jì)算 N 點(diǎn)頻率處的頻譜值，然后估計(jì)出所撥號(hào)碼。但 FFT計(jì)

11、算了許多不需要的值，計(jì)算量太大，而且為保證頻率分辨率，F(xiàn)FT的點(diǎn)數(shù)較大，不利于實(shí)時(shí)實(shí)現(xiàn)。因此，F(xiàn)FT 不適合于 DTMF 信號(hào)解碼的應(yīng)用。³ 由于只需要知道 8 個(gè)特定點(diǎn)的頻譜值，因此采用一種稱為 Goertzel 算法的 IIR 濾波器可以有效地提高計(jì)算效率。其傳遞函數(shù)為：³2、實(shí)驗(yàn)?zāi)康?a)復(fù)習(xí)和鞏固 IIR 數(shù)字濾波器的基本概念；(b)掌握 IIR 數(shù)字濾波器的設(shè)計(jì)方法；(c)掌握 IIR 數(shù)字濾波器的實(shí)現(xiàn)結(jié)構(gòu)；(d)能夠由濾波器的實(shí)現(xiàn)結(jié)構(gòu)分析濾波器的性能（字長效應(yīng)）；(e)了解通信系統(tǒng)電話 DTMF 撥號(hào)的基本原理和 IIR 濾波器實(shí)現(xiàn)方法。3.設(shè)計(jì)思路及實(shí)驗(yàn)步驟1

12、）編碼：DTMF 信號(hào)是將撥號(hào)盤上的 0F 共16 個(gè)數(shù)字，用音頻范圍的 8 個(gè)頻率來表示的一種編碼方式。8 個(gè)頻率分為高頻群和低頻群兩組，分別作為列頻和行頻。 每個(gè)字符的信號(hào)由來自列頻和行頻的兩個(gè)頻率的正弦信號(hào)疊加而成，組合方式如表格所示1209hz1336hz1477hz16333hz697hz123A770hz456B852hz789C941hz*0#D根據(jù)表格即可得到各個(gè)數(shù)字對(duì)應(yīng)的DTFM信號(hào)。通過zeros全零矩陣來設(shè)置占空比，以達(dá)到題目要求。得到信號(hào)后，使用sound函數(shù)來播放撥號(hào)音，writewave

13、將信號(hào)寫入聲音文件。2）解碼：在解碼時(shí)，使用Goertzel算法。濾波器調(diào)諧到這8個(gè)頻率之上后，在相應(yīng)的頻率上的頻譜值最大，通過與標(biāo)準(zhǔn)值的對(duì)比找出在DTMF圖中的行和列，再對(duì)應(yīng)出相應(yīng)的撥號(hào)數(shù)字。由Mock確定的Goertzel算法參數(shù)，假設(shè)一次諧波的DFT長度為205，此時(shí)對(duì)應(yīng)的離散頻率點(diǎn)k的值為18、20、22、24、31、34、38、42。依次對(duì)應(yīng)頻率697hz,770hz,852hz,941hz,1209hz,1336hz,1447hz,1633hz。4.實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果1）實(shí)驗(yàn)代碼：（我的電話號(hào)碼N=400;fs=8000; %每個(gè)號(hào)碼100mstm=49,

14、50,51,65;52,53,54,66;55,56,57,67;42,48,35,68; % DTMF表中鍵的16個(gè)ASCII碼n=1:N;%取樣點(diǎn)fl=697 770 852 941;%低頻fh=1209 1336 1477 1633;%高頻x01=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(1)*n/fs); %1編碼過程x02=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(3)*n/fs); %3x03=sin(2*pi*fl(2)*n/fs)+sin(2*pi*fh(2)*n/fs); %5x04=sin(2*pi*fl(1)*n/

15、fs)+sin(2*pi*fh(2)*n/fs); %2x05=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(1)*n/fs); %1x06=sin(2*pi*fl(3)*n/fs)+sin(2*pi*fh(3)*n/fs); %9x07=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(1)*n/fs);%1x08=sin(2*pi*fl(2)*n/fs)+sin(2*pi*fh(3)*n/fs); %6x09=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(3)*n/fs); %3x10=sin(2*pi*fl(4)*n/fs)+sin

16、(2*pi*fh(2)*n/fs); %0x11=sin(2*pi*fl(3)*n/fs)+sin(2*pi*fh(1)*n/fs); %7x=x01,x02,x03,x04,x05,x06,x07,x08,x09,x10,x11;%組成矩陣x01_z=x01,zeros(1,400);%補(bǔ)零滿足占空比要求x02_z=x02,zeros(1,400);x03_z=x03,zeros(1,400);x04_z=x04,zeros(1,400);x05_z=x05,zeros(1,400); %補(bǔ)零滿足占空比要求x06_z=x06,zeros(1,400);x07_z=x07,zeros(1,40

17、0); %補(bǔ)零滿足占空比要求x08_z=x08,zeros(1,400);x09_z=x09,zeros(1,400);x10_z=x10,zeros(1,400);x11_z=x11,zeros(1,400);x_z=x01_z,x02_z,x03_z,x04_z,x05_z,x06_z,x07_z,x08_z,x09_z,x10_z,x11_z;%組成矩陣x_z=x_z/max(abs(x_z);subplot(2,1,1);%將圖示屏幕分為兩行一列plot(x_z);%畫圖函數(shù)sound(x_z);%發(fā)出聲音filename='phone.wav'%設(shè)置文件名audiow

18、rite(filename,x_z,fs);%寫入聲音文件k=18 20 22 24 31 34 38 42;%解碼N=205;%取依次諧波的DFT點(diǎn)數(shù)subplot(2,1,2);disp('解碼得到的號(hào)碼是：')%在命令行窗口顯示for i=1:11 %對(duì)11位的已經(jīng)編好碼的電話號(hào)碼進(jìn)行解碼 m=400*(i-1); X=goertzel(x(m+1:m+N),k+1); %goertzel算法做變換 v=abs(X); %求模 stem(k,v,'.'); %畫脈沖圖 grid; xlabel('k'); ylabel('x(k)&

19、#39;);%作出8個(gè)點(diǎn)處X(K)的模值，比較最大的兩個(gè)組合得到的數(shù)字即解碼結(jié)果。 set(gcf,'color','w');%設(shè)置圖片顯示背景為白色 shg;%顯示圖形窗口 pause; limit=80; for s=5:8; if v(s)>limit,break,end%查找列號(hào)endfor r=1:4; if v(r)>limit,break,end%查找行號(hào)enddisp(setstr(tm(r,s-4)%顯示解碼的字符end2）實(shí)驗(yàn)結(jié)果 K值為18和31，對(duì)應(yīng)數(shù)字1 K值為18和38，對(duì)應(yīng)數(shù)字3 K值為20和34，對(duì)應(yīng)數(shù)字5 K值為18

20、和34，對(duì)應(yīng)數(shù)字2 K值為18和31，對(duì)應(yīng)數(shù)字1 K值為22和38，對(duì)應(yīng)數(shù)字9 K值為18和31，對(duì)應(yīng)數(shù)字1 K值為20和38，對(duì)應(yīng)數(shù)字6 K值為18和38，對(duì)應(yīng)數(shù)字3 K值為24和34，對(duì)應(yīng)數(shù)字0 K值為22和31，對(duì)應(yīng)數(shù)字7在命令窗口顯示對(duì)應(yīng)號(hào)碼實(shí)驗(yàn)三：FIR 數(shù)字濾波器的設(shè)計(jì)和實(shí)現(xiàn)1、實(shí)驗(yàn)內(nèi)容及要求：³ 錄制自己的一段聲音，長度為 10秒，取樣頻率 32kHz，然后疊加一個(gè)高斯白噪聲，使得信噪比為 20dB。請采用窗口法設(shè)計(jì)一個(gè) FIR 帶通濾波器，濾除噪聲提高質(zhì)量。n 提示：³ 濾波器指標(biāo)參考：通帶邊緣頻率為 4kHz，阻帶邊緣頻率為4.5kHz，阻帶衰減大于 50

21、dB；³ Matlab 函數(shù) y = awgn(x,snr,'measured') ，首先測量輸入信號(hào) x 的功率，然后對(duì)其疊加高斯白噪聲；2、實(shí)驗(yàn)?zāi)康?#179; 通過本次實(shí)驗(yàn)，掌握以下知識(shí)：± FIR 數(shù)字濾波器窗口設(shè)計(jì)法的原理和設(shè)計(jì)步驟；± Gibbs 效應(yīng)發(fā)生的原因和影響；± 不同類型的窗函數(shù)對(duì)濾波效果的影響，以及窗函數(shù)和長度 N 的選擇。3.設(shè)計(jì)思路及實(shí)驗(yàn)步驟1）fir濾波器設(shè)計(jì)原理：設(shè)計(jì)fir濾波器的方法有很多，窗函數(shù)設(shè)計(jì)法，頻率采樣法，最優(yōu)化設(shè)計(jì)法等，該實(shí)驗(yàn)我們采用窗函數(shù)設(shè)計(jì)法，其基本原理是用一定寬度的窗函數(shù)截取無限脈沖響應(yīng)

22、序列獲得有限長的脈沖響應(yīng)序列。其設(shè)計(jì)步驟主要有·通過傅里葉逆變換獲得理想濾波器的單位脈沖響應(yīng)。·由性能指標(biāo)確定窗函數(shù)的種類和窗口長度·求得實(shí)際濾波器的單位脈沖響應(yīng)h(n),即為所設(shè)計(jì)的FIR濾波器系數(shù)向量·檢驗(yàn)濾波器性能。2）窗函數(shù)的選?。簽V波器的指標(biāo)要求為 通帶邊緣為4KHZ，阻帶邊緣頻率為4.5Khz，阻帶衰減大于50db。接下來應(yīng)選取合適的窗函數(shù)實(shí)現(xiàn)濾波指標(biāo)，然后計(jì)算出單位脈沖響應(yīng)，再利用函數(shù)freqz()求出其幅頻特性看是否滿足濾波器要求。由于直接對(duì)信號(hào)加矩形窗截?cái)鄷?huì)產(chǎn)生頻譜泄露，為了改善頻率泄露的問題，通常加非矩形窗，一般加漢明窗，因

23、為漢明窗的幅頻特性是旁瓣衰減較大，主瓣峰值與第一個(gè)旁瓣衰減可達(dá)50db。在實(shí)驗(yàn)中使用Fir1()設(shè)計(jì)滿足濾波器指標(biāo)要求的濾波器。3）其次需要實(shí)現(xiàn)在聲音上面疊加噪聲，使用awgn函數(shù)。使用作圖函數(shù)分別畫出原始聲音信號(hào)和加噪聲后的聲音信號(hào)的時(shí)域和頻域波形。4）濾波操作：使用matlab庫函數(shù)hamming來進(jìn)行濾波。再作出濾波后的圖形。最后作出濾波器的幅頻特性圖。4.實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果1）實(shí)驗(yàn)代碼close all;%關(guān)閉所有圖形窗口 %原始聲音信號(hào)頻譜 s1,fs=audioread('ashi.wav');%在當(dāng)前文件夾下讀取音頻文件 fs=32000; %采樣頻率32khz

24、N=4096; %fft點(diǎn)數(shù) k=1:4096; %一個(gè)K 對(duì)應(yīng)一個(gè)數(shù)字頻率 s1k=fft(s1,N); %對(duì)原始音頻信號(hào)進(jìn)行頻譜分析figure(1);%在圖形文件1中顯示原始音頻信號(hào)的時(shí)域和頻譜 subplot(2,1,1); plot(s1); title('原始信號(hào)（時(shí)域');subplot(2,1,2); plot(32/4096*k,abs(s1k); %頻率分辨率=fs/N axis(0 6 0 7);%設(shè)定橫縱坐標(biāo)范圍 xlabel('f(kHz)'); ylabel('s1(k)'); title('原聲音信號(hào)（頻域）

25、'); %疊加噪聲 s2=awgn(s1,20,'measured','db');%計(jì)算功率并疊加噪聲 s2k=fft(s2,N);%對(duì)疊加噪聲過后的信號(hào)進(jìn)行fft變換 figure(2);%在圖形文件2中繪制疊加噪音后的時(shí)域頻域幅頻特性 subplot(2,1,1); plot(s2); title('時(shí)域（疊加噪聲）'); filename='noise.wav'audiowrite(filename,s2,48000);%將噪聲文件寫入 subplot(2,1,2); plot(32/4096*k,abs(s2k);

26、 %頻率分辨率=fs/Ntitle('頻域（疊加噪聲后）');xlabel('f(kHz)');ylabel('s2(k)');axis(0 6 0 7);%設(shè)計(jì)濾波器過程fp=4000;%通帶起始頻率fr=4500;%阻帶頻率wp=2*pi*fp/fs;%將濾波器指標(biāo)轉(zhuǎn)換到數(shù)字域wr=2*pi*fr/fs;width=wr-wp;%帶寬N1=ceil(6.6*pi/width)+1;%取整數(shù)n=0:1:N1;wc=(wr+wp)/2;%計(jì)算wcalpha=(N1-1)/2;n=0:1:N1-1;m=n-alpha+eps;hd=sin(wc*m

27、)./(pi*m);b=fir1(N1,wc/(4*pi/3.65);%返回截止頻率為wc的N1階濾波器hn，漢明窗%通過Fir函數(shù)設(shè)計(jì)出滿足要求的濾波器w_hamming=(hamming(N1)'h=hd.*w_hamming;%加窗，即時(shí)域進(jìn)行卷積H,w=freqz(h,1,1000,'whole');%求幅頻響應(yīng)(fir濾波器的系統(tǒng)函數(shù)分子項(xiàng)為h，分母項(xiàng)為1)H=(H(1:1:501)'w=(w(1:1:501)'mag=abs(H);%求絕對(duì)值幅頻響應(yīng)db=20*log10(mag+eps)/max(mag);%采用衰減特性，db做單位pha=angle(H);%相位響應(yīng)delta_w=2*pi/1000;s3=filter(b,1,s2);%對(duì)s2進(jìn)行濾波，s3為濾波后的信號(hào)filename='ashifil.wav'audiowrite(filename,s3,48000);%將濾波后的聲音文件寫入%進(jìn)行濾波處理figure(3);在圖形文件3中繪制濾波后的波形文件的時(shí)域與頻域響應(yīng)subplot(2,1,1);plot(s3);title('時(shí)域（濾波后）');s3k=fft(s3,N);%對(duì)濾波后的信號(hào)進(jìn)行頻