利用时频分析识别DTMF信号中的数字

你几乎可以把任何时变信号分成足够短的时间间隔，使信号在每个部分基本上是静止的。时频分析最常用的方法是将信号分割成这些短周期并在滑动窗口上估计频谱。的pspectrum函数与的谱图选项在每个滑动窗口上计算基于fft的频谱估计，并让您可视化信号的频率内容如何随时间变化。

以数字电话拨号的信号系统为例。这种系统产生的信号称为双音多频(DTMF)信号。每个拨打的号码所产生的声音由两个正弦信号的和组成 $-$或音调 $-$从两个相互排斥的组中提取频率。每一对音调包含一个低组的频率(697 Hz, 770 Hz, 852 Hz，或941 Hz)和一个高组的频率(1209 Hz, 1336 Hz，或1477Hz)，代表一个唯一的符号。以下是分配给电话簿按键的频率:

生成一个DTMF信号并监听它。

[tones, Fs] = helperDTMFToneGenerator();p = audioplayer(tone,Fs,16);玩(p)

听信号，你可以知道拨的是一个三位数的号码。然而，你不知道是哪个数字。接下来，将信号在650到1500 Hz波段的时域和频域上可视化。设置“漏”参数pspectrum函数到1可以使用矩形窗口并提高频率分辨率。

N = numel(音调);t = (0:N-1)/Fs;Subplot (2,1,1) plot(1e3*t,tones) xlabel(“时间(ms)”) ylabel (“振幅”)标题(DTMF信号的) subplot(2,1,2) pspectrum(tone,Fs，“漏”，1，“FrequencyLimits”(650、1500))

信号的时域图证实了三次能量爆发的存在，对应于三个按下的按钮。要测量爆发的长度，可以取RMS包络的脉冲宽度。

