双音多频通信设计的Matlab仿真

双音多频(DTMF)通信设计的MATLAB仿真摘要：讨论以MATLAB作为仿真工具产生DTMF信号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码。关键词：FFT；DFT；频谱分析；卷积；滤波；差分方程；MATLAB0引言双音多频(DTMF:DoubleToneMulti-Frequency)是按键电话通信，也广泛应用于电子邮件和银行系统中。用户可从电话发送DTMF信号来选择菜单进行操作。DTMF信号容易用软件产生和解码。MATLAB是一个高度集成的软件系统，通过交互式的命令(语句)可以十分简便地实现许多复杂的数值计算。本文采用MATLAB作为仿真工具产生DTMF信号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码，由此得出：时域和频域是研究信号的两个窗口，其中信号处理大都在时域中进行，而信号分析往往在频域中比较方便直观。且数字信号处理技术中的DFT、FFT、卷积、滤波、差分方程这几个概念之间有内在联系。1DTMF信号的产生DTMF是数字音频信号，在DTMF通信系统中共有8个频率，分为4个高频音和4个低频音，用一个高频音和一个低频音的组合表示一个信号，这样共有16种组合，分别代表16种信号，如表1所示：7^)f(Hz)\1209133614771633L697123A770456B852789C941*0#D表1DTMF信号组合表例如，当按下数字键“1”时，则产生低频697Hz和高频1209Hz这两个正弦信号的迭加。由于语音信号的最高频率为4KHz，根据奈奎斯特取样定理，取样频率fs应大于或等于原信号最高频率fc的两倍，即fsN2fc (1)才能保证取样后的信号不失真，所以电话音频信号在数字信号处理时，取样频率fs为2X4k=8kHz,这里，每个数字信号持续时间为100ms,后面加上100ms的间隔时间(用0表示)。上述DTMF信号产生方法如下：建立拨号数字的表矩阵，用查表法(查表1)求用户所按数字键对应的高、低频音。为简化起见，仅允许选择“0-9”这十个键，在开始时还可拨空信号。产生相应的DTMF信号及间隔时间。由于fs=8kHz,各信号持续时间为100ms,因此在程序中每个信号取800点，间隔时间也取800点，结果存入数组中。画图并监听产生的DTMF信号。程序如下：%[程序1],DTMF信号的产生clear 清%除内存TAB=[9411336;6971209;6971336;6971477;7701209;7701336;7701477;8521209;8521336;8521477];%拨号数字表矩阵n=input('n=');%DTMF信号的个数l=input('n0=');%空信号点数fori=1:nk=input('0〜9');%输入的数字键fL=TAB(k+1,1);%k对应的低频音fH=TAB(k+1,2);%k对应的高频音n1=800;fs=8000;%产生k对应的DTMF信号，取样频率8kHz,每个信号共1600点，其中，j=0:1:n1-1; 前800点为信号持续时间，后800点为间隔时间，结果存入数组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; %数组中每个数据除以2subplot(211);plot(out);%绘图并监听DTMF信号。sound(out,fs)wavwrite(out,fs,'D2.wav');%out另存入声音文件程序结果：例如：若要产生数字键“1”所对应的DTMF信号，且开始时有200点空信号，则程序运行如下：询问：输入：n=1n0= 2000-9 1得图1：图1DTMF信号图图1中，前200点为空信号，第201-1000点为键“1”对应的DTMF信号，第1001-1800点为间隔时间0.5。-0.5-10 200 400 600 800 1000 12001400 1600 1800另外，在程序中为使软件设计更接近于实际硬件的开发应用，可用求解差分方程的方法来代替正弦函数的调用。设正弦序列为h(n)=sin(3kn)u(n),为实时实现h(n)，必须找到其满足的差分方程。正弦序列h(n)的z变换为：H(z)= zsi「3k ,令H(z)分子/分母系数：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两边进行反变换，得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)用迭代法解此差分方程，即得数字频率为3k的正弦序列h(n)。在硬件中，该差分方程是由加法器、乘法器和单位延时单元构成的系统。将x(n)=6(n)输入该系统后，输出的就是h(n)。本文中，每个DTMF信号h(n)是两个频率的正弦序列相迭加，设为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)。相关程序如下：%[程序2]，用解差分方程方法产生DTMF信号。x=zeros(1,800);x(2)=1;%x(n)=6(n-1)(x(n)取800点)n=input('n=')； %DTM信号的个数l=input('n0=')； 空信号的点数out=zeros(1,1600*n+l);%out数组初始化w=2*pi/8000*[9411336;6971209;6971336;6971477;7701209;7701336;7701477;8521209;8521336;8521477];%各信号对应的数字频率tab=[2*cos(w)sin(w)];%各信号满足的差分方程h〔(n)和hH(n)中系数a,b的矩阵。fori=1:nk=input'0-9-'); 输入的数字键hL=zeros(1,3);hH=zeros(1,3);%]h(n)和hH(n)初始化forj=1:800 迭代法求解hL(n)和hH(n)产生相应的DTMF信号，存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信号的解码解码就是对接收到的DTMF信号进行频谱分析，从中找出代表各信号的的特征字，由此获知用户按下的数字键。为在频谱图中分辨出不同的频率分量，于是对信号取200点为一帧，则频谱分辨率F=f/N=8000/200=40Hz<73Hz（表1中任意两频率的最小间隔），这样，即可满足s频谱分析的要求。本文采用数字信号处理技术中的FFT算法、DFT算法、卷积法及迭代法这4种算法实现对DTMF信号的解码。一、 快速傅里叶变换（FFT）算法FFT是有限长序列离散傅里叶变换（DFT）的快速算法，其基本运算是蝶形算法，它使DFT计算时间缩短了几个数量级，在信号处理中占有极重要的地位。这里采用基2-FFT算法对DTMF信号进行频谱分析。解码过程如下：（1） 接收DTMF信号，并画图。（2） 对信号作FFT，画频谱图，从中找出代表各信号的频率分量。这部分，信号将完成从时域到频域的转换。1） 每帧信号（200点）做一次N=256点的FFT，从中取64点画频谱图。在MATLAB中,FFT可由语句“y=fft（x,N）”来实现。而FFT中，要求序列长度N=2e（E为整数），所以N=28=256，频谱分辨率F=fs/N^31.25Hzo因为信号x为实数序列，所其幅频谱|y|具有偶对称性，于是，幅频谱可以仅画N/2点，其中第N/2点对应实际频率为fs/2=4KHz,而DTMF信号中最高频率为1633Hz，小于2KHz（fs/4），因此，这里只画N/4=64点。DTMF>号是两个正弦波的迭加，它的幅频谱就是两根谱线，谱线的横坐标就是该信号的两个频率分量点kl和kh。2） ；肖除频谱泄漏现象。由于信号x是有限长的，这就相当于对无限长的信号加矩形窗，所以在频谱图中必然会出现频谱泄漏现象，使信号能量散布到其他谱线位置。为此，应选择一适当阀值，将出现在这两条谱线周围的幅度较小的谱线消除（置）），从而解决了这一问题。最后，将处理后的幅频谱数据存入数组c中。3） 将各DTMF信号还原为相应的数字键。在幅频谱图中，频率轴的定标方式为频率点K而不是实际频率f,转换关系为：K=f/F,因此，数字键0-9对应频率点如下表所示："...404448231232645628789310表2数字键对应频率点组合表数组c的不等于0的下标就是各信号的频率点，查表2,即可将各DTMF信号还原为相应的数字键。相关程序如下：%[程序3],FFT算法解码程序。A=wavread('D2.wav');%接收到的DTMF信号subplot(212);plot(A);%绘图N=256;fors=1:8*n 对%每帧信号作N=256点的FFTR=out(200*(s-1)+1:200*s);y=fft(R,N);c(s,:)=abs(y(1:64));%幅频谱取64点，存入。r(s,:)=c(s,:); %r=cz=find(c(s,:)<40);消除频谱泄漏现象(阀值=40)，结果再存入。c(s,z)=zeros(size(z));endsm=[3144;2340;2344;2348;2640;2644;2648;2840;2844;2848];数字键0-9对应的频率点表矩阵smfori3=1:8*nb=nnz(c(i3,:)); 命令nnz 确定数组中N0数据的个数；ifb==2 若b=2,则c为信号幅频谱，其N0的下标q1即为频率点。q1=find(c(i3,:));fori4=1:10 查表矩阵sm,将q1还原成相应的数字键，存Aanifq1==sm(i4,:)AN(i3)=i4-1;break;endendelseAN(i3)=NaN; 若%)^2，则c为间隔时间，^^AN=NaN(空信号标志)endendAN=AN 显示解码结果AN程序结果：解码：询问 输入n= 1 (一个DTMF信号)n0= 200 (信号前有200点空信号)0-9 1 (输入数字键“1”)显示：AN=NaN1111NaNNaNNaN(8帧信号)上一行中，第一个“NaN”表示接收到一个空帧(200点)，后4帧均为“1”，是解码得到的结果，与按下的数字键相符，且表示信号持续时间为200X4=800点，最后3帧为间隔时间，因为第一帧采到了空信号，所以这里只有三帧。图形A接收到的DTMF信号同图1B输入plot(r(2,:)),得图2

