北京郵電大學(xué)DSP軟件實驗報告

北京郵電大學(xué)DSP軟件實驗報告學(xué)院：xxx學(xué)院姓名：xxx班級：xxxxxxxxxx學(xué)號：xxxxxxxxxx目錄實驗一：數(shù)字信號的FFT分析 21、實驗內(nèi)容及要求 22、實驗?zāi)康?23、設(shè)計思路 24、實驗代碼與實驗結(jié)果 3實驗二：DTMF信號的編碼 61、實驗內(nèi)容及要求 62、實驗?zāi)康?73、設(shè)計思路 74、實驗代碼與實驗結(jié)果 7實驗三：FIR數(shù)字濾波器的設(shè)計和實現(xiàn) 121、實驗內(nèi)容及要求： 122、實驗?zāi)康?133、設(shè)計思路 134、實驗代碼與運行結(jié)果 13實驗總結(jié) 17實驗一：數(shù)字信號的FFT分析1、實驗內(nèi)容及要求(1)離散信號的頻譜分析：設(shè)信號此信號的0.3pi和0.302pi兩根譜線相距很近，譜線0.45pi的幅度很小，請選擇合適的序列長度N和窗函數(shù)，用DFT分析其頻譜，要求得到清楚的三根譜線。(2)DTMF信號頻譜分析用計算機(jī)聲卡采用一段通信系統(tǒng)中電話雙音多頻（DTMF）撥號數(shù)字0～9的數(shù)據(jù)，采用快速傅立葉變換（FFT）分析這10個號碼DTMF撥號時的頻譜。2、實驗?zāi)康耐ㄟ^本次實驗，應(yīng)該掌握：(a)用傅立葉變換進(jìn)行信號分析時基本參數(shù)的選擇。(b)經(jīng)過離散時間傅立葉變換（DTFT）和有限長度離散傅立葉變換（DFT）后信號頻譜上的區(qū)別，前者DTFT時間域是離散信號，頻率域還是連續(xù)的，而DFT在兩個域中都是離散的。(c)離散傅立葉變換的基本原理、特性，以及經(jīng)典的快速算法（基2時間抽選法），體會快速算法的效率。(d)獲得一個高密度頻譜和高分辨率頻譜的概念和方法，建立頻率分辨率和時間分辨率的概念，為將來進(jìn)一步進(jìn)行時頻分析（例如小波）的學(xué)習(xí)和研究打下基礎(chǔ)。(e)建立DFT從整體上可看成是由窄帶相鄰濾波器組成的濾波器組的概念，此概念的一個典型應(yīng)用是數(shù)字音頻壓縮中的分析濾波器，例如DVDAC3和MPEGAudio。3、設(shè)計思路（1）由信號xn=0.001*cos0.45nπ+sin0.3nπ-cos?(0.302nπ-π4)函數(shù)限定坐標(biāo)軸范圍。（2）雙音多頻信號DTMF每個數(shù)字由兩個不同頻率的正弦波組成，低頻有：697Hz,770Hz,852Hz,941Hz，高頻有：1209Hz,1336Hz,1477Hz,1633Hz，0-9這十個數(shù)字每個數(shù)字對應(yīng)一個低頻信號和一個高頻信號疊加。分別用兩個數(shù)組裝載高頻和低頻，再產(chǎn)生由兩個正弦波疊加成的DTMF信號，最后利用plot和fft函數(shù)畫出對應(yīng)的頻譜圖。4、實驗代碼與實驗結(jié)果離散信號的頻譜分析【實驗代碼】clearcloseallN=20000;n=1:1:N;x=0.001*cos(0.45*pi*n)+sin(0.3*pi*n)-cos(0.302*pi*n-pi/4);y=fft(x,N);magy=abs(y(1:1:N/2+1));k=0:1:N/2;w=2*pi/N*k;stem(w/pi,magy);axis([0.25,0.5,0,60])【實驗結(jié)果】頻譜圖如下所示：圖SEQ圖\*ARABIC1DTMF信號頻譜分析【實驗代碼】clearcloseallcolumn=[1209,1336,1477,1633];line=[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);key(1,:)=cos(2*pi*column(1)*ts*n)+cos(2*pi*line(1)*ts*n);key(2,:)=cos(2*pi*column(2)*ts*n)+cos(2*pi*line(1)*ts*n);key(3,:)=cos(2*pi*column(3)*ts*n)+cos(2*pi*line(1)*ts*n);key(4,:)=cos(2*pi*column(1)*ts*n)+cos(2*pi*line(2)*ts*n);key(5,:)=cos(2*pi*column(2)*ts*n)+cos(2*pi*line(2)*ts*n);key(6,:)=cos(2*pi*column(3)*ts*n)+cos(2*pi*line(2)*ts*n);key(7,:)=cos(2*pi*column(1)*ts*n)+cos(2*pi*line(3)*ts*n);key(8,:)=cos(2*pi*column(2)*ts*n)+cos(2*pi*line(3)*ts*n);key(9,:)=cos(2*pi*column(3)*ts*n)+cos(2*pi*line(3)*ts*n);key(10,:)=cos(2*pi*column(2)*ts*n)+cos(2*pi*line(4)*ts*n);figure;fori=1:9subplot(3,4,i)plot(f,abs(fft(key(i,:))));xlabel('頻率（Hz）');ylabel('幅值');title(i);grid;endsubplot(3,4,10)plot(f,abs(fft(key(10,:))));xlabel('頻率（Hz）');ylabel('幅值');title(0);grid;【實驗結(jié)果】頻譜圖如下所示：圖SEQ圖\*ARABIC2實驗二：DTMF信號的編碼1、實驗內(nèi)容及要求（1）把您的聯(lián)系電話號碼通過DTMF編碼生成為一個.wav文件。技術(shù)指標(biāo)：根據(jù)ITUQ.23建議，DTMF信號的技術(shù)指標(biāo)是：傳送/接收率為每秒10個號碼，或每個號碼100ms。每個號碼傳送過程中，信號存在時間至少45ms，且不多于55ms，100ms的其余時間是靜音。在每個頻率點上允許有不超過±1.5%的頻率誤差。任何超過給定頻率±3.5%的信號，均被認(rèn)為是無效的，拒絕接收。（其中關(guān)鍵是不同頻率的正弦波的產(chǎn)生。可以使用查表方式模擬產(chǎn)生兩個不同頻率的正弦波。正弦表的制定要保證合成信號的頻率誤差在±1.5%以內(nèi)，同時使取樣點數(shù)盡量少）（2）對所生成的DTMF文件進(jìn)行解碼。DTMF信號解碼可以采用FFT計算N點頻率處的頻譜值，然后估計出所撥號碼。但FFT計算了許多不需要的值，計算量太大，而且為保證頻率分辨率，F(xiàn)FT的點數(shù)較大，不利于實時實現(xiàn)。因此，F(xiàn)FT不適合于DTMF信號解碼的應(yīng)用。由于只需要知道8個特定點的頻譜值，因此采用一種稱為Goertzel算法的IIR濾波器可以有效地提高計算效率。其傳遞函數(shù)為：2、實驗?zāi)康?a)復(fù)習(xí)和鞏固IIR數(shù)字濾波器的基本概念；(b)掌握IIR數(shù)字濾波器的設(shè)計方法；(c)掌握IIR數(shù)字濾波器的實現(xiàn)結(jié)構(gòu)；(d)能夠由濾波器的實現(xiàn)結(jié)構(gòu)分析濾波器的性能（字長效應(yīng)）；(e)了解通信系統(tǒng)電話DTMF撥號的基本原理和IIR濾波器實現(xiàn)方法。3、設(shè)計思路編碼：DTMF撥號鍵盤由一個4*4行列構(gòu)成，每列代表一個高頻信號，每行代表一個低頻信號，每當(dāng)按下一個鍵時，產(chǎn)生高、低頻率的兩個正弦信號，代表一個特定的數(shù)字或符號，根據(jù)ITUQ.23頒布的國際標(biāo)準(zhǔn)，DTMF傳送或接受每個號碼的時間為100ms，其中每個號碼傳送的過程中，信號存在時間至少45ms，其余時間靜音。用一個字符串變量來接受輸入的電話號碼，并將各個數(shù)字和符號的ASCII碼用一個4*4矩陣表示，每接收到一個數(shù)字就對應(yīng)兩個頻率，并產(chǎn)生由兩個正弦波疊加的信號，完成DTMF編碼，利用MATLAB提供的fft函數(shù)畫出其頻譜，用sound函數(shù)發(fā)出聲音。解碼：采用Goertzel算法來檢測DTMF信號，它是用IIR濾波器實現(xiàn)DFT算法的一種特殊方法，在實際DTMF解碼中，只需要知道輸入信號即DTMF信號的離散傅里葉變換X(k)的幅度信息，忽略相位信息，因為只要能得到8個特定頻率點的幅度值，看哪兩個頻率對應(yīng)的幅度最大，就能知道對應(yīng)的是哪個數(shù)字，達(dá)到解碼的目的，可以利用MATLAB提供的goertzel函數(shù)來對信號解碼。4、實驗代碼與實驗結(jié)果把聯(lián)系電話號碼通過DTMF編碼生成為一個.wav文件【實驗代碼】d=input('請輸入電話號碼：','s');sum=length(d);total_x=[];sum_x=[];sum_x=[sum_x,zeros(1,800)];fora=1:sumsymbol=abs(d(a));tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68];forp=1:4;forq=1:4;iftm(p,q)==abs(d(a));break,endendiftm(p,q)==abs(d(a));break,endendf1=[697,770,852,941];f2=[1209,1336,1477,1633];n=1:400;x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000);x=[x,zeros(1,400)];sum_x=sum_x+x;total_x=[total_x,x];endaudiowrite('soundwave.wav',total_x,8000);sound(total_x);t=(1:8800)/8000;subplot(2,1,1);plot(t,total_x);axis([0,1.2,-2,2]);xlabel('時間/s')title('DTMF信號時域波形')xk=fft(x);mxk=abs(xk);subplot(2,1,2);k=(1:800)*sum*8000/800;plot(k,mxk);xlabel('頻率');title('DTMF信號頻譜');disp('雙頻信號已生成并發(fā)出')【實驗結(jié)果】如下圖圖3所示，輸入本人電話號碼按下回車鍵之后能聽到撥碼聲音，并有圖4所示圖像生成。圖SEQ圖\*ARABIC3圖SEQ圖\*ARABIC4對所生成的DTMF文件進(jìn)行解碼【實驗代碼】k=[1820222431343842];N=205;disp(['接收端檢測到的號碼為'])fora=1:summ=800*(a-1);X=goertzel(total_x(m+1:m+N),k+1);val=abs(X);xk2=val.^2;limit=80;fors=5:8ifval(s)>limit,break,endendforr=1:4ifval(r)>limit,break,endenddisp([setstr(tm(r,s-4))])xk2end【實驗結(jié)果】在命令行窗口可以看到圖5-圖8的解碼信息，檢測到的號碼依次所得結(jié)果正確。圖SEQ圖\*ARABIC5圖SEQ圖\*ARABIC6圖SEQ圖\*ARABIC7圖SEQ圖\*ARABIC8實驗三：FIR數(shù)字濾波器的設(shè)計和實現(xiàn)1、實驗內(nèi)容及要求：錄制自己的一段聲音，長度為45秒，取樣頻率32kHz，然后疊加一個高斯白噪聲，使得信噪比為20dB。請采用窗口法設(shè)計一個FIR帶通濾波器，濾除噪聲提高質(zhì)量。提示：濾波器指標(biāo)參考：通帶邊緣頻率為4kHz，阻帶邊緣頻率為4.5kHz，阻帶衰減大于50dB；Matlab函數(shù)y=awgn(x,snr,'measured')，首先測量輸入信號x的功率，然后對其疊加高斯白噪聲；2、實驗?zāi)康耐ㄟ^本次實驗，掌握以下知識：FIR數(shù)字濾波器窗口設(shè)計法的原理和設(shè)計步驟；Gibbs效應(yīng)發(fā)生的原因和影響；不同類型的窗函數(shù)對濾波效果的影響，以及窗函數(shù)和長度N的選擇。3、設(shè)計思路首先通過audioread命令讀取指定地址的音頻文件，再利用MATLAB提供的awgn()函數(shù)為錄制的聲音信號加上高斯白噪聲，得到合成信號，通過audiowrite命令在指定地址生成新的音頻文件。要得到濾波器的沖激響應(yīng)或頻率響應(yīng)，首先根據(jù)濾波器的性能指標(biāo)得到窗函數(shù)的寬度N，再用布萊克曼窗得到濾波器的頻率響應(yīng)，經(jīng)過傅里葉反變換可得到濾波器的沖激相應(yīng)，把合成信號通過濾波器還原出原來的聲音信號。在figure(1)窗口中顯示濾波器的幅頻、相頻特性，在figure(2)窗口中顯示濾波前后的時域波形，在figure(3)窗口中顯示濾波前后的頻譜。4、實驗代碼與運行結(jié)果【實驗代碼】fs=32000;[x,fs]=audioread('mysound.wav');snr=20;x2=awgn(x,snr,'measured','db');audiowrite('sound_noise.wav',x2,fs);t=0:1/fs:(size(x2)-1)/fs;wp=8000*pi/32000;ws=9000*pi/32000;wdelta=ws-wp;N=ceil(11*pi/wdelta);wn=(ws+wp)/2;b=fir1(N,wn/pi,blackman(N+1));figure(1)freqz(b,1,512)f2=filter(b,1,x2);title('濾波器幅頻、相頻特性');figure(2)subplot(2,1,1)plot(t,x2)title('濾波前時域波形');subplot(2,1,2)plot(t,f2)title('濾波后時域波形');F0=fft(f2,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);s