一、实验目的
1. 了解数字信号处理当今应用的基本情况。
2. 对该课程做系统地总结。
3．将所学知识运用到实践中，能够学以致用。
二、实验原理
双音多频 DTMF（ Dual Tone Multi-Frequency ）信号，是用两个特定的单音频率信号的组合来代表数字或功能。在 DTMF 电话机中有 16 个按键，其中 10 个数字键 0—9 ， 6 个功能键 * 、 # 、 A 、 B 、 C 、 D 。其中 12 个按键是我们比较熟悉的按键另外由第 4 列确定的按键作为保留，作为功能键留为今后他用。 根据 CCITT建议，国际上采用 697Hz 、770Hz 、852Hz 、94lHz 低频群及 1209Hz 、1336Hz 、1477Hz 、1633Hz高频群。对应的 DTMF 信号用sin(2πf₁t) + sin(2πf₂t) 表示。从低频群和高频群任意各抽出一种频率进行组合，共有 16 种组合，代表 16 种不同的数字键或功能，每个按键唯一地由一组行频和列频组成，如表1 所示。[10]
（1）DTMF信号的产生合成
一般用数字方法产生DTMF信号：规定用8kHz对DTMF信号进行采样，采样后得到时域离散信号为：x(n) = sin(2πf1n / 8000) + sin(2πf2n / 8000)（式4-1）
有计算法和查表法两种方法可形成上述序列，计算法运行速度慢，查表法速度快但要占用一定的存储空间。因为采样频率是8kHz，既要求每125ms输出一个样本，得到的序列再送到D/A转换器。输出经过平滑滤波便是连续的DTMF信号， DTMF信号通过电话线路再送到交换机。
（2）DTMF信号的检测
要实现电话拨号音（ DTMF ）信号的检测识别，可以先将收到的连续信号经过 A/D 变换，变成数字信号再进行检测。可用DFT(FFT)对双音多频信号进行频谱分析，由信号的幅度谱判断信号的两个频率，最后确定对应的数字和符号。当检测的频率数目较少时，用滤波器组实现更适合。FFT 是 DFT的快速算法，但当要计算的频率点数目远小于 DFT 的区间长度时，用 FFT 快速算法的效果并不明显，而且还要占用很多内存，因此不如直接用 DFT 合适。戈泽尔（Goertzel）算法是计算 DFT的一种线性滤波方法，用该算法实现软件实验。对于戈泽尔算法建议同学们自己查阅相关资料，掌握算法的基本思想。
（3） 检测 DTMF 信号的DFT 参数选择
用 DFT 检测模拟信号 DTMF 信号所含的两个音频频率，即为用 DFT 对模拟信号进行频谱分析的问题。根据数字信号处理所学知识，要确定三个参数：采样频率 Fs ，DFT 的变换点数 N，需要对信号的观察时间的长度Tp 。这三个参数不能随意选取，要根据对信号频谱分析要求确定。对信号频谱分析有三个要求：
1．频谱分析的分辨率。观察要检测到的 8 个频率，相邻间隔最小的是第一和第二个频率，间隔是 73Hz，即要求 DFT 至少能分辨相隔 73Hz 的两个信号。所以对信号的观察时间应为Tp = 1/ F = 1/ 73 = 13.7ms 。考虑到可靠性，应留有富裕量，要求按键时间在 40ms 以上。
2. 频谱分析的频率范围 要检测信号的频率范围是 697—1633Hz，但考虑到存在语音干扰，除检测这 8 个频率外，还要检测他们的二次倍频的幅度大小。这样频谱分析的频率范围是697-3266Hz。按照采样定理，则有最小采样频率为Fs min = 6.35kHz ，系统规定 Fs = 8kHz ，这样可算出对信号的最少采样点数为Nmin = Tp min Fs ≈ 110 。
3. 检测频率的准确性 这是一个用 DFT 检测正弦波频率是否准确的问题。序列的 DFT 是序列频域函数在0 − 2π 区间的等间隔采样，如果是周期序列，截取周期序列的整个周期进行 DFT，其采样点刚好在周期信号的频率上，DFT 的幅度最大处准确的是信号频率。分析 DTMF信号不可能经过采样得到周期序列，因此存在一个检测频率的准确性问题。
DFT 的频率采样点频率为wk = 2πk / N (k = 0,1,2,..., N − 1) ，相应的在模拟域的采样点频率为 fk = Fsk / N (k = 0,1,2,..., N − 1) ，要选择一个合适的 N 值使得用该公式计算出的 fk 能接近要检测的实际频率。或者将实际的频率值带入能得到的 k值接近整数值。经过分析研究，发现 N=205 是最好的。随意将实际频率带入上面的公式可得到最接近整数的 K 值。例如：
通过一系列的分析，可确定： Fs = 8kHz, N = 205,Tp ≥ 40ms
工具箱函数 Goerztel 的调用格式为：Xgk=goertzel(xn,K+1)
其中 xn 是被变换的时域序列，用于 DTMF 信号检测时就是 DTMF 信号的 205 个采样值。K 是要求计算的 DFT 频点的序号向量。用 N 表示 xn 的长度，则要求 0≤ K ≤ N − 1。由以上的分析可知只计算 DTMF 信号 8 个基频时K=[18,20,22,24,31,34,38,42]
Xgk 是变换结果向量，存放的是由 K 指定的频率点的 DFT[x(n)]的值[10]。
三、 实验内容
要求输入 6 位或 8 位电话号码，利用戈泽尔算法求出各数字信号的频谱并绘出频谱图形，再检测输出电话号码。
参考流程图：
四、实验报告要求
1.简述实验目的及原理。明白具体的问题解决的思路是什么。
2.整理好经过运行并证明是正确的实验程序并加上注释，绘出相应的图形。
f1=697;f2=770;f3=852;f4=941;F1=1209;F2=1336;F3=1477;F4=1633;N=205;tmin=40/1000;fs=8000;dt=1/fs;N1=floor(tmin/dt);if N<=N1t=[0:N-1]*dt;elseerror('The time is not enough!');endk1=sin(2*pi*f1*t)+sin(2*pi*F1*t);k2=sin(2*pi*f1*t)+sin(2*pi*F2*t);k3=sin(2*pi*f1*t)+sin(2*pi*F3*t);ka=sin(2*pi*f1*t)+sin(2*pi*F4*t);k4=sin(2*pi*f2*t)+sin(2*pi*F1*t);k5=sin(2*pi*f2*t)+sin(2*pi*F2*t);k6=sin(2*pi*f2*t)+sin(2*pi*F3*t);kb=sin(2*pi*f2*t)+sin(2*pi*F4*t);k7=sin(2*pi*f3*t)+sin(2*pi*F1*t);k8=sin(2*pi*f3*t)+sin(2*pi*F2*t);k9=sin(2*pi*f3*t)+sin(2*pi*F3*t);kc=sin(2*pi*f3*t)+sin(2*pi*F4*t);km=sin(2*pi*f4*t)+sin(2*pi*F1*t);k0=sin(2*pi*f4*t)+sin(2*pi*F2*t);kj=sin(2*pi*f4*t)+sin(2*pi*F3*t);kd=sin(2*pi*f4*t)+sin(2*pi*F4*t);key=['1','2','3','a';'4','5','6','b';'7','8','9','c';'*','0','#','d'];k=[18,20,22,24,31,34,38,42];num=input('please enter the key:','s');num=num-48;nn=length(num);disp('The number of the key is: ');disp(nn);number=zeros(nn,length(t));for i=1:nnswitch num(i)case 1number(i,1:N)=k1;case 2number(i,1:N)=k2;case 3number(i,1:N)=k3;case 4number(i,1:N)=k4;case 5number(i,1:N)=k5;case 6number(i,1:N)=k6;case 7number(i,1:N)=k7;case 8number(i,1:N)=k8;case 9number(i,1:N)=k9;case 0number(i,1:N)=k0;case 49number(i,1:N)=ka;case 50number(i,1:N)=kb;case 51number(i,1:N)=kc;case 52number(i,1:N)=kd;case -6number(i,1:N)=km;case -13number(i,1:N)=kj;otherwiseerror('The key is not right!');endenddisp('The key is: ');for i=1:nnxgk=goertzel(number(i,1:N),k+1);figure;x=[697,770,852,941,1209,1336,1477,1633];stem(x,abs(xgk),'.');xlabel('Hz');title('频谱分析');zb=find(abs(xgk)>50);disp(key(zb(1),zb(2)-4));end
结果: 随机输入数字:12378
