双音多频通信设计的Matlab仿真.docx
- 文档编号:24307788
- 上传时间:2023-05-26
- 格式:DOCX
- 页数:26
- 大小:110.59KB
双音多频通信设计的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-(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 双音 通信 设计 Matlab 仿真