双音多频通信设计的Matlab仿真.docx
《双音多频通信设计的Matlab仿真.docx》由会员分享，可在线阅读，更多相关《双音多频通信设计的Matlab仿真.docx（26页珍藏版）》请在冰点文库上搜索。
双音多频通信设计的Matlab仿真
双音多频（DTMF通信设计的MATLAB仿真
摘要：
讨论以MATLAB乍为仿真工具产生DTMF言号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码。
关键词：
FFT;DFT；频谱分析；卷积；滤波；差分方程；MATLAB
0引言
双音多频（DTMF:
DoubleToneMulti-Frequency）是按键电话通信，也广泛应用于电子邮件和银行系统中。
用户可从电话发送DTMF信号来选择菜单进行操作。
DTMF信号容易用软件
产生和解码。
MATLAB是一个高度集成的软件系统，通过交互式的命令（语句）可以十分简便地实现许多复杂的数值计算。
本文采用MATLAB乍为仿真工具产生DTMF信号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码，由此得出：
时域和频域是研究信号的两个窗口，其中信号处理大都在时域中进行，而信号分析往往在频域中比较方便直观。
且数字信号处理技术中的DFTFFT、卷积、滤波、差分方程这几个概念之间有内在联系。
1DTMF信号的产生
DTMF是数字音频信号，在DTMFI信系统中共有8个频率，分为4个高频音和4个低频音，用一个高频音和一个低频音的组合表示一个信号，这样共有16种组合，分别代表16
种信号，如表1所示：
=〈h（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
表1DTMF信号组合表
例如，当按下数字键“1”时，则产生低频697Hz和高频1209Hz这两个正弦信号的迭加。
由于语音信号的最高频率为4KHz,根据奈奎斯特取样定理，取样频率fs应大于或等于原信号最高频率fc的两倍，即
fs>2fc
（1）
才能保证取样后的信号不失真，所以电话音频信号在数字信号处理时，取样频率fs为
2X4k=8kHz,这里，每个数字信号持续时间为100ms,后面加上100ms的间隔时间（用0表示）。
上述DTMF信号产生方法如下：
（1）建立拨号数字的表矩阵，用查表法（查表1）求用户所按数字键对应的高、低频音。
为简化起见，仅允许选择“0-9”这十个键，在开始时还可拨空信号。
（2）产生相应的DTMFft号及间隔时间。
由于fs=8kHz,各信号持续时间为100ms,因此在程序中每个信号取800点，间隔时间也取800点，结果存入数组中。
（3）画图并监听产生的DTMF信号。
程序如下：
%程序1],DTMF信号的产生
clear%清除内存
TAB=[9411336;6971209;6971336;6971477;7701209;7701336;7701477;852
1209;8521336;8521477];%拨号数字表矩阵
n=input（
l=input（
fori=1:
n
'n='
'n0='
）;%DTMF信号的个数
）;%
空信号点数
k=input（fL=TAB（k+1,1）;%kfH=TAB（k+1,2）;%kn1=800;fs=8000;%j=0:
1:
n1-1;
中。
'0~9'
）;%
输入的数字键
对应的低频音
对应的高频音
产生k对应的DTMF信号，取样频率8kHz,每个信号共1600点，其中，前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;
end
out=out./2;%out数组中每个数据除以2
subplot（211）;plot（out）;%绘图并监听DTMF信号。
sound（out,fs）
wavwrite（out,fs,'D2.wav'）;%out另存入声音文件
程序结果：
例如：
若要产生数字键“1”所对应的DTMF信号，且开始时有200点空信号，则程序运行如
下：
询问：
输入：
n=1
n0=200
0-91
_图1DTMF信号图
1
图1中，前200点为空信号，第201-1000点为键“为间隔时间。
得图1：
1”对应的DTMFft号，第1001-1800点
-0.5
另外，在程序中为使软件设计更接近于实际硬件的开发应用，可用求解差分方程的方法来代
替正弦函数的调用。
设正弦序列为h（n）=sin（w小）u（n）,为实时实现h（n），必须找到其满足
的差分方程。
正弦序列h（n）的z变换为：
H（z）=—zsinWk,令H（z）分子/分母系数：
sinWk=b,2coswk=az—2coswkz+1
则H（z）嵋
zb
~2
z-az1
bz」
1-az二z之
Y（z）（1
121
-az_z）二bzX（z）,
两边进行反变换，得y（n）-ay（n-1）+y（n-2）=bx（n-1）,式中，若令x（n）=S（n）,贝U得到h（n）
的差分方程为h（n）-ah（n-1）+h（n-2）=bS（n-1）,即：
h（n）=ah（n_1）_h（n_2）+bS（n-1）
（2）
用迭代法解此差分方程，即得数字频率为Wk的正弦序列h（n）。
在硬件中，该差分方程是由
加法器、乘法器和单位延时单元构成的系统。
将x（n）=S（n）输入该系统后，输出的就是h（n）。
本文中，每个DTMF言号h（n）是两个频率的正弦序列相迭加，设为hL（n）和环（n）,为此，分
别求得hL（n）和hH（n）所满足的差分方程：
hL（n）=aLhL（n-1）-hL（n-2）+bls（n-1）
;hH（n）=aHhH（n-1）-hH（n-2）+bhS（n-1）
则h（n）=hl（n）+hH（n）。
相关程序如下：
%程序2］，用解差分方程方法产生
x=zeros（1,800）;x
（2）=1;%x（n）=
n=input（'n='）;%DTMF
DTMF言号。
S（n-1）（x（n）
信号的个数
取800点）
l=input（'n0='）;%
out=zeros（1,1600*n+l）;%out
空信号的点数
数组初始化
w=2*pi/8000*[9411336;6971209;6971336;6971477;7701209;7701336;770
1477;8521209;8521336;8521477];%各信号对应的数字频率
tab=[2*cos（w）sin（w）];%各信号满足的差分方程hL（n）和hH（n）中系数a,b的矩阵。
fori=1:
n
k=input（'0-9-'）;%输入的数字键
hL=zeros（1,3）;hH=zeros（1,3）;%hL（n）和hH（n）初始化
forj=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;
end
end
2DTMF信号的解码
解码就是对接收到的DTM信号进行频谱分析，从中找出代表各信号的的特征字，由此
获知用户按下的数字键。
为在频谱图中分辨出不同的频率分量，于是对信号取200点为一帧，
则频谱分辨率F=fs/N=8000/200=40Hz<73Hz（表1中任意两频率的最小间隔），这样，即可满足频谱分析的要求。
本文采用数字信号处理技术中的FFT算法、DFT算法、卷积法及迭代法这4种算法实现对
DTMI信号的解码。
一、快速傅里叶变换（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=2^（E为整数），
所以N=2=256,频谱分辨率F=fs/N-31.25Hz。
因为信号x为实数序列，所其幅频谱|y|具有偶对称性，于是，幅频谱可以仅画N/2点，其中第N/2点对应实际频率为fs/2=4KHz,而DTMF
信号中最高频率为1633Hz,小于2KHz（fs/4）,因此，这里只画N/4=64点。
DTMF!
号是两个正弦波的迭加，它的幅频谱就是两根谱线，谱线的横坐标就是该信号的两个频率分量点K.
和Kh。
2）消除频谱泄漏现象。
由于信号x是有限长的，这就相当于对无限长的信号加矩形窗，所以在频谱图中必然会出现频谱泄漏现象，使信号能量散布到其他谱线位置。
为此，应选择一适当阀值，将出现在这两条谱线周围的幅度较小的谱线消除（置0）,从而解决了这一问题。
最后，将处理后的幅频谱
数据存入数组c中。
3）将各DTMF!
号还原为相应的数字键。
在幅频谱图中，频率轴的定标方式为频率点K而不是实际频率f,转换关系为：
K=f/F,因此，数字键0-9对应频率点如下表所示：
H
40
44
48
23
1
2
3
26
4
5
6
28
7
8
9
31
0
表2数字键对应频率点组合表
数组c的不等于0的下标就是各信号的频率点，查表2,即可将各DTM信号还原为相应的数字
键。
相关程序如下：
%程序3],FFT算法解码程序。
A=wavread（'D2.wav'）;%subplot（212）;plot（A）;%N=256;
fors=1:
8*n%R=out（200*（s-1）+1:
200*s）;y=fft（R,N）;
c（s,:
）=abs（y（1:
64））;%r（s,:
）=c（s,:
）;%r=cz=find（c（s,:
）<40）;%c（s,z）=zeros（size（z））;
%
fori3=1:
8*n
数字键0-9对应的频率点表矩阵sm
b=nnz（c（i3,:
））;%
ifb==2%
q1=find（c（i3,:
））;
fori4=1:
10%ifq1==sm（i4,:
）AN（i3）=i4-1;break;
命令nnz---确定数组中工0数据的个数；
若b=2,则c为信号幅频谱，其工0的下标q1即为频率点。
查表矩阵sm，将q1还原成相应的数字键，存入AN
end
end
else
AN（i3）=NaN;%
end
end
AN=AN%
程序结果：
1）解码：
询问输入
若b丰2，则c为间隔时间，贝UAN=NaN（空信号标志）
显示解码结果AN
n=1（一个DTMFB号）
n0=200（信号前有200点空信号）
0-91（输入数字键“1”）
显示：
AN=NaN1111NaNNaNNaN（8帧信号）
上一行中，第一个“NaN1表示接收到一个空帧（200点），后4帧均为“1”，是解码得到的结果，与按下的数字键相符，且表示信号持续时间为200X4=800点，最后3帧为间隔时间，因为第一帧采到了空信号，所以这里只有三帧。
2）图形
A接收到的DTMFB号同图1
B输入plot（r（2,:
））,得图2
接收到的DTMF信号
绘图
对每帧信号作N=256点的FFT
幅频谱取64点，存入c
消除频谱泄漏现象（阀值=40），结果再存入cendsm=[3144;2340;2344;2348;2640;2644;2648;2840;2844;2848];
图2为数字键““1”对应DTMF言号的幅频谱，可见，它有两条谱线，在频率点但周围发生了频谱泄漏现象。
C输入plot（c（2,:
））,得消除频谱泄漏后的幅频谱图（图3）
20010203040506070
23、40处。
10
图3消除频谱泄漏后的幅频谱图
二、有限长序列离散傅里叶变换（DFT算法
用FFT算法解码每帧信号共涉及256个频率分量，故每帧信号要算N=256点FFT，但实际上，组成所有DTMF信号只用到8个频率分量（见表1的fL和J）,于是，可直接利用
DFT定义式进行频谱分析，且每帧信号只算8点DFT，以避开FFT中许多无意义的计算，且同样达到解码的目的。
DFT算法解码过程如下：
（1）接收DTMF信号，并画图。
（2）对信号作DFT（每帧信号作8点），画频谱图，从中找出代表各信号的特征字。
长度为N的序列x（n），其DFT仍是一个长度为N的序列X（k）,它们的关系是：
N4
X（k）-DFT[x（n）]=6x（n）k=0,1,...N-1（3）
n=0
其中WN=N称为旋转因子，则
k-j^k-fj
Wn=eN=es二e"1叫=cos3k-jsing（4）
这里，fk频率点k对应的实际频率（Hz）
fs取样频率（Hz）
3kk对应的数字频率（rad）
证：
幕
，即
2兀fk
~~3k
Ts
-（4）式成立。
与它在其它频率上DFT的幅值相比，是最大的。
于是，以200点为一帧，对x（n）在8
个特定频率（见表1的fL和fQ上作DFT,并画幅频谱图，从中找出幅值最大的两条谱线，在解决频谱泄漏现象之后，这两条谱线的横坐标就是代表各信号的特征字。
（3）查表3,将各DTMF信号还原为相应的数字键。
H
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=[6977708529411209133614771633];
a1=2*pi/8000;w=a1*w;%DTMF信号用到的八个数字频率存入w数组
for11=1:
8*n%每200点为一帧，在8个特定频率上作DFT，幅频谱数据存入r
fork1=1:
8
s1=0;m=0:
199;
d=cos（w（k1）*m）-j*sin（w（k1）*m）;
a2=out（（l1-1）*200+m+1）;
s仁sum（d.*a2）;
r（l1,k1）=abs（s1）;
end
c（l1,：
）=r（l1,:
）;
z=find（c（l1,:
）<40）;%消除频谱泄漏现象，数据存入c
c（l1,z）=zeros（size（z））;
end
sm=[46;15;16;17;25;26;27;35;36;37];
%sm---数字键0-9对应特征字表矩阵
fori3=1:
8*n%查sm,将特征字还原为相应的数字键，过程与FFT法相似。
b=nnz（c（i3,:
））;
ifb==2
q1=find（c（i3,:
））;
fori4=1:
10
ifq1==sm（i4,:
）
AN（i3）=i4-1;break;
end
end
else
AN（i3）=NaN;
end
end
程序结果：
0.5
1）解码-0.5
输入和同FFT算法所得结杲。
2）图形
-1
0
200400
600
80010001200140016001800
5040302010
图4DTMF信号作8点DFT所得幅频谱图
图4为数字键“1”对应DTM信号的8点DFT的幅频谱。
两条谱线的横坐标“1、5”即为数字键
“1”的特征字。
程序结果:
0.5
0
3）解码
输入和同
4）图形
-0.5
FFT算法所得结杲。
-1
020040060080010001200140016001800
50
40302010
图4DTMF信号作8点DFT所得幅频谱图
图4为数字键“1”对应DTM信号的8点DFT的幅频谱。
两条谱线的横坐标“1、5”即为数字键“1”的特征字。
三、卷积法
当用DFT算法求x（n）的频谱X（k）时，我们发现这也是一卷积过程，下面就具体说明怎样用卷积法解码。
（1）接收DTMF言号，并绘图。
（2）用卷积法检测DTMF言号，找出代表各信号的特征字。
（一）卷积法解码的框图如下图所示：
—竺thk（n）=WN^u（n）】>yk（n）=x（n）*hk（n）—忖tX（k）
图5卷积法解码框图1
图5的方框是一个离散时间系统，广义来说，它就是一个数字滤波器，其单位脉冲响应为:
hk（n）二Wgu（n）,长度为N的因杲序列x（n）通过该系统后，输出为yk（n），而x（n）在某一
频率点k上的DFT:
X（k）=yk（N）（5）
证：
由于线性时不变系统对任意输入信号的响应
综上所述，X（k）就是x（n）与单位脉冲响应为hk（n）=WNknu（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）
先构造一滤波器，其单位脉冲响应为：
』WNNsn=0,1...N-1
h（n八0其它
则卷积法框图变为：
x（n）>h（n厂也（n）.—y（n）：
y（n）>X（k）=y（N-1）：
y（N-1）（N■-1）
图6卷积法框图2
N丄
X（k）Jx（m）W『=x（0）W（x
（1）W,...-x（N—1）W『N」）
m-0-
NJ
又y（n）=x（n）*h（n）-'x（m）h（n-m）
mz0
NJ
.y（N-1）x（m）h（N-1_m）=x（0）h（N-1）x
（1）h（N-2）…x（N-1）h（0）
m-0
=x（0）W^NJXN；）]k■x
（1）Wr；N1XN2）]k■...-x（N-1）wNN1）k=
x（0）WN，x
（1）W,...x（N一1）W,（7,.X（k）=y（N-1）
当N.1时，h（n）=wNNs）k，wNN』）k=W,kWN』k=1.W/k=WN』k=hk（n）
.此时，x（n）通过h（n）相当于通过hk（n）,输出y（n）：
、yk（n）,则
X（k）=y（N-1）：
yk（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---DTMF信号用到的数字频率
forl1=1:
8*n%每200点为一帧，在8个特定频率w
（1）~w（8）上作卷积y=x1*hk
fork1=1:
8%x1---信号；hk---滤波器（hk=cos（wk1m）+jsin（wk1m））
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）满足的差分方程。
1）由hk（n）—Z变换一；Hk（z）——-hk（n）的差分方程
由卷积法得hk（n）二Wn~au（n）,则其z变换为：
Hk（z）'Wn~az"'（Wn~zJ）n
n卫n卫
1Yk（z）
1-Wn'z1_X（z）
，（1-Wn上z丄）
Yk（z）=X（z）,即Yk（z）—Wfz」Yk（z）=X（z）,
对上式两边进行z反变换，得hk（n）满足的差分方程为：
yk（n）一WN上yk（n—1）=x（n）
即:
-k
yk（n）=WNyk（n-1）+x（n）（7）
（7）式对应的系统结构图如下:
x（n）
yk（n）
图8hk（n）系统结构图1
O-
-k
w
则yk（n）|n=N=X（k）o
2）Hk（z）分母有理化。
（7）式中，由于Wk是复数，因此在计算yk（n）时，会碰到复数运算，不方便，为避免复数运算，首先可对H（z）进行分母有理化，其基本原理为：
用一对复共轭极点代替单极点滤波器的思想。
则：
Hk（z）
1
1-Wnz
1-WNzJ
k1k1
（1_Wnz）（1_WnZ）
kk_1_2
1-（
升级会员 