雙音多頻(DTMF)通信設(shè)計(jì)的MATLAB仿真摘要：討論以MATLAB作為仿真工具產(chǎn)生DTMF信號(hào)，并用FFT算法、DFT算法、卷積法及迭代法來對(duì)DTMF信號(hào)進(jìn)行解碼。關(guān)鍵詞：FFT；DFT；頻譜分析；卷積；濾波；差分方程；MATLAB0引言雙音多頻(DTMF:DoubleToneMulti-Frequency)是按鍵電話通信，也廣泛應(yīng)用于電子郵件和銀行系統(tǒng)中。用戶可從電話發(fā)送DTMF信號(hào)來選擇菜單進(jìn)行操作。DTMF信號(hào)容易用軟件產(chǎn)生和解碼。MATLAB是一個(gè)高度集成的軟件系統(tǒng)，通過交互式的命令(語句)可以十分簡(jiǎn)便地實(shí)現(xiàn)許多復(fù)雜的數(shù)值計(jì)算。本文采用MATLAB作為仿真工具產(chǎn)生DTMF信號(hào)，并用FFT算法、DFT算法、卷積法及迭代法來對(duì)DTMF信號(hào)進(jìn)行解碼，由此得出：時(shí)域和頻域是研究信號(hào)的兩個(gè)窗口，其中信號(hào)處理大都在時(shí)域中進(jìn)行，而信號(hào)分析往往在頻域中比較方便直觀。且數(shù)字信號(hào)處理技術(shù)中的DFT、FFT、卷積、濾波、差分方程這幾個(gè)概念之間有內(nèi)在聯(lián)系。1DTMF信號(hào)的產(chǎn)生DTMF是數(shù)字音頻信號(hào)，在DTMF通信系統(tǒng)中共有8個(gè)頻率，分為4個(gè)高頻音和4個(gè)低頻音，用一個(gè)高頻音和一個(gè)低頻音的組合表示一個(gè)信號(hào)，這樣共有16種組合，分別代表16種信號(hào)，如表1所示：7^)f(Hz)\1209133614771633L697123A770456B852789C941*0#D表1DTMF信號(hào)組合表例如，當(dāng)按下數(shù)字鍵“1”時(shí)，則產(chǎn)生低頻697Hz和高頻1209Hz這兩個(gè)正弦信號(hào)的迭加。由于語音信號(hào)的最高頻率為4KHz，根據(jù)奈奎斯特取樣定理，取樣頻率fs應(yīng)大于或等于原信號(hào)最高頻率fc的兩倍，即fsN2fc (1)才能保證取樣后的信號(hào)不失真，所以電話音頻信號(hào)在數(shù)字信號(hào)處理時(shí)，取樣頻率fs為2X4k=8kHz,這里，每個(gè)數(shù)字信號(hào)持續(xù)時(shí)間為100ms,后面加上100ms的間隔時(shí)間(用0表示)。上述DTMF信號(hào)產(chǎn)生方法如下：建立撥號(hào)數(shù)字的表矩陣，用查表法(查表1)求用戶所按數(shù)字鍵對(duì)應(yīng)的高、低頻音。為簡(jiǎn)化起見，僅允許選擇“0-9”這十個(gè)鍵，在開始時(shí)還可撥空信號(hào)。產(chǎn)生相應(yīng)的DTMF信號(hào)及間隔時(shí)間。由于fs=8kHz,各信號(hào)持續(xù)時(shí)間為100ms,因此在程序中每個(gè)信號(hào)取800點(diǎn)，間隔時(shí)間也取800點(diǎn)，結(jié)果存入數(shù)組中。畫圖并監(jiān)聽產(chǎn)生的DTMF信號(hào)。程序如下：%[程序1],DTMF信號(hào)的產(chǎn)生clear 清%除內(nèi)存TAB=[9411336;6971209;6971336;6971477;7701209;7701336;7701477;8521209;8521336;8521477];%撥號(hào)數(shù)字表矩陣n=input('n=');%DTMF信號(hào)的個(gè)數(shù)l=input('n0=');%空信號(hào)點(diǎn)數(shù)fori=1:nk=input('0?9');%輸入的數(shù)字鍵fL=TAB(k+1,1);%k對(duì)應(yīng)的低頻音fH=TAB(k+1,2);%k對(duì)應(yīng)的高頻音n1=800;fs=8000;%產(chǎn)生k對(duì)應(yīng)的DTMF信號(hào)，取樣頻率8kHz,每個(gè)信號(hào)共1600點(diǎn)，其中，j=0:1:n1-1; 前800點(diǎn)為信號(hào)持續(xù)時(shí)間，后800點(diǎn)為間隔時(shí)間，結(jié)果存入數(shù)組out中。x=sin(2*pi*fL*j/fs)+sin(2*pi*fH*j/fs);out(1600*(i-1)+1+l:1600*i-800+l)=x;out(1600*i-799+l:1600*i+l)=0;endout=out./2; %數(shù)組中每個(gè)數(shù)據(jù)除以2subplot(211);plot(out);%繪圖并監(jiān)聽DTMF信號(hào)。sound(out,fs)wavwrite(out,fs,'D2.wav');%out另存入聲音文件程序結(jié)果：例如：若要產(chǎn)生數(shù)字鍵“1”所對(duì)應(yīng)的DTMF信號(hào)，且開始時(shí)有200點(diǎn)空信號(hào)，則程序運(yùn)行如下：詢問：輸入：n=1n0= 2000-9 1得圖1：圖1DTMF信號(hào)圖圖1中，前200點(diǎn)為空信號(hào)，第201-1000點(diǎn)為鍵“1”對(duì)應(yīng)的DTMF信號(hào)，第1001-1800點(diǎn)為間隔時(shí)間0.5。-0.5-10 200 400 600 800 1000 12001400 1600 1800另外，在程序中為使軟件設(shè)計(jì)更接近于實(shí)際硬件的開發(fā)應(yīng)用，可用求解差分方程的方法來代替正弦函數(shù)的調(diào)用。設(shè)正弦序列為h(n)=sin(3kn)u(n),為實(shí)時(shí)實(shí)現(xiàn)h(n)，必須找到其滿足的差分方程。正弦序列h(n)的z變換為：H(z)= zsi「3k ,令H(z)分子/分母系數(shù)：sin3=b,2cos3=aZ2一2cos3z+1 k k則H(z)=Y(z)=—zb——=——^^——,Y(z)(1-az-1+z2)=bz-iX(z),X(z)z2-az+1 1-az-i+z-2兩邊進(jìn)行反變換，得y(n)-ay(n-1)+y(n-2)=bx(n-1)，式中，若令x(n)=6(n)，則得到h(n)的差分方程為h(n)-ah(n-1)+h(n-2)=b6(n-1),即：h(n)=ah(n-1)-h(n-2)+bd(n-1) (2)用迭代法解此差分方程，即得數(shù)字頻率為3k的正弦序列h(n)。在硬件中，該差分方程是由加法器、乘法器和單位延時(shí)單元構(gòu)成的系統(tǒng)。將x(n)=6(n)輸入該系統(tǒng)后，輸出的就是h(n)。本文中，每個(gè)DTMF信號(hào)h(n)是兩個(gè)頻率的正弦序列相迭加，設(shè)為hL(n)和h「n),為此，分別求得hL(n)和hH(n)所滿足的差分方程： LHh(n)=ah(n-1)-h(n-2)+b6(n-1)；h(n)=ah(n-1)-h(n-2)+b6(n-1)L LL L L H HH H H則h(n)=hL(n)+hH(n)。相關(guān)程序如下：%[程序2]，用解差分方程方法產(chǎn)生DTMF信號(hào)。x=zeros(1,800);x(2)=1;%x(n)=6(n-1)(x(n)取800點(diǎn))n=input('n=')； %DTM信號(hào)的個(gè)數(shù)l=input('n0=')； 空信號(hào)的點(diǎn)數(shù)out=zeros(1,1600*n+l);%out數(shù)組初始化w=2*pi/8000*[9411336;6971209;6971336;6971477;7701209;7701336;7701477;8521209;8521336;8521477];%各信號(hào)對(duì)應(yīng)的數(shù)字頻率tab=[2*cos(w)sin(w)];%各信號(hào)滿足的差分方程h〔(n)和hH(n)中系數(shù)a,b的矩陣。fori=1:nk=input'0-9-'); 輸入的數(shù)字鍵hL=zeros(1,3);hH=zeros(1,3);%]h(n)和hH(n)初始化forj=1:800 迭代法求解hL(n)和hH(n)產(chǎn)生相應(yīng)的DTMF信號(hào)，存AouthL(3)=tab(k+1,1)*hL(2)-hL(1)+tab(k+1,3)*x(j);hL(1)=hL(2);hL(2)=hL(3);hH(3)=tab(k+1,2)*hH(2)-hH(1)+tab(k+1,4)*x(j);hH(1)=hH(2);hH(2)=hH(3);out(1600*(i-1)+j+l)=hL(3)+hH(3);out(1600*(i-1)+j+800+l)=0;endend2DTMF信號(hào)的解碼解碼就是對(duì)接收到的DTMF信號(hào)進(jìn)行頻譜分析，從中找出代表各信號(hào)的的特征字，由此獲知用戶按下的數(shù)字鍵。為在頻譜圖中分辨出不同的頻率分量，于是對(duì)信號(hào)取200點(diǎn)為一幀，則頻譜分辨率F=f/N=8000/200=40Hz<73Hz（表1中任意兩頻率的最小間隔），這樣，即可滿足s頻譜分析的要求。本文采用數(shù)字信號(hào)處理技術(shù)中的FFT算法、DFT算法、卷積法及迭代法這4種算法實(shí)現(xiàn)對(duì)DTMF信號(hào)的解碼。一、 快速傅里葉變換（FFT）算法FFT是有限長(zhǎng)序列離散傅里葉變換（DFT）的快速算法，其基本運(yùn)算是蝶形算法，它使DFT計(jì)算時(shí)間縮短了幾個(gè)數(shù)量級(jí)，在信號(hào)處理中占有極重要的地位。這里采用基2-FFT算法對(duì)DTMF信號(hào)進(jìn)行頻譜分析。解碼過程如下：（1） 接收DTMF信號(hào)，并畫圖。（2） 對(duì)信號(hào)作FFT，畫頻譜圖，從中找出代表各信號(hào)的頻率分量。這部分，信號(hào)將完成從時(shí)域到頻域的轉(zhuǎn)換。1） 每幀信號(hào)（200點(diǎn)）做一次N=256點(diǎn)的FFT，從中取64點(diǎn)畫頻譜圖。在MATLAB中,FFT可由語句“y=fft（x,N）”來實(shí)現(xiàn)。而FFT中，要求序列長(zhǎng)度N=2e（E為整數(shù)），所以N=28=256，頻譜分辨率F=fs/N^31.25Hzo因?yàn)樾盘?hào)x為實(shí)數(shù)序列，所其幅頻譜|y|具有偶對(duì)稱性，于是，幅頻譜可以僅畫N/2點(diǎn)，其中第N/2點(diǎn)對(duì)應(yīng)實(shí)際頻率為fs/2=4KHz,而DTMF信號(hào)中最高頻率為1633Hz，小于2KHz（fs/4），因此，這里只畫N/4=64點(diǎn)。DTMF>號(hào)是兩個(gè)正弦波的迭加，它的幅頻譜就是兩根譜線，譜線的橫坐標(biāo)就是該信號(hào)的兩個(gè)頻率分量點(diǎn)kl和kh。2） ；肖除頻譜泄漏現(xiàn)象。由于信號(hào)x是有限長(zhǎng)的，這就相當(dāng)于對(duì)無限長(zhǎng)的信號(hào)加矩形窗，所以在頻譜圖中必然會(huì)出現(xiàn)頻譜泄漏現(xiàn)象，使信號(hào)能量散布到其他譜線位置。為此，應(yīng)選擇一適當(dāng)閥值，將出現(xiàn)在這兩條譜線周圍的幅度較小的譜線消除（置）），從而解決了這一問題。最后，將處理后的幅頻譜數(shù)據(jù)存入數(shù)組c中。3） 將各DTMF信號(hào)還原為相應(yīng)的數(shù)字鍵。在幅頻譜圖中，頻率軸的定標(biāo)方式為頻率點(diǎn)K而不是實(shí)際頻率f,轉(zhuǎn)換關(guān)系為：K=f/F,因此，數(shù)字鍵0-9對(duì)應(yīng)頻率點(diǎn)如下表所示："...404448231232645628789310表2數(shù)字鍵對(duì)應(yīng)頻率點(diǎn)組合表數(shù)組c的不等于0的下標(biāo)就是各信號(hào)的頻率點(diǎn)，查表2,即可將各DTMF信號(hào)還原為相應(yīng)的數(shù)字鍵。相關(guān)程序如下：%[程序3],FFT算法解碼程序。A=wavread('D2.wav');%接收到的DTMF信號(hào)subplot(212);plot(A);%繪圖N=256;fors=1:8*n 對(duì)%每幀信號(hào)作N=256點(diǎn)的FFTR=out(200*(s-1)+1:200*s);y=fft(R,N);c(s,:)=abs(y(1:64));%幅頻譜取64點(diǎn)，存入。r(s,:)=c(s,:); %r=cz=find(c(s,:)<40);消除頻譜泄漏現(xiàn)象(閥值=40)，結(jié)果再存入。c(s,z)=zeros(size(z));endsm=[3144;2340;2344;2348;2640;2644;2648;2840;2844;2848];數(shù)字鍵0-9對(duì)應(yīng)的頻率點(diǎn)表矩陣smfori3=1:8*nb=nnz(c(i3,:)); 命令nnz 確定數(shù)組中N0數(shù)據(jù)的個(gè)數(shù)；ifb==2 若b=2,則c為信號(hào)幅頻譜，其N0的下標(biāo)q1即為頻率點(diǎn)。q1=find(c(i3,:));fori4=1:10 查表矩陣sm,將q1還原成相應(yīng)的數(shù)字鍵，存Aanifq1==sm(i4,:)AN(i3)=i4-1;break;endendelseAN(i3)=NaN; 若%)^2，則c為間隔時(shí)間，^^AN=NaN(空信號(hào)標(biāo)志)endendAN=AN 顯示解碼結(jié)果AN程序結(jié)果：解碼：詢問 輸入n= 1 (一個(gè)DTMF信號(hào))n0= 200 (信號(hào)前有200點(diǎn)空信號(hào))0-9 1 (輸入數(shù)字鍵“1”)顯示：AN=NaN1111NaNNaNNaN(8幀信號(hào))上一行中，第一個(gè)“NaN”表示接收到一個(gè)空幀(200點(diǎn))，后4幀均為“1”，是解碼得到的結(jié)果，與按下的數(shù)字鍵相符，且表示信號(hào)持續(xù)時(shí)間為200X4=800點(diǎn)，最后3幀為間隔時(shí)間，因?yàn)榈谝粠傻搅丝招盘?hào)，所以這里只有三幀。圖形A接收到的DTMF信號(hào)同圖1B輸入plot(r(2,:)),得圖2

