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

双音多频（DTMF）通信设计的MATLAB仿真摘要：讨论以MATLAB作为仿真工具产生DTMF信号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码。关键词：FFT；DFT；频谱分析；卷积；滤波；差分方程；MATLAB0 引 言 双音多频(DTMF:Double Tone Multi-Frequency)是按键电话通信，也广泛应用于电子邮件和银行系统中。用户可从电话发送DTMF信号来选择菜单进行操作。DTMF信号容易用软件产生和解码。MATLAB是一个高度集成的软件系统，通过交互式的命令（语句）可以十分简便地实现许多复杂的数值计算。 本文采用MATLAB作为仿真工具产生DTMF信号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码，由此得出：时域和频域是研究信号的两个窗口，其中信号处理大都在时域中进行，而信号分析往往在频域中比较方便直观。且数字信号处理技术中的DFT、FFT、卷积、滤波、差分方程这几个概念之间有内在联系。1 DTMF信号的产生DTMF是数字音频信号，在DTMF通信系统中共有8个频率，分为4个高频音和4个低频音，用一个高频音和一个低频音的组合表示一个信号，这样共有16种组合，分别代表16种信号，如表1所示： fH(Hz)fL(Hz) 1209 1336 1477 1633 697 1 2 3 A 770 4 5 6 B 852 7 8 9 C 941 * 0 # D 表1 DTMF信号组合表例如，当按下数字键“1”时，则产生低频697Hz和高频1209Hz这两个正弦信号的迭加。由于语音信号的最高频率为4KHz,根据奈奎斯特取样定理，取样频率fs应大于或等于原信号最高频率fc的两倍，即 fs2fc (1)才能保证取样后的信号不失真，所以电话音频信号在数字信号处理时，取样频率fs为2×4k=8kHz,这里，每个数字信号持续时间为100ms,后面加上100ms的间隔时间（用0表示）。上述DTMF信号产生方法如下：（1） 建立拨号数字的表矩阵，用查表法（查表1）求用户所按数字键对应的高、低频音。为简化起见，仅允许选择“0-9”这十个键，在开始时还可拨空信号。（2） 产生相应的DTMF信号及间隔时间。由于fs=8kHz,各信号持续时间为100ms,因此在程序中每个信号取800点，间隔时间也取800点，结果存入数组中。（3）画图并监听产生的DTMF信号。程序如下：%程序1， DTMF信号的产生clear %清除内存TAB=941 1336;697 1209;697 1336;697 1477;770 1209;770 1336;770 1477;852 1209;852 1336;852 1477;%拨号数字表矩阵n=input('n='); %DTMF信号的个数l=input('n0=');%空信号点数for i=1:nk=input('09');%输入的数字键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; %out数组中每个数据除以2subplot(211); plot(out); %绘图并监听DTMF信号。sound(out,fs)wavwrite(out,fs,'D2.wav');%out另存入声音文件程序结果：例如：若要产生数字键“1”所对应的DTMF信号，且开始时有200点空信号，则程序运行如下： 询问： 输入： n= 1 n0= 2000-9 1得图1： 图1 DTMF信号图图1中，前200点为空信号，第201-1000点为键“1”对应的DTMF信号，第1001-1800点为间隔时间。另外，在程序中为使软件设计更接近于实际硬件的开发应用，可用求解差分方程的方法来代替正弦函数的调用。设正弦序列为h(n)=sin(kn)u(n),为实时实现h(n)，必须找到其满足的差分方程。正弦序列h(n)的z变换为：两边进行反变换，得 y(n)-ay(n-1)+y(n-2)=bx(n-1),式中，若令x(n)=(n)，则得到h(n)的差分方程为h(n)-ah(n-1)+h(n-2)=b(n-1),即： h(n)= ah(n-1)-h(n-2)+b(n-1) （2）用迭代法解此差分方程，即得数字频率为k的正弦序列h(n)。在硬件中，该差分方程是由加法器、乘法器和单位延时单元构成的系统。将x(n)=(n)输入该系统后，输出的就是h(n)。本文中，每个DTMF信号h(n)是两个频率的正弦序列相迭加，设为hL(n)和hH(n),为此，分别求得hL(n)和hH(n)所满足的差分方程：hL(n)= aLhL(n-1)-hL(n-2)+bL(n-1) ；hH(n)= aHhH(n-1)-hH(n-2)+bH(n-1)则h(n)= hL(n)+hH(n) 。相关程序如下：%程序2，用解差分方程方法产生DTMF信号。x=zeros(1,800);x(2)=1;%x(n)=(n-1)(x(n)取800点)n=input('n=')； %DTMF信号的个数l=input('n0=')； %空信号的点数out=zeros(1,1600*n+l);%out数组初始化w=2*pi/8000*941 1336;697 1209;697 1336;697 1477; 770 1209;770 1336;770 1477;852 1209;852 1336;852 1477;%各信号对应的数字频率tab=2*cos(w) sin(w);%各信号满足的差分方程hL(n)和hH(n)中系数a,b的矩阵。 for i=1:n k=input('0-9-'); %输入的数字键 hL=zeros(1,3);hH=zeros(1,3);% hL(n)和hH(n)初始化 for j=1:800 %迭代法求解hL(n)和hH(n)产生相应的DTMF信号，存入out hL(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; endend2 DTMF信号的解码 解码就是对接收到的DTMF信号进行频谱分析，从中找出代表各信号的的特征字，由此获知用户按下的数字键。为在频谱图中分辨出不同的频率分量，于是对信号取200点为一帧，则频谱分辨率F=fs/N=8000/200=40Hz<73Hz(表1中任意两频率的最小间隔)，这样，即可满足频谱分析的要求。 本文采用数字信号处理技术中的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/N31.25Hz。因为信号x为实数序列，所其幅频谱|y|具有偶对称性，于是，幅频谱可以仅画N/2点，其中第N/2点对应实际频率为fs/2=4KHz,而DTMF信号中最高频率为1633Hz，小于2KHz（fs/4），因此，这里只画N/4=64点。DTMF信号是两个正弦波的迭加，它的幅频谱就是两根谱线，谱线的横坐标就是该信号的两个频率分量点KL和KH 。2） 消除频谱泄漏现象。由于信号x是有限长的，这就相当于对无限长的信号加矩形窗，所以在频谱图中必然会出现频谱泄漏现象,使信号能量散布到其他谱线位置。为此，应选择一适当阀值，将出现在这两条谱线周围的幅度较小的谱线消除（置0），从而解决了这一问题。最后，将处理后的幅频谱数据存入数组c中。3） 将各DTMF信号还原为相应的数字键。在幅频谱图中，频率轴的定标方式为频率点K而不是实际频率f,转换关系为：K=f/F,因此，数字键0-9对应频率点如下表所示： KH KL 40 44 48 23 1 2 3 26 4 5 6 28 7 8 9 31 0 表 2 数字键对应频率点组合表数组c的不等于0的下标就是各信号的频率点，查表2，即可将各DTMF信号还原为相应的数字键。相关程序如下：%程序3，FFT算法解码程序。A=wavread('D2.wav'); %接收到的DTMF信号subplot(212); plot(A);%绘图N=256;for s=1:8*n %对每帧信号作N=256点的FFT R=out(200*(s-1)+1:200*s); y=fft(R,N); c(s,:)=abs(y(1:64);%幅频谱取64点，存入c r(s,:)=c(s,:); %r=c z=find(c(s,:)<40); %消除频谱泄漏现象（阀值=40），结果再存入c c(s,z)=zeros(size(z);endsm=31 44;23 40;23 44;23 48;26 40;26 44;26 48;28 40;28 44;28 48; %数字键0-9对应的频率点表矩阵sm for i3=1:8*n b=nnz(c(i3,:); %命令nnz-确定数组中0数据的个数； if b=2 %若b=2,则c为信号幅频谱，其0的下标q1即为频率点。 q1=find(c(i3,:); for i4=1:10 %查表矩阵sm，将q1还原成相应的数字键，存入AN if q1=sm(i4,:) AN(i3)=i4-1;break; end end else AN(i3)=NaN; %若b2，则c为间隔时间，则AN=NaN(空信号标志) endendAN=AN %显示解码结果AN程序结果：1) 解码： 询问 输入 n= 1 （一个DTMF信号） n0= 200（信号前有200点空信号）0-9 1 （输入数字键“1”） 显示：AN=NaN 1 1 1 1 NaN NaN NaN （8帧信号）上一行中，第一个“NaN”表示接收到一个空帧（200点），后4帧均为“1”，是解码得到的结果，与按下的数字键相符,且表示信号持续时间为200×4=800点,最后3帧为间隔时间，因为第一帧采到了空信号，所以这里只有三帧。2）图形A 接收到的DTMF信号同图1B 输入 plot(r(2,:),得图2 图2 DTMF信号作FFT所得的幅频谱图图2为数字键“1”对应DTMF信号的幅频谱，可见，它有两条谱线，在频率点23、40处。但周围发生了频谱泄漏现象。C 输入 plot(c(2,:),得消除频谱泄漏后的幅频谱图（图3） 图3 消除频谱泄漏后的幅频谱图二、 有限长序列离散傅里叶变换（DFT）算法 用FFT算法解码每帧信号共涉及256个频率分量，故每帧信号要算N=256点FFT,但实际上，组成所有DTMF信号只用到8个频率分量（见表1的fL和fH）,于是，可直接利用DFT定义式进行频谱分析，且每帧信号只算8点DFT,以避开FFT中许多无意义的计算，且同样达到解码的目的。DFT算法解码过程如下：（1） 接收DTMF信号，并画图。（2） 对信号作DFT（每帧信号作8点），画频谱图，从中找出代表各信号的特征字。长度为N的序列x(n)，其DFT仍是一个长度为N的序列X(k),它们的关系是：(3) 式具有相关意义，即当x(n)中包含频率为fk的分量时，则x(n)在fk上DFT的幅值|X(k)|与它在其它频率上DFT的幅值相比，是最大的。于是，以200点为一帧，对x(n)在8个特定频率（见表1的fL和fH）上作DFT,并画幅频谱图，从中找出幅值最大的两条谱线，在解决频谱泄漏现象之后，这两条谱线的横坐标就是代表各信号的特征字。（3） 查表3，将各DTMF信号还原为相应的数字键。 KH KL 5 6 7 1 1 2 3 2 4 5 6 3 7 8 9 4 0 表3 数字键对应特征字组合表相关程序如下：A=wavread('D2.wav');subplot(212); plot(A);%画出接收到的DTMF信号w=697 770 852 941 1209 1336 1477 1633;a1=2*pi/8000;w=a1*w;%DTMF信号用到的八个数字频率存入w数组for l1=1:8*n %每200点为一帧，在8个特定频率上作DFT,幅频谱数据存入r for k1=1:8 s1=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); end c(l1,:)=r(l1,:);z=find(c(l1,:)<40); %消除频谱泄漏现象，数据存入cc(l1,z)=zeros(size(z);endsm=4 6;1 5;1 6;1 7;2 5;2 6;2 7;3 5;3 6;3 7; %sm-数字键0-9对应特征字表矩阵for i3=1:8*n %查sm,将特征字还原为相应的数字键,过程与FFT法相似。 b=nnz(c(i3,:); if b=2 q1=find(c(i3,:); for i4=1:10 if q1=sm(i4,:) AN(i3)=i4-1;break; end end else AN(i3)=NaN; endend程序结果：1) 解码输入和同FFT算法所得结果。2) 图形 图4 DTMF信号作8点DFT所得幅频谱图图4为数字键“1”对应DTMF信号的8点DFT的幅频谱。两条谱线的横坐标“1、5”即为数字键“1”的特征字。程序结果：3) 解码输入和同FFT算法所得结果。4) 图形 图4 DTMF信号作8点DFT所得幅频谱图图4为数字键“1”对应DTMF信号的8点DFT的幅频谱。两条谱线的横坐标“1、5”即为数字键“1”的特征字。三、卷积法当用DFT算法求x(n)的频谱X(k)时，我们发现这也是一卷积过程，下面就具体说明怎样用卷积法解码。(1)接收DTMF信号，并绘图。(2)用卷积法检测DTMF信号，找出代表各信号的特征字。（一） 卷积法解码的框图如下图所示： 图 5 卷积法解码框图1图5的方框是一个离散时间系统，广义来说，它就是一个数字滤波器，其单位脉冲响应为：,长度为N的因果序列x(n)通过该系统后，输出为yk(n)，而x(n)在某一频率点k上的DFT: X(k)=yk(N) (5)证：由于线性时不变系统对任意输入信号的响应综上所述，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中用“conv”函数实现。（二）程序实现中的问题由于序列x(n)下标从“0”开始，而MATLAB中数组下标从“1”开始，为了编程方便，于是提出： 当N>>1时，X(k)= yk(N)yk(N-1) (6)先构造一滤波器，其单位脉冲响应为： 图 6 卷积法框图2相关程序如下：N=200; %接收DTMF信号，并绘图A=wavread('D2.wav');subplot(212); plot(A);w=697 770 852 941 1209 1336 1477 1633;a1=2*pi/8000;w=a1*w; %w-DTMF信号用到的数字频率for l1=1:8*n %每200点为一帧，在8个特定频率w(1)w(8)上作卷积y=x1*hk for k1=1:8 %x1-信号；hk-滤波器（hk=cos(k1m)+jsin(k1m)） x1=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);%将x1幅频谱|X(k1)|y(N-1)|存入r数组 end c(l1,:)=r(l1,:); z=find(c(l1,:)<40); %消除频谱泄漏现象 c(l1,z)=zeros(size(z);end（3）查表3，将各DTMF信号还原为相应的数字键。程序同DFT法相应部分。卷积法程序结果与DFT算法的结果相同。四、迭代法为使解码过程更接近于硬件的实现，可在卷积法的基础上，找到hk(n)满足的差分方程，这样，以接收到的DTMF信号x(n)为输入，用迭代法解此差分方程，则在n=N时刻的输出yk(n)就是X(k),再画|X(k)|2,即可对信号进行频谱分析，从而解码。迭代法解码过程如下：(1) 接收DTMF信号，并绘图。(2) 用迭代法对DTMF信号进行检测，找出代表各信号的特征字。（一） 求hk(n)满足的差分方程。即： yk(n)=WN-kyk(n-1)+x(n) (7)yk(n)z-1x(n)(7)式对应的系统结构图如下：WN-k 图 8 hk(n)系统结构图1则yk(n)|n=N =X(k)。2)Hk(z)分母有理化。(7)式中，由于WN-k是复数，因此在计算yk(n)时，会碰到复数运算，不方便，为避免复数运算，首先可对Hk(z)进行分母有理化，其基本原理为：用一对复共轭极点代替单极点滤波器的思想。d2=WNk=c2d1=WN-kHk(z)分母有理化d1=WN-kHk(z)分母有理化前后在z平面上的零极点分布如下图所示： 图9 Hk(z)零极点分布图从图9中可得：Hk(z)原来只有一个极点d1=WN-k（在单位圆上），经分母有理化后，Hk(z)的极点变为2个，即d1=WN-k、 d2=WNk,且又出现一个零点c2=WNk,因为d2和c2相互抵消，所以Hk(z)的表达式不变，且Hk(z)分母中没有出现复数。由（8）式得：yk(n)此时，hk(n)系统结构图如下：2coskWNK z-1 z-1 z-1 图10 hk(n)系统结构图2式（9）为二阶差分方程，而图10中却有三个单位延时环节，所以图10不是标准形式的二阶系统结构图，可通过正则化，将它变换为只有两个延时环节的标准形式（称为正准型）。正则化步骤如下：WNK2cosk x(n) v(n) yk(n) z-1 z-1 图11 hk(n)系统结构图3（正准型）显然，图11是正准型的二阶系统结构图。综上所述，hk(n)满足的差分方程为 (9)和（9）式，即：（二） 画|X(k)|2,对信号进行频谱分析。 (9)式中，n=0,1N-1,共进行N次实数乘法，而(9)中，仅在n=N时刻计算1次乘法，但WNK还是复数。由于解码算法中，不关心X(k),只要求出|X(k)|或|X(k)|2即可进行频谱分析，这里，算|X(k)|2。由式（11）得，这样解码过程中就完全避免了复数运算。余下解码过程同DFT算法相应部分。相关程序如下：%程序6迭代法程序N=200; %接收DTMF信号，并绘图a=wavread('D2.wav');subplot(212);plot(a);w=697 770 852 941 1209 1336 1477 1633;a1=2*pi/8000;w=a1*w;%w-8个特定的数字频率for l=1:8*n %200点为一帧，在8个特定的数字频率上将x(n)输入hk(n)满足的 for k1=1:8 %差分方程，以迭代法解此方程，求得v(N)和v(N-1) v=zeros(1,3); for m=1:200 v(3)=2*cos(w(k1)*v(2)-v(1)+a(l-1)*200+m); v(1)=v(2); v(2)=v(3); end r(l,k1)=v(2).2+v(1).2-2*cos(w(k1)*v(2)*v(1); %|x(k)|2= end c(l,:)=r(l,:);%消除频谱泄漏现象 q=find(c(l,:)<2000); c(l,q)=zeros(size(q);end(3) 将各DTMF信号还原为相应的数字键。程序同DFT算法。本程序结果同前三种算法。3 四种解码算法的比较下面，以处理1帧信号（200点）为单位，对这四种算法作一比较，如下表所示： 计算量占用内存量计算复杂度FFT法N/2log2N=1024次复数乘法（N=256k取256点）（算|X(k)|）（最少） 256点（|X(k)|）（最大）在MATLAB中极易，但算法不“透明”DFT法N×8=1600次复数乘法（N=200k取8点）（算|X(k)|） 8点（|X(k)|） 较难卷积法同DFT法同DFT法与DFT法相当迭代法(N+4)×8=1632次实数乘法（算|X(k)|2）同DFT法（|X(k)|2）最难，但更贴近硬件实现4 结论 本文采用四种算法对仿真的DTMF信号解码，它们的共同点在于对信号作DFT,从而进行频谱分析。DFT的快速算法为FFT,而DFT又是一种运算，可通过将信号输入一数字滤波器（离散时间系统），取其在某点上的输出实现，上述方法既是一卷积过程，又是一解差分方程的过程。由此得出：DFT、FFT、卷积、滤波、差分方程这几个概念之间有内在联系。 MATLAB的应用，大大减少了计算工作量，且使比较抽象的概念得以直观显示，它不愧为功能强大的计算机辅助教学工具。参考文献1 戴悟僧，冯筱林，安邦健.数字信号处理导论.M.上海：科学技术出版社，2000.2吴湘淇.信号、系统与信号处理（下）.M.北京：电子工业出版社，2000.