5040504030201050“0 10 20 30 40 50 60 704030图2DTMF信号作FFT所得的幅频谱图图2为数字键0“1”对应DTMF信号的幅频谱，可见，它有两条谱线，在频率点23、40处。但周围发生了频谱泄漏现象。C输入plot(ft(2,:)),得消除频谱泄漏后的幅频谱图(图3)2001010203040506070504001101201130*1501607030一20一10一1111.II"l1111-0010203040506070图3消除频谱泄漏后的幅频谱图二、 有限长序列离散傅里叶变换(DFT)算法用FFT算法解码每帧信号共涉及256个频率分量，故每帧信号要算N=256点FFT,但实际上，组成所有DTMF信号只用到8个频率分量(见表1的fL和fH)，于是，可直接利用DFT定义式进行频谱分析，且每帧信号只算8点DFT,以避开FfT中许多无意义的计算，且同样达到解码的目的。DFT算法解码过程如下：接收DTMF信号，并画图。对信号作DFT(每帧信号作8点)，画频谱图，从中找出代表各信号的特征字。长度为N的序列x(n)，其DFT仍是一个长度为N的序列X(k),它们的关系是：X(k)=DFT[x(n)]=Tx(n)Wknk=0,1,...N-1(3)n=0…. ,2K 其中WN=Cn称为旋转因子，则% _另N k k这里，f 频率点k对应的实际频率(Hz)kfs 取样频率(Hz)Wk=eJn=efN k k这里，f 频率点k对应的实际频率(Hz)kfs 取样频率(Hz)① k对应的数字频率(^刁)kkfskf 2兀k 2兀f证:f=kF=N=矣,即一^=fk=3k证:・•.(4)式成立。・•.(4)式成立。⑶式具有相关意义，即当x(n)中包含频率为［的分量时，则x(n)在fk上DFT的幅值|X(k)|与它在其它频率上DFT的幅值相比，是最大的。于是，以200点为一帧，对x(n)在8个特定频率(见表1的fL和fH)上作DFT，并画幅频谱图，从中找出幅值最大的两条谱线，在解决频谱泄漏现象之后，这两条谱线的横坐标就是代表各信号的特征字。查表3，将各DTMF信号还原为相应的数字键。56711232456378940表3数字键对应特征字组合表相关程序如下：A=wavread('D2.wav');subplot(212);plot(A);%画出接收到的DTMF信号w=[6977708529411209133614771633];a1=2*pi/8000;w=a1*w;%DTMF>号用到的八个数字频率存入加数组forl1=1:8*n 每200点为一帧，在8个特定频率上作DFT，幅频谱数据存入「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); %消除频谱泄漏现象，数据存入。c(l1,z)=zeros(size(z));endsm=[46;15;16;17;25;26;27;35;36;37];%敏字键0-9对应特征字表矩阵fori3=1:8*n 壹sm,将特征字还原为相应的数字键,过程与FFT法相似。b=nnz(c(i3,:));ifb==2q1=find(c(i3,:));fori4=1:10ifq1==sm(i4,:)AN(i3)=i4-1;break;endendelseAN(i3)=NaN;endend程序结果:

0.501）解码输入和同F0而算法所得结果。200 400600800 1000 12001400 1600 18005040302010200 400600800 1000 12001400 1600 18005040302010图4DTMF信号作8点DFT所得幅频谱图即为数字键图4为数字键“1”对应DTMF信号的8点DFT的幅频谱。两条谱线的横坐标“1、5”“1”的特征字。即为数字键程序结果:

0.50解码输入和同F0而算法所得结果。4）图形-10200 400600800 1000 12001400 1600 180050403020100 14）图形-10200 400600800 1000 12001400 1600 180050403020100 1234图4DTMF信号作8点DFT所得幅频谱图图4为数字键“1”对应DTMF信号的8点DFT的幅频谱。两条谱线的横坐标“1、5”即为数字键“1”的特征字。三、卷积法当用DFT算法求x(n)的频谱X(k)时，我们发现这也是一卷积过程，下面就具体说明怎样用卷积法解码。接收DTMF信号，并绘图。用卷积法检测DTMF信号，找出代表各信号的特征字。(一) 卷积法解码的框图如下图所示：—k(n)=W-knu(n)J >y(n)=x(n)*h(n)—n=N>X(k)图5卷积法解码框图1图5的方框是一个离散时间系统，广义来说，它就是一个数字滤波器，其单位脉冲响应为:hk(n)=Wnknu(n),长度为N的因果序列x(n)通过该系统后，输出为yk(n)，而x(n)在某一频率点k上的DFT:X(k)=yk(N) ⑸证：由于线性时不变系统对任意输入信号的响应综上所述，X(k)就是x(n)与单位脉冲响应为hk(n)=WN-knu(n)的滤波器(系统)的卷积在