5040504030201050“0 10 20 30 40 50 60 704030圖2DTMF信號(hào)作FFT所得的幅頻譜圖圖2為數(shù)字鍵0“1”對(duì)應(yīng)DTMF信號(hào)的幅頻譜，可見，它有兩條譜線，在頻率點(diǎn)23、40處。但周圍發(fā)生了頻譜泄漏現(xiàn)象。C輸入plot(ft(2,:)),得消除頻譜泄漏后的幅頻譜圖(圖3)2001010203040506070504001101201130*1501607030一20一10一1111.II"l1111-0010203040506070圖3消除頻譜泄漏后的幅頻譜圖二、 有限長(zhǎng)序列離散傅里葉變換(DFT)算法用FFT算法解碼每幀信號(hào)共涉及256個(gè)頻率分量，故每幀信號(hào)要算N=256點(diǎn)FFT,但實(shí)際上，組成所有DTMF信號(hào)只用到8個(gè)頻率分量(見表1的fL和fH)，于是，可直接利用DFT定義式進(jìn)行頻譜分析，且每幀信號(hào)只算8點(diǎn)DFT,以避開FfT中許多無意義的計(jì)算，且同樣達(dá)到解碼的目的。DFT算法解碼過程如下：接收DTMF信號(hào)，并畫圖。對(duì)信號(hào)作DFT(每幀信號(hào)作8點(diǎn))，畫頻譜圖，從中找出代表各信號(hào)的特征字。長(zhǎng)度為N的序列x(n)，其DFT仍是一個(gè)長(zhǎng)度為N的序列X(k),它們的關(guān)系是：X(k)=DFT[x(n)]=Tx(n)Wknk=0,1,...N-1(3)n=0…. ,2K 其中WN=Cn稱為旋轉(zhuǎn)因子，則% _另N k k這里，f 頻率點(diǎn)k對(duì)應(yīng)的實(shí)際頻率(Hz)kfs 取樣頻率(Hz)Wk=eJn=efN k k這里，f 頻率點(diǎn)k對(duì)應(yīng)的實(shí)際頻率(Hz)kfs 取樣頻率(Hz)① k對(duì)應(yīng)的數(shù)字頻率(^刁)kkfskf 2兀k 2兀f證:f=kF=N=矣,即一^=fk=3k證:??.(4)式成立。??.(4)式成立。⑶式具有相關(guān)意義，即當(dāng)x(n)中包含頻率為［的分量時(shí)，則x(n)在fk上DFT的幅值|X(k)|與它在其它頻率上DFT的幅值相比，是最大的。于是，以200點(diǎn)為一幀，對(duì)x(n)在8個(gè)特定頻率(見表1的fL和fH)上作DFT，并畫幅頻譜圖，從中找出幅值最大的兩條譜線，在解決頻譜泄漏現(xiàn)象之后，這兩條譜線的橫坐標(biāo)就是代表各信號(hào)的特征字。查表3，將各DTMF信號(hào)還原為相應(yīng)的數(shù)字鍵。56711232456378940表3數(shù)字鍵對(duì)應(yīng)特征字組合表相關(guān)程序如下：A=wavread('D2.wav');subplot(212);plot(A);%畫出接收到的DTMF信號(hào)w=[6977708529411209133614771633];a1=2*pi/8000;w=a1*w;%DTMF>號(hào)用到的八個(gè)數(shù)字頻率存入加數(shù)組forl1=1:8*n 每200點(diǎn)為一幀，在8個(gè)特定頻率上作DFT，幅頻譜數(shù)據(jù)存入「fork1=1:8s1=0;m=0:199;d=cos(w(k1)*m)-j*sin(w(k1)*m);a2=out((l1-1)*200+m+1);s1=sum(d.*a2);r(l1,k1)=abs(s1);endc(l1,:)=r(l1,:);z=find(c(l1,:)<40); %消除頻譜泄漏現(xiàn)象，數(shù)據(jù)存入。c(l1,z)=zeros(size(z));endsm=[46;15;16;17;25;26;27;35;36;37];%敏字鍵0-9對(duì)應(yīng)特征字表矩陣fori3=1:8*n 壹sm,將特征字還原為相應(yīng)的數(shù)字鍵,過程與FFT法相似。b=nnz(c(i3,:));ifb==2q1=find(c(i3,:));fori4=1:10ifq1==sm(i4,:)AN(i3)=i4-1;break;endendelseAN(i3)=NaN;endend程序結(jié)果:

