\(\begin{align*} \newcommand{\dif}{\mathop{}\!\mathrm{d}} \newcommand{\belowarrow}[1]{\mathop{#1}\limits_{\uparrow}} \newcommand{\bd}{\boldsymbol} \newcommand{\L}{\mathscr{L}} \end{align*}\)
实验目的
- 了解电话按键音形成的原理,理解DTMF音频产生软件和DTMF解码算法;
- 利用FFT算法识别按键音。
实验要求
- 设计音频产生函数,音频信号见下图,每个数据信号持续半秒;
- 实现解码函数:接收(1)产生的DTMF信号,识别信号的频率,并生成包含拨号数字的序列。
实验过程
什么是DTMF?DTMF(Dual Tone Multi Frequency),双音多频,由高频群和低频群组成,高低频群各包含4个频率。一个高频信号和一个低频信号叠加组成一个组合信号,代表一个数字。高频、低频及对应的数字如下图所示:
频率 | 1209Hz | 1336Hz | 1477Hz | 1633Hz |
---|---|---|---|---|
697Hz | 1 | 2 | 3 | A |
770Hz | 4 | 5 | 6 | B |
852Hz | 7 | 8 | 9 | C |
941Hz | * | 0 | # | D |
CCITT规定每秒最多按10个键,即每个键时隙最短为100MS,其中音频实际持续时间至少为45MS,不大于55MS,时隙的其他时间内保持静默,因此按键产生双音频信号时,相继的两个信号间隔一段时间;解码器利用这个时间识别出双音频信号,并转换成对应的数字信息,而且要识别出间隙信息。因此流程包含音频任务和静默任务,前者是产生双音频采样值,后者产生静默样值,每个任务结束时,要重置定时器和下一个任务。其中静默任务还要加上一个任务:从数字缓冲区取出数字并解包。
function compound =dtmfDial(numString, Fs)
% DTMF Dial: generate a sound according to the numString
% numString: the telephone numbers
% Fs: sampling rate
fl=[697 770 852 941];%低频频率
fh=[1209 1336 1477 1633];%高频频率
last_time=0.5;%单个按键声音持续时间
compound=[];
for i=1:length(numString)
switch numString(i)
case'1'
freq_low=fl(1);freq_hgh=fh(1);
case'2'
freq_low=fl(1);freq_hgh=fh(2);
case'3'
freq_low=fl(1);freq_hgh=fh(3);
case'4'
freq_low=fl(2);freq_hgh=fh(1);
case'5'
freq_low=fl(2);freq_hgh=fh(2);
case'6'
freq_low=fl(2);freq_hgh=fh(3);
case'7'
freq_low=fl(3);freq_hgh=fh(1);
case'8'
freq_low=fl(3);freq_hgh=fh(2);
case'9'
freq_low=fl(3);freq_hgh=fh(3);
case'*'
freq_low=fl(4);freq_hgh=fh(1);
case'0'
freq_low=fl(4);freq_hgh=fh(2);
case'#'
freq_low=fl(4);freq_hgh=fh(3);
case'A'
freq_low=fl(1);freq_hgh=fh(4);
case'a'
freq_low=fl(1);freq_hgh=fh(4);
case'B'
freq_low=fl(2);freq_hgh=fh(4);
case'b'
freq_low=fl(2);freq_hgh=fh(4);
case'C'
freq_low=fl(3);freq_hgh=fh(4);
case'c'
freq_low=fl(3);freq_hgh=fh(4);
case'D'
freq_low=fl(4);freq_hgh=fh(4);
case'd'
freq_low=fl(4);freq_hgh=fh(4);
otherwise
error('The n!');
end
single=0.25*sin(2*pi*freq_low*[1/Fs:1/Fs:last_time])+... %低频
0.25*sin(2*pi*freq_hgh*[1/Fs:1/Fs:last_time]); %高频
single=[single,zeros(1,0.3*Fs)]; %加上间隔
compound=[compound,single];%将每个按键串在一起
end
end
numString='123456789*0#';
Fs=8000;
compound = dtmfDial(numString, Fs);
plot(compound);
sound(compound,Fs);%播放声音
%audiowrite('test.wav',compound,Fs);%保存声音文件
产生的音频的时域图如下图所示:
为了检测每个数据,我们需要将音频切片。若音频中一段时间内的平均幅度过小,则可认为是静默。因此,我们需要对整个音频进行滑动平均滤波,即:
\[y(n)=\frac{1}{windowSize}\big[ x(n)+x(n-1)+\cdots+x(n-(windowSize-1)) \big]\]窗口不能取太大,也不能取太小,否则会导致滤波效果不好。由于最小静默时间为 45ms,所以窗口可以取 5ms. 可以用 filter()
函数进行滤波
function [nstart, nstop]=dtmfCut(compound, Fs)
%return the start and end index of each number
%nstarts: a vector of the start index
%nends: a vector of the end index
compound = compound(:)'/max(abs(compound)); %归一化
len = length(compound);
windowSize = round(0.01*Fs);
setpoint = 0.02 %低于 2% 的点视为静默
compound = filter(ones(1,windowSize)/windowSize, 1, abs(compound)); %滑动平均滤波
compound = diff(compound>setpoint); %差分
jmp = find(compound~=0)'; %找出跳变点
if compound(jmp(1))<0 %第1个跳变点为负,说明信号开始时是拨号音
jmp=[1 jmp];
end
if compound(jmp(end))>0 %最后跳变点为正,说明信号结束时是拨号音
jmp=[jmp len];
end
index=[];
while length(jmp>1)
if jmp(2)>(jmp(1)+10*windowSize) %若两个跳变点时间太短则忽略
index = [index; jmp(1:2)];
end
jmp(1:2)=[];
end
nstart = index(1,:);
nstop = index(2,:);
end
经过上面的处理后,可以得到每个片段的起始点与终止点,如下图所示:
然后对每个片段进行 DFT 变换,找出峰值。这里我们以第 1 个片段为例子。
for kk=1:length(nstart)
indx=[nstart(kk):nstop(kk)];
x_seg(kk,[1:length(indx)])=compound(indx);
end
x_fft = abs(fft(x_seg(1,:)));
[pks, locs]=findpeaks(x_fft,'MINPEAKHEIGHT',max(x_fft)/4);
f=[1:length(x_fft)];
plot(f, x_fft, f(locs), pks, 'or')
得到的 DFT 和峰值如下图所示。
找出峰值后,即可根据下面的公式求出频率:
\[\Omega_k = 2 \cdot \frac{k}{N} \cdot f_s\\ 其中,k<\frac{N}{2}\]fk=locs/length(x_fft)*Fs
求出频率为 698.6 和 1210.7。CCITT 规定频率偏移在 1.5% 内为正常,因此可以认为这两个频率近似于 697 和 1209,从而识别出该片段代表 “1”
将上面过程用函数实现,如下。
function numString = dtmfRun(compound, Fs)
fl=[697 770 852 941];%低频频率
fh=[1209 1336 1477 1633];%高频频率
num=['1','2','3','A';
'4','5','6','B';
'7','8','9','C';
'*','0','#','D'];
[nstart, nstop]=dtmfCut(compound, Fs); %获取每个片段的开始与结束位置
for kk=1:length(nstart) %获取每个片段
indx=[nstart(kk):nstop(kk)];
x_seg(kk,[1:length(indx)])=compound(indx);
end
x_freq=[];
for kk=1:length(nstart) %寻峰
x_fft=abs(fft(x_seg(kk,:)));
[~, locs]=findpeaks(x_fft,'MINPEAKHEIGHT',max(x_fft)/4);
fk=locs/length(x_fft)*Fs;
x_freq=[x_freq;fk(1:2)];
end
numString=[];
for kk=1:length(nstart) %频率转数字
lnum=0;hum=0;
f1=x_freq(kk,1);
f2=x_freq(kk,2);
lnum=find(abs(fl-f1)./fl<0.015);
hnum=find(abs(fh-f2)./fh<0.015);
if lnum~=0 && hnum~=0
numString=[numString,num(lnum,hnum)];
end
end
end
最后验证一下代码的正确性:
numString='123456789*0#';
Fs=8000;
compound = dtmfDial(numString, Fs);
dtmfRun(compound, Fs)
输出:
ans =
'123456789*0#'
结果正确。