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

双音多频（DTMF）通信设计的MATLAB仿真

摘要：讨论以MATLAB作为仿真工具产生DTMF信号，并用FFT算法、DFT算法、卷积法及迭代法来对DTMF信号进行解码。

关键词：FFT；DFT；频谱分析；卷积；滤波；差分方程；MATLAB

0 引 言

双音多频(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) 1209 1336 1477 1633 fL(Hz) 697 1 2 3 A 770 852 941 4 7 * 5 8 0 6 9 # B C 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:n

k=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; 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-9 1 得图1：

10.50-0.5-1020040060080010001200140016001800 图1 DTMF信号图 1图1中，前200点为空信号，第201-1000点为键“1”对应的DTMF信号，第1001-1800点

0.5为间隔时间。 0

-0.5

-1 020040060080010001200140016001800

另外，在程序中为使软件设计更接近于实际硬件的开发应用，可用求解差分方程的方法来代替正弦函数的调用。设正弦序列为h(n)=sin(ωkn)u(n),为实时实现h(n)，必须找到其满足的差分方程。正弦序列h(n)的z变换为：

H(z)?zsinωk ,令H(z)分子/分母系数：sinωk?b, 2cosωk?a2z?2cosωkz?1Y(z)zbbz?1则H(z)??2?,Y(z)(1?az?1?z2)?bz?1X(z), ?1?2X(z)z?az?11?az?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; end end

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点画频谱图。

E

在MATLAB中,FFT可由语句“y=fft(x,N)”来实现。而FFT中，要求序列长度N=2(E为整数)，

8

所以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信号是两个正弦波的迭加，它的幅频谱就是两根谱线，谱线的横坐标就是该信号的两个频率分量点KL和KH 。

2） 消除频谱泄漏现象。

由于信号x是有限长的，这就相当于对无限长的信号加矩形窗，所以在频谱图中必然会出现 频谱泄漏现象,使信号能量散布到其他谱线位置。为此，应选择一适当阀值，将出现在这两条谱线周围的幅度较小的谱线消除（置0），从而解决了这一问题。最后，将处理后的幅频谱数据存入数组c中。

3） 将各DTMF信号还原为相应的数字键。

在幅频谱图中，频率轴的定标方式为频率点K而不是实际频率f,转换关系为： K=f/F,因此，数字键0-9对应频率点如下表所示： KH KL 40 44 48 23 26 28 31 1 4 7 2 5 8 0 3 6 9 表 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)); end

sm=[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; %若b≠2，则c为间隔时间，则AN=NaN(空信号标志) end end

AN=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信号同图1 B 输入 plot(r(2,:)),得图2

504030201050040010203040506070 图2 DTMF信号作FFT所得的幅频谱图 3050图2为数字键“1”对应DTMF信号的幅频谱，可见，它有两条谱线，在频率点23、40处。 2040但周围发生了频谱泄漏现象。 1030C 输入 plot(c(2,:)),得消除频谱泄漏后的幅频谱图（图3） 020010203040506070

10500403020100010203040506070010203040506070 图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),它们的关系是：

knX(k)?DFT[x(n)]??x(n)WN k?0,1,...N-1 （3）n?0?j2?NN?1其中WN?ekWN?e?j2πkN称为旋转因子，则?j2ω kfs?e?e?jωk?cosωk?jsinωk (4)这里，fk?????频率点k对应的实际频率(Hz) fs?????取样频率(Hz) ωk?????k对应的数字频率（rad)证：?fk?kF??(4)式成立。

fkfsk2? k2? fk??k,即??ωkNNfsNfs(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 2 3 4 1 4 7 2 5 8 0 3 6 9 表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); %消除频谱泄漏现象，数据存入c c(l1,z)=zeros(size(z)); end

sm=[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; end end 程序结果：

10.501)解码 -0.5输入和同FFT算法所得结果。 2)图形 -102004006005040302010001238001000120014001600180045678 图4 DTMF信号作8点DFT所得幅频谱图 图4为数字键“1”对应DTMF信号的8点DFT的幅频谱。两条谱线的横坐标“1、5”即为数字键“1”的特征字。

程序结果：

10.503)解码 -0.5输入和同FFT算法所得结果。 4)图形 -102004006005040302010001238001000120014001600180045678 图4 DTMF信号作8点DFT所得幅频谱图 图4为数字键“1”对应DTMF信号的8点DFT的幅频谱。两条谱线的横坐标“1、5”即为数字键“1”的特征字。