0.501）解碼輸入和同F(xiàn)0而算法所得結(jié)果。200 400600800 1000 12001400 1600 18005040302010200 400600800 1000 12001400 1600 18005040302010圖4DTMF信號(hào)作8點(diǎn)DFT所得幅頻譜圖即為數(shù)字鍵圖4為數(shù)字鍵“1”對(duì)應(yīng)DTMF信號(hào)的8點(diǎn)DFT的幅頻譜。兩條譜線的橫坐標(biāo)“1、5”“1”的特征字。即為數(shù)字鍵程序結(jié)果:

0.50解碼輸入和同F(xiàn)0而算法所得結(jié)果。4）圖形-10200 400600800 1000 12001400 1600 180050403020100 14）圖形-10200 400600800 1000 12001400 1600 180050403020100 1234圖4DTMF信號(hào)作8點(diǎn)DFT所得幅頻譜圖圖4為數(shù)字鍵“1”對(duì)應(yīng)DTMF信號(hào)的8點(diǎn)DFT的幅頻譜。兩條譜線的橫坐標(biāo)“1、5”即為數(shù)字鍵“1”的特征字。三、卷積法當(dāng)用DFT算法求x(n)的頻譜X(k)時(shí)，我們發(fā)現(xiàn)這也是一卷積過程，下面就具體說明怎樣用卷積法解碼。接收DTMF信號(hào)，并繪圖。用卷積法檢測(cè)DTMF信號(hào)，找出代表各信號(hào)的特征字。(一) 卷積法解碼的框圖如下圖所示：—k(n)=W-knu(n)J >y(n)=x(n)*h(n)—n=N>X(k)圖5卷積法解碼框圖1圖5的方框是一個(gè)離散時(shí)間系統(tǒng)，廣義來說，它就是一個(gè)數(shù)字濾波器，其單位脈沖響應(yīng)為:hk(n)=Wnknu(n),長(zhǎng)度為N的因果序列x(n)通過該系統(tǒng)后，輸出為yk(n)，而x(n)在某一頻率點(diǎn)k上的DFT:X(k)=yk(N) ⑸證：由于線性時(shí)不變系統(tǒng)對(duì)任意輸入信號(hào)的響應(yīng)綜上所述，X(k)就是x(n)與單位脈沖響應(yīng)為hk(n)=WN-knu(n)的濾波器(系統(tǒng))的卷積在

