北郵數字信號處理MATLAB實驗報告

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

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

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

4、，一個低頻，每個數字由兩個不同頻率的正弦波組成，高頻為：1209HZ,1336HZ,1633HZ,1477HZ,低頻為697HZ,770hz，852hz,941hz。分別用兩個數組來裝載高頻頻率和低頻頻率，然后再分別產生兩路頻率混疊的信號。再利用fft頻譜分析做出頻譜圖。4.實驗代碼及實驗結果1)離散信號的頻譜分析代碼clf;close all;%關閉所有圖形窗口N=1000;%DFT點數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);%對FFT結果求模w=2*pi/

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

6、) DTMF 信號頻譜分析代碼clearclose all%關閉所有窗口fh=1209,1336,1477,1633;%一個數組存放高頻fl=697,770,852,941;%一個數組存放低頻fs=8000;%采樣頻率N=1024;%采樣點數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);%數字鍵1key(2,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(1)*ts*n); %數字鍵2key

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

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

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

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

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

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

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

14、50,51,65;52,53,54,66;55,56,57,67;42,48,35,68; % DTMF表中鍵的16個ASCII碼n=1: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);%補零滿足占空比要求x02_z=x02,zeros(1,400);x03_z=x03,zeros(1,400);x04_z=x04,zeros(1,400);x05_z=x05,zeros(1,400); %補零滿足占空比要求x06_z=x06,zeros(1,400);x07_z=x07,zeros(1,40

17、0); %補零滿足占空比要求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);%畫圖函數sound(x_z);%發(fā)出聲音filename='phone.wav'%設置文件名audiow

18、rite(filename,x_z,fs);%寫入聲音文件k=18 20 22 24 31 34 38 42;%解碼N=205;%取依次諧波的DFT點數subplot(2,1,2);disp('解碼得到的號碼是：')%在命令行窗口顯示for i=1:11 %對11位的已經編好碼的電話號碼進行解碼 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個點處X(K)的模值，比較最大的兩個組合得到的數字即解碼結果。 set(gcf,'color','w');%設置圖片顯示背景為白色 shg;%顯示圖形窗口 pause; limit=80; for s=5:8; if v(s)>limit,break,end%查找列號endfor r=1:4; if v(r)>limit,break,end%查找行號enddisp(setstr(tm(r,s-4)%顯示解碼的字符end2）實驗結果 K值為18和31，對應數字1 K值為18和38，對應數字3 K值為20和34，對應數字5 K值為18

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

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

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

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

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

25、'); %疊加噪聲 s2=awgn(s1,20,'measured','db');%計算功率并疊加噪聲 s2k=fft(s2,N);%對疊加噪聲過后的信號進行fft變換 figure(2);%在圖形文件2中繪制疊加噪音后的時域頻域幅頻特性 subplot(2,1,1); plot(s2); title('時域（疊加噪聲）'); 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);%設計濾波器過程fp=4000;%通帶起始頻率fr=4500;%阻帶頻率wp=2*pi*fp/fs;%將濾波器指標轉換到數字域wr=2*pi*fr/fs;width=wr-wp;%帶寬N1=ceil(6.6*pi/width)+1;%取整數n=0:1:N1;wc=(wr+wp)/2;%計算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函數設計出滿足要求的濾波器w_hamming=(hamming(N1)'h=hd.*w_hamming;%加窗，即時域進行卷積H,w=freqz(h,1,1000,'whole');%求幅頻響應(fir濾波器的系統(tǒng)函數分子項為h，分母項為1)H=(H(1:1:501)'w=(w(1:1:501)'mag=abs(H);%求絕對值幅頻響應db=20*log10(mag+eps)/max(mag);%采用衰減特性，db做單位pha=angle(H);%相位響應delta_w=2*pi/1000;s3=filter(b,1,s2);%對s2進行濾波，s3為濾波后的信號filename='ashifil.wav'audiowrite(filename,s3,48000);%將濾波后的聲音文件寫入%進行濾波處理figure(3);在圖形文件3中繪制濾波后的波形文件的時域與頻域響應subplot(2,1,1);plot(s3);title('時域（濾波后）');s3k=fft(s3,N);%對濾波后的信號進行頻