Env =包络(音，80，“rms”）;脉冲宽度(env, Fs)

ans =3×10.1041 0.1042 0.1047

标题(“RMS包络脉宽”）

这里你可以看到三个脉冲，每个脉冲大约100毫秒长。但是，你不知道拨的是哪个号码。频域图可以帮助您解决这个问题，因为它显示了信号中出现的频率。

通过估计四个不同频段的平均频率来定位频率峰值。

f = [meanfreq(tones,Fs，[700 800])，.．.meanfreq(音调,Fs 900 [800]),.．.meanfreq(音调,Fs 1000 [900]),.．.meanfreq(音调,Fs 1400 [1300])];轮(f)

ans =1×4770 852 941 1336

通过将估计的频率与电话簿的图表相匹配，您可以说拨打的按钮是“5”、“8”和“0”。然而，频域图不提供任何类型的时间信息，使您能够找出拨号的顺序。组合可以是'580'、'508'、'805'、'850'、'085'或'058'。为了解决这个难题，使用pspectrum函数计算谱图，观察信号的频率含量如何随时间变化。

计算650到1500赫兹波段的光谱图，并删除以下的内容 $-$10分贝的功率水平，只可视化的主要频率组件。要查看音调持续时间和它们在时间上的位置，请使用0%重叠。

pspectrum(音调,Fs,的谱图，“漏”，1，“OverlapPercent”0,.．.“MinThreshold”, -10,“FrequencyLimits”(650、1500));

光谱图的颜色编码频率功率等级。黄色表示功率较高的频率内容;蓝色表示极低功率的频率内容。一条强烈的黄色水平线表示存在一个特定频率的音调。图中清楚地显示了在所有三个拨号数字中都存在1336赫兹的音，告诉您它们都在键盘的第二列上。从图中你可以看到，最低的频率，770赫兹，是先拨的。其次是最高频率941赫兹。中频852hz排在最后。因此，拨打的号码是508。

权衡时间和频率分辨率，以获得信号的最佳表现

的pspectrum函数将信号分成若干段。较长的片段提供更好的频率分辨率;较短的片段提供更好的时间分辨率。控件可以控制段长度“FrequencyResolution”而且“TimeResolution”参数。当没有指定频率分辨率或时间分辨率值时，pspectrum试图在基于输入信号长度的时间和频率分辨率之间找到一个良好的平衡。

考虑以下4千赫采样的信号，它包含了太平洋蓝鲸歌曲的颤音部分:

负载whaleTrillp = audioplayer(whaleTrill,Fs,16);玩(p)

颤音信号由一串音调脉冲组成。看得到的时间信号和谱图pspectrum当没有指定分辨率且时间分辨率设置为10毫秒时。设置“漏”参数设为1以使用矩形窗口。因为我们想要定位脉冲的时间位置，所以将重叠百分比设置为0。最后，使用“MinThreshold”的 $-$60 dB从光谱图视图中去除背景噪声。

t = (0:length(whaleTrill)-1)/Fs;图ax1 = subplot(3,1,1);plot(t,whaleTrill) ax2 = subplot(3,1,2);pspectrum (whaleTrill Fs,的谱图，“OverlapPercent”0,.．.“漏”，1，“MinThreshold”, -60) colorbar (ax2,“关闭”) ax3 = subplot(3,1,3);pspectrum (whaleTrill Fs,的谱图，“OverlapPercent”0,.．.“漏”，1，“MinThreshold”, -60,“TimeResolution”colorbar(ax3;“关闭”) linkaxes (ax₁,ax2 ax3],“x”）

选择47毫秒的时间分辨率pspectrum不足以定位光谱图中的所有颤音脉冲。另一方面，10毫秒的时间分辨率足以及时定位每个颤音脉冲。如果我们放大几个脉冲，就会更清楚:

xlim ([0.3 - 0.55])

现在加载一个信号，由一只大棕色蝙蝠(Eptesicus fuscus)发出的回声定位脉冲组成。用7微秒的采样间隔测量信号。分析信号的频谱图。

负载batsignalFs = 1/DT;图pspectrum (batsignal Fs,的谱图）

默认参数值的谱图显示四个粗时频脊。将频率分辨率值降低到3khz，以获得每个脊的频率变化的更多细节。

pspectrum (batsignal Fs,的谱图，“FrequencyResolution”, 3 e3)

观察到现在频率脊在频率上被更好地定位了。然而，由于频率和时间分辨率成反比，谱图的时间分辨率是相当小的。设置99%的重叠来平滑时间窗口。使用一个“MinThreshold”的 $-$60分贝删除不需要的背景内容。

pspectrum (batsignal Fs,的谱图，“FrequencyResolution”3 e3,.．.“OverlapPercent”, 99,“MinTHreshold”, -60)

新的设置产生一个光谱图，清楚地显示回声定位信号的四个频率脊。

时频再赋值

即使我们已经能够识别四个频率脊，我们仍然可以看到每个脊分布在几个相邻的频率箱上。这是由于在时间和频率上使用的加窗方法的泄漏。

的pspectrum函数能够在时间和频率上对每个谱估计的能量中心进行估计。如果您将每个估计的能量重新分配到最接近新的时间和频率中心的容器，您可以纠正窗口的一些泄漏。可以通过使用“再分配”参数。将此参数设置为真正的计算信号的重新分配谱图。

pspectrum (batsignal Fs,的谱图，“FrequencyResolution”3 e3,.．.“OverlapPercent”, 99,“MinTHreshold”, -60,“再分配”,真正的)

现在频率脊更清晰，时间定位也更好。您还可以使用函数定位信号能量fsst，这将在下一节中讨论。

重构时频脊

考虑下面的录音，由频率随时间下降的啁啾信号和最后的啪啪声组成。

负载长条木板p = audioplayer(y,Fs,16);玩(p) pspectrum (y, Fs,的谱图）

让我们通过提取时频平面上的脊来重建“啪啦啪啦”声音的一部分。我们使用fsst为了锐化杂波信号的频谱，tfridge为了识别啁啾声的脊，和ifsst重建啁啾。该过程对重构信号进行降噪处理。

将高斯噪声添加到“飞溅”声音的啁啾部分。添加的噪音模拟了用廉价麦克风录制的音频。检查时频谱内容。

rng (“默认”t =(0:长度(y)-1)/Fs;yNoise = y + 0.1*randn(size(y));yChirp = yNoise(t<0.35);pspectrum (yChirp Fs,的谱图，“MinThreshold”, -70)

使用傅里叶同步压缩变换锐化频谱，fsst．fsst通过在固定时间内在频率上重新分配能量，在时频平面上定位能量。计算并绘制噪声啁啾的同步压缩变换。

fsst (yChirp Fs,“桠溪”）

啁啾在时频平面上表现为局部的脊。用以下方法识别山脊tfridge．沿着变换画出山脊。

[sst,f] = fsst(yChirp,Fs);[冰箱，iridge] = tfridge(sst,f,10);helperPlotRidge (yChirp、Fs、冰箱);

接下来，利用脊指数向量重构啁啾信号iridge．在山脊两侧各设一个垃圾箱。绘制重构信号的谱图。

Yrec = ifsst(sst,kaiser(256,10)，iridge，“NumFrequencyBins”1);pspectrum (yrec Fs,的谱图，“MinThreshold”, -70)

重建脊去掉了信号中的噪声。连续播放噪声信号和去噪信号来听出区别。

p = audioplayer([yChirp;zero (size(yChirp));yrec]，Fs,16);玩(p);

测量能力

考虑一个复杂的线性调频(LFM)脉冲，这是一种常见的雷达波形。使用1.27微秒的时间分辨率和90%的重叠计算信号的谱图。

Fs = 1e8;Bw = 60e6;t = 0:1/Fs:10e-6;IComp = chirp(t，-bw/2,t(end)， bw/2，“线性”, 90) + 0.15 * randn(大小(t));QComp = chirp(t，-bw/2,t(end)， bw/2，“线性”, 0) + 0.15 * randn(大小(t));IQData = IComp + 1i*QComp;segmentLength = 128;pspectrum (IQData Fs,的谱图，“TimeResolution”1.27 e-6“OverlapPercent”, 90)

用于计算谱图的参数给出了LFM信号的清晰时频表示。pspectrum计算功率谱图，这意味着颜色值对应真正的功率水平，单位为dB。颜色条表示信号的功率水平在附近 $-$4 dB。

对数频率刻度可视化

在某些应用中，最好是在对数频率标度上可视化信号的谱图。可以通过更改YScale的属性y设在。例如，考虑在1khz采样的对数啁啾。啁啾的频率在10秒内从10hz增加到400hz。

Fs = 1e3;t = 0:1/Fs:10;Fo = 10;F1 = 400;Y =啁啾(t,fo,10,f1，“对数”）;pspectrum (y, Fs,的谱图，“FrequencyResolution”，1，.．.“OverlapPercent”, 90,“漏”, 0.85,“FrequencyLimits”, 1 f / 2)

当频率尺度为对数时，啁啾的谱图变为一条直线。

Ax = gca;斧子。YScale =“日志”；

三维瀑布可视化

与视图命令，您可以将信号的光谱图可视化为三维瀑布图。控件还可以更改显示颜色colormap函数。

Fs = 10e3;t = 0:1/Fs:2;x1 = vco(锯齿(2*pi*t,0.5)，[0.1 0.4]*Fs,Fs);pspectrum (x1, Fs,的谱图，“漏”, 0.8)

colormap视图(-45、65)骨

利用持续谱发现干扰

信号的持续频谱是一种时间-频率视图，显示给定频率在信号中出现的时间百分比。持久谱是工频空间的直方图。在信号演化过程中，一个特定频率在信号中持续的时间越长，它的时间百分比就越高，因此它在显示器中的颜色就越亮或“热”。使用持续频谱来识别隐藏在其他信号中的信号。

考虑一个嵌入宽带信号中的窄带干扰信号。产生一个啁啾采样在1千赫500秒。在测量过程中，啁啾的频率从180赫兹增加到220赫兹。

Fs = 1000;T = (0:1/fs:500)';x =唧唧喳喳(220 t, 180 t(结束),)+ 0.15 * randn(大小(t));

该信号还包含210hz的干扰，振幅为0.05，仅在总信号持续时间的1/6中出现。

Idx =地板(长度(x)/6);(1: idx) = x (1: idx) + 0.05 * cos(2 *π* t (1: idx) * 210);

计算信号在100到290 Hz区间内的功率谱。微弱的正弦信号被啁啾掩盖了。

pspectrum (x, fs,“FrequencyLimits”290年[100])

计算信号的持续频谱。现在两个信号成分都清晰可见了。

图colormapparulapspectrum (x, fs,“坚持不懈”，“FrequencyLimits”(100 290),“TimeResolution”, 1)

结论

在本例中，您学习了如何使用pspectrum函数执行时频分析，以及如何解释谱图数据和功率级别。您学习了如何改变时间和频率分辨率以提高您对信号的理解，以及如何锐化光谱和提取时频脊fsst，ifsst,tfridge．您学习了如何配置谱图图，以获得对数频率刻度和三维可视化。最后，您学习了如何通过计算持久性谱来查找干扰信号。

附录

本例中使用了以下helper函数。