n=N上的輸出值。而對(duì)不同的k,hk（n）就不同。于是本文中k取DTMF信號(hào)用到的8個(gè)頻率點(diǎn)，且每次卷積x（n）和hk（n）均取N=200點(diǎn)，余下的解碼過程同DFT算法。另外，卷積在MATLAB中用“。時(shí)^函數(shù)實(shí)現(xiàn)。"（二）程序?qū)崿F(xiàn)中的問題由于序列x（n）下標(biāo)從“0”開始，而MATLAB中數(shù)組下標(biāo)從“1”開始，為了編程方便，于是提出： 當(dāng)N>>1時(shí)，X（k）=yk（N）eyk（N-1） （6）先構(gòu)造一濾波器，其單位脈沖響應(yīng)為：h(n)=<W(N-1-n)kn=0,1...N-1h(n)=<0其它則卷積法框圖變?yōu)椋骸獂(n) >h(n)牝h (n)] >j(n)牝 j (n)—n=N 1>X(k) =y(N一1)牝j (N一1) (N>>1)k 圖k6卷積法框圖2 kX(k)=祝x(m)WNm=x(0)嗎+x⑴嗎+...+x(N-1)Wn(n一i)m=0又j(n)=x(n)*h(n)=Sx(m)h(n-m)m=0y(N-1)=Sx(m)h(N-1-m)=x(0)h(N-1)+x⑴仞N一2)+...+x(N-1)h(0)m=0=x(0)W[N-1-(N-1)]k+x(1)W[N-1-(N-2)]k+...+x(N-1)W(N-1)k=NN Nx(0)W0+x⑴Wk+...+x(N-1)Wk(N-1),X(k)=y(N-1)NN N?.?當(dāng)N>>1時(shí)，h(n)=W(N-1-n)k牝W(N-n)k=WNk?W-nk=1?W-nk=W-nk=h(n)此時(shí)，x(n)通過h(n)相當(dāng)于通過七(n)，輸出y(n)牝y^(n),則X(k)=y(N-1)注y(N-1)(N>>1),證明成立。相關(guān)程序如下：N=200; 接收DTMF信號(hào)，并繪圖A=wavread('D2.wav');subplot(212);plot(A);w=[6977708529411209133614771633];a1=2*pi/8000;w=a1*w; %w——DTfS號(hào)用到的數(shù)字頻率forl1=1:8*n 每200點(diǎn)為一幀，在8個(gè)特定頻率w(1)?w(8)上作卷積y=x1*hkfork1=1:8 %x--信號(hào)；h 濾波器(h=cos(Wm)+jsin(3m))1 k k k1 k1x1=out((l1-1)*200+1:l1*200);m=0:1:N-1;hk=cos(w(k1)*m)+j*sin(w(k1)*m);y=conv(x1,hk);r(l1,k1)=abs(y(N));將X]幅頻譜|X(、)|^|y(N-1)|存入工數(shù)組endc(l1,:)=r(l1,:);z=find(c(l1,:)<40); 消除頻譜泄漏現(xiàn)象c(l1,z)=zeros(size(z));end查表3，將各DTMF信號(hào)還原為相應(yīng)的數(shù)字鍵。