三、卷积法

当用DFT算法求x(n)的频谱X(k)时，我们发现这也是一卷积过程，下面就具体说明怎样用卷积法解码。

(1)接收DTMF信号，并绘图。

(2)用卷积法检测DTMF信号，找出代表各信号的特征字。 （一） 卷积法解码的框图如下图所示：

?N?kn?x(n)???hk(n)?WNu(n)???yk(n)?x(n)*hk(n)?n???X(k)?? 图 5 卷积法解码框图1

图5的方框是一个离散时间系统，广义来说，它就是一个数字滤波器，其单位脉冲响应为：

?knhk(n)?WNu(n),长度为N的因果序列x(n)通过该系统后，输出为yk(n)，而x(n)在某一

频率点k上的DFT:

X(k)=yk(N) (5) 证：

由于线性时不变系统对任意输入信号的响应

-kn

综上所述，X(k)就是x(n)与单位脉冲响应为hk(n)=WNu(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) 先构造一滤波器，其单位脉冲响应为：

(N?1?n)k?WN n?0,1...N-1h(n)?? 其它?0 则卷积法框图变为：?N?1?x(n)????h(n)?hk(n)????y(n)?yk(n)?n???X(k)?y(N?1)?yk(N?1) (N??1) 图 6 卷积法框图2

km0kk(N?1)?X(k)??x(m)WN?x(0)WN?x(1)WN?...?x(N?1)WNm?0N?1又?y(n)?x(n)*h(n)??x(m)h(n?m)m?0N?1?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[N?1?(N?1)]k[N?1?(N?2)]k(N?1)k?x(0)WN?x(1)WN?...?x(N?1)WN?0kk(N?1)x(0)WN?x(1)WN?...?x(N?1)WN,?X(k)?y(N?1)(N?1?n)k(N?n)k?n k?n k?n k?当N??1时，h(n)?WN?WN?WNN k?WN?1?WN?WN?hk(n)N?1?此时，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=[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时刻的输出

2

yk(n)就是X(k),再画|X(k)|,即可对信号进行频谱分析，从而解码。迭代法解码过程如下： (1) 接收DTMF信号，并绘图。

(2) 用迭代法对DTMF信号进行检测，找出代表各信号的特征字。 （一） 求hk(n)满足的差分方程。

变换1)由hk(n)?Z???Hk(z)???hk(n)的差分方程??由卷积法得hk(n)?W??knNu(n),则其z变换为：Hk(z)??Wn?0?knNz?n?k?1n??(WNz)n?0Yk(z)1?k?1?k?1? ，?（1?WNz）Yk(z)?X(z), 即Yk(z)?WNzYk(z)?X(z),?k?1X(z)1?WNz?k对上式两边进行z反变换，得hk(n)满足的差分方程为：yk(n)?WNyk(n?1)?x(n)即： yk(n)=WNyk(n-1)+x(n) (7)

(7)式对应的系统结构图如下：

x(n) yk(n) z-1 -k

图 8 hk(n)系统结构图1 WN

-k

则yk(n)|n=N =X(k)。 2)Hk(z)分母有理化。

-k

(7)式中，由于WN是复数，因此在计算yk(n)时，会碰到复数运算，不方便，为避免复数运算，首先可对Hk(z)进行分母有理化，其基本原理为：用一对复共轭极点代替单极点滤波器的思想。

k?1k?11?WNz1?WNz1则：Hk(z)????k?1?k?1k?1?kk1?WNz(1?WNz)(1?WNz)1?(WN?WN)z?1?z?2? W?kN?W?ekNj2πkN?e-j2πkN2π?2cos( k)?2cosωk ,?Nk?11?WNzHk(z)? (8)1?2cosωk z-1?z?2

Hk(z)分母有理化前后在z平面上的零极点分布如下图所示：

d1=WN

-k

Hk(z)分母有理化

d1=WN

-k

d2=WN=c2

k

图9 Hk(z)零极点分布图

-k

从图9中可得：Hk(z)原来只有一个极点d1=WN（在单位圆上），经分母有理化后，Hk(z)的极

-kkk

点变为2个，即d1=WN、 d2=WN,且又出现一个零点c2=WN,因为d2和c2相互抵消，所以Hk(z) 的表达式不变，且Hk(z)分母中没有出现复数。由（8）式得：