n=N上的输出值。而对不同的k,hk（n）就不同。于是本文中k取DTMF信号用到的8个频率点，且每次卷积x（n）和hk（n）均取N=200点，余下的解码过程同DFT算法。另外，卷积在MATLAB中用“。时^函数实现。"（二）程序实现中的问题由于序列x（n）下标从“0”开始，而MATLAB中数组下标从“1”开始，为了编程方便，于是提出： 当N>>1时，X（k）=yk（N）eyk（N-1） （6）先构造一滤波器，其单位脉冲响应为：h(n)=<W(N-1-n)kn=0,1...N-1h(n)=<0其它则卷积法框图变为：—x(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•.•当N>>1时，h(n)=W(N-1-n)k牝W(N-n)k=WNk•W-nk=1•W-nk=W-nk=h(n)此时，x(n)通过h(n)相当于通过七(n)，输出y(n)牝y^(n),则X(k)=y(N-1)注y(N-1)(N>>1),证明成立。相关程序如下：N=200; 接收DTMF信号，并绘图A=wavread('D2.wav');subplot(212);plot(A);w=[6977708529411209133614771633];a1=2*pi/8000;w=a1*w; %w——DTfS号用到的数字频率forl1=1:8*n 每200点为一帧，在8个特定频率w(1)〜w(8)上作卷积y=x1*hkfork1=1:8 %x--信号；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)|存入工数组endc(l1,:)=r(l1,:);z=find(c(l1,:)<40); 消除频谱泄漏现象c(l1,z)=zeros(size(z));end查表3，将各DTMF信号还原为相应的数字键。

程序同DFT法相应部分。卷积法程序结果与DFT算法的结果相同。四、迭代法为使解码过程更接近于硬件的实现，可在卷积法的基础上，找到hk(n)满足的差分方程，这样，以接收到的DTMF信号x(n)为输入，用迭代法解此差分方程，则在n=N时刻的输出yk(n)就是X(k),再画|X(k)|2,即可对信号进行频谱分析，从而解码。迭代法解码过程如下:接收DTMF信号，并绘图。用迭代法对DTMF信号进行检测，找出代表各信号的特征字。(一) 求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对上式两边进行z反变换，得七(n)满足的差分方程为：*(n)-W-ky^(n-1)—x(n)即：(7)式对应的系统结构图如下:yk(n)=WN-kyk(n-1)+x(n) (7)即：(7)式对应的系统结构图如下:因此在计算yk(n)因此在计算yk(n)时，会碰到复数运算，不方便，为避免复数运Hk(z)—1项：Hk(z)分母有理化 K*kl djWN「kW-k

N则yk(n)ln=N=X(k)。2)Hk(z)分母有理化。(7)式中，由于WN-k是复数算，首先可对Hk(z)进行分母有理化，其基本原理为：用一对复共轭极点代替单极点滤波器的思想。则：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平面上的零极点分布如下图所示:d2=WNk=C2(9)即：阶系统结构图，可通过正则化，将它变换为只有两个延时环节的标准形式(称为正准型)。正则化步骤如下：(9)即：阶系统结构图，可通过正则化，将它变换为只有两个延时环节的标准形式(称为正准型)。正则化步骤如下：图9Hk(z)零极点分布图从图9中可得：Hk(z)原来只有一个极点d1=WN-k(在单位圆上)，经分母有理化后，Hk(z)的极点变为2个，即d=W-k、d=Wk,且又出现一个零点c=Wk，因为d和c相互抵消，所以H(z)TOC\o"1-5"\h\z1N 2N 2N 2 2 k的表达式不变，且Hk(z)分母中没有出现复数。由(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对上式两边进行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'')由卷积定理：频域相乘对应于时域相卷得：h(n)=h(n)*h(n)(10)k 1 2因此可将h(n)拆成h(n)和h(n)相串联的形式，再合并相同的单位延时环节，k 1 2x(n)v(n)yk(n)x(n)v(n)yk(n)图11hk(n)系统结构图3(正准型)Z-1TOC\o"1-5"\h\z显然，图11是正准型的二阶系统结构图。综上所述，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,对信号进行频谱分析。(9’)式中，n=0,1...N-1,共进行N次实数乘法，而(9’’)中，仅在n=N时刻计算1次乘法，但wnk还是复数。由于解码算法中，不关心X(k),只要求出|X(k)|或|X(k)|2即可进行频谱分析，这里，算|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(