程序同DFT法相應(yīng)部分。卷積法程序結(jié)果與DFT算法的結(jié)果相同。四、迭代法為使解碼過程更接近于硬件的實(shí)現(xiàn)，可在卷積法的基礎(chǔ)上，找到hk(n)滿足的差分方程，這樣，以接收到的DTMF信號(hào)x(n)為輸入，用迭代法解此差分方程，則在n=N時(shí)刻的輸出yk(n)就是X(k),再畫|X(k)|2,即可對(duì)信號(hào)進(jìn)行頻譜分析，從而解碼。迭代法解碼過程如下:接收DTMF信號(hào)，并繪圖。用迭代法對(duì)DTMF信號(hào)進(jìn)行檢測(cè)，找出代表各信號(hào)的特征字。(一) 求hk(n)滿足的差分方程。1)由、(n)—Hk(z)——七(n)的差分方程由卷積法得h(n)=W-knu(n)，則其z變換為：H(z)^^W~knZ-nJE(W-kz-1)nn—0 n—0二——1———^^，-.(1-W-kz-1)Y(z)—X(z),即Y(z)—W-kz-1Y(z)—X(z),1一W-kz-1X(z) Nk kNkN對(duì)上式兩邊進(jìn)行z反變換，得七(n)滿足的差分方程為：*(n)-W-ky^(n-1)—x(n)即：(7)式對(duì)應(yīng)的系統(tǒng)結(jié)構(gòu)圖如下:yk(n)=WN-kyk(n-1)+x(n) (7)即：(7)式對(duì)應(yīng)的系統(tǒng)結(jié)構(gòu)圖如下:因此在計(jì)算yk(n)因此在計(jì)算yk(n)時(shí)，會(huì)碰到復(fù)數(shù)運(yùn)算，不方便，為避免復(fù)數(shù)運(yùn)Hk(z)—1項(xiàng)：Hk(z)分母有理化 K*kl djWN「kW-k