k?11?WNzYk(z)k?1Hk(z)? ?, ?（1?2cosωk z-1?z?2）Yk(z)?（1?WNz）X(z)-1?2X(z)1?2cosωk z?zk对上式两边进行z反变换，得yk(n)?2cosωkyk(n?1)?yk(n?2)?x(n)?WNx(n?1)k即：yk(n)?x(n)?WNx(n?1)?2cosωkyk(n?1)?yk(n?2) (9)此时，hk(n)系统结构图如下：

yk(n) z-1 z -1WNK 2cosωk z -1 图10 hk(n)系统结构图2 式（9）为二阶差分方程，而图10中却有三个单位延时环节，所以图10不是标准形式的二阶系统结构图，可通过正则化，将它变换为只有两个延时环节的标准形式（称为正准型）。 正则化步骤如下：

(8)式中，令H1(z)?k?1H2(z)?1?WNz?1V(z)? (8')?1?21?2cosωkz?zX(z)Yk(z) (8''),则 Hk(z)?H1(z)?H2(z) V(z)则H1(z)满足的差分方程为v(n)?x(n)?2cosωkv(n?1)?v(n?2)　　(9')kH2(z)满足的差分方程为yk(n)?v(n)?WNv(n?1)　　(9'')由卷积定理：频域相乘对应于时域相卷得：hk(n)?h1(n)*h2(n) (10)因此可将hk(n)拆成h1(n)和h2(n)相串联的形式，再合并相同的单位延时环节，得下图：

x(n) v(n) yk(n)

2cosωk WNK 图11 hk(n)系统结构图3（正准型）

z

z-1 -1 显然，图11是正准型的二阶系统结构图。综上所述，hk(n)满足的差分方程为 (9’)和（9’’）式，即：

(9')??v(n)?x(n)?2cosωkv(n?1)?v(n?2)　　??ky(n)?v(n)?Wv(n?1)　　 (9'')N?k?k且X(k)?yk(N)?v(N)?WNv(N?1) (11)（二） 画|X(k)|,对信号进行频谱分析。

(9’)式中，n=0,1…N-1,共进行N次实数乘法，而(9’’)中，仅在n=N时刻计算1次乘法，但 K2

WN还是复数。由于解码算法中，不关心X(k),只要求出|X(k)|或|X(k)|即可进行频谱分析，

2

这里，算|X(k)|。 由式（11）得，

kk|X(k)|2?|v(N?WNv(N?1) |2 (其中WN?cos?k?jsin?k)2

?|v(N)?cos?kv(N?1)?jsin?kv(N?1)|2?[v(N)?cos?kv(N?1)]2?sin2?kv2(N?1)?v2(N)?2cos?kv(N)v(N?1)?v2(N?1) (12)这样解码过程中就完全避免了复数运算。 余下解码过程同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 0点为一帧，在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);

22 %|x(k)|2= v(N)?2cos?kv(N)v(N?1)?v(N?1)

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次 256点 在MATLAB中 复数乘法（N=256 （|X(k)|） 极易，但算法k取256点） （最大） 不“透明” （算|X(k)|） （最少） DFT法 N×8=1600次 8点 较难 复数乘法（N=200 （|X(k)|） k取8点） （算|X(k)|） 卷积法 同DFT法 同DFT法 与DFT法 相当 迭代法 (N+4)×8=1632次同DFT法 最难，但更贴2实数乘法 （|X(k)|） 近硬件实现 2（算|X(k)|）

4 结论

本文采用四种算法对仿真的DTMF信号解码，它们的共同点在于对信号作DFT,从而进行

频谱分析。DFT的快速算法为FFT,而DFT又是一种运算，可通过将信号输入一数字滤波器 （离散时间系统），取其在某点上的输出实现，上述方法既是一卷积过程，又是一解差分方程的过程。由此得出：DFT、FFT、卷积、滤波、差分方程这几个概念之间有内在联系。

MATLAB的应用，大大减少了计算工作量，且使比较抽象的概念得以直观显示，它不愧为功能强大的计算机辅助教学工具。

参考文献 [1] 戴悟僧，冯筱林，安邦健.数字信号处理导论.[M].上海：科学技术出版社，2000. [2]吴湘淇.信号、系统与信号处理（下）.[M].北京：电子工业出版社，2000.

本文来源：https://www.bwwdw.com/article/r94v.html