双音多频通信设计的matlab仿真.doc

双音多频（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的两倍，即 fs≥2fc (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('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; %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; endend 2 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/N≈31.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]; 。