N則yk(n)ln=N=X(k)。2)Hk(z)分母有理化。(7)式中，由于WN-k是復(fù)數(shù)算，首先可對(duì)Hk(z)進(jìn)行分母有理化，其基本原理為：用一對(duì)復(fù)共軛極點(diǎn)代替單極點(diǎn)濾波器的思想。則：H(z)- 1 -―上M——―匕M—k 1-WNkZ-1 (1-WNkZ-1)(1-Wkz-1)1-(WNk+Wjk)z-1+Z-2.2n 2n1入,2兀一入W-k+Wk=ejNk+ejNk=2cos(^k)=2cosro,/.1-Wkz-1 N (8)1-2cosroz-1+z-2Hk(z)分母有理化前后在z平面上的零極點(diǎn)分布如下圖所示:d2=WNk=C2(9)即：階系統(tǒng)結(jié)構(gòu)圖，可通過正則化，將它變換為只有兩個(gè)延時(shí)環(huán)節(jié)的標(biāo)準(zhǔn)形式(稱為正準(zhǔn)型)。正則化步驟如下：(9)即：階系統(tǒng)結(jié)構(gòu)圖，可通過正則化，將它變換為只有兩個(gè)延時(shí)環(huán)節(jié)的標(biāo)準(zhǔn)形式(稱為正準(zhǔn)型)。正則化步驟如下：圖9Hk(z)零極點(diǎn)分布圖從圖9中可得：Hk(z)原來只有一個(gè)極點(diǎn)d1=WN-k(在單位圓上)，經(jīng)分母有理化后，Hk(z)的極點(diǎn)變?yōu)?個(gè)，即d=W-k、d=Wk,且又出現(xiàn)一個(gè)零點(diǎn)c=Wk，因?yàn)閐和c相互抵消，所以H(z)TOC\o"1-5"\h\z1N 2N 2N 2 2 k的表達(dá)式不變，且Hk(z)分母中沒有出現(xiàn)復(fù)數(shù)。由(8)式得：H(z)=——1—WNzT =Xk^,(1-2cos?z-1+z-2)Y(Z)二(1—Wkz-1)X(z)k 1一2cos?z-1+z-2 X(z) k k N對(duì)上式兩邊進(jìn)行z反變換，得)(n)-2cos?y(n-1)+y(n-2)=x(n)-Wkx(n-1)y(n)=x(n)-Wkx(n-1)+2cos?y(n-1)-y(n-2)k N kk k(8)式中，令丑(z)= 1 =^^(8')\o"CurrentDocument"1 1-2coswkz-1+z-2 X(z)H(z)=1-Wkz-1= (8''),則H(z)=H(z)?H(z)\o"CurrentDocument"2 N V(z) k1 2則戲(z)滿足的差分方程為v(n)=x(n)+2coswkv(n-1)-v(n-2)(9')H(z)滿足的差分方程為y(n)=v(n)-Wkv(n-1)(9'')由卷積定理：頻域相乘對(duì)應(yīng)于時(shí)域相卷得：h(n)=h(n)*h(n)(10)k 1 2因此可將h(n)拆成h(n)和h(n)相串聯(lián)的形式，再合并相同的單位延時(shí)環(huán)節(jié)，k 1 2x(n)v(n)yk(n)x(n)v(n)yk(n)圖11hk(n)系統(tǒng)結(jié)構(gòu)圖3(正準(zhǔn)型)Z-1TOC\o"1-5"\h\z顯然，圖11是正準(zhǔn)型的二階系統(tǒng)結(jié)構(gòu)圖。綜上所述，hk(n)滿足的差分方程為(9’)和(9’’)式，即： "(v(n)=x(n)+2coscov(n-1)-v(n-2)(9')1Iy(n)=v(n)-Wkv(n-1) (9'')I\o"CurrentDocument"k N J且x(k)=y(N)=v(N)-Wkv(N-1) (11)kN(二)畫|X(k)|2,對(duì)信號(hào)進(jìn)行頻譜分析。(9’)式中，n=0,1...N-1,共進(jìn)行N次實(shí)數(shù)乘法，而(9’’)中，僅在n=N時(shí)刻計(jì)算1次乘法，但wnk還是復(fù)數(shù)。由于解碼算法中，不關(guān)心X(k),只要求出|X(k)|或|X(k)|2即可進(jìn)行頻譜分析，這里，算|X(k)|2。由式(11)得，IX(k)|2=|v(N一Wkv(N一1)|2(其中Wk=cos①一jsin①)=Iv(N)一cos①v(N一1)+jsin①v(N一1)|2=[v(N)一cos①v(