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

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

考虑一下数字电话拨号的信号系统。这种系统产生的信号被称为双音多频(DTMF)信号。每个拨出的号码所发出的声音是两个正弦波的和 $-$或音调 $-$频率来自两个互斥的组。每对音调包含一个低频率组(697 Hz, 770 Hz, 852 Hz，或941 Hz)和一个高频率组(1209 Hz, 1336 Hz，或1477Hz)，代表一个唯一的符号。以下是分配给电话板按钮的频率:

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

[tones, Fs] = helperDTMFToneGenerator();p = audioplayer(音调，Fs,16);玩(p)

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

N =数字(音调);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

标题(“均方根包络的脉冲宽度”）

这里你可以看到三个脉冲，每个脉冲大约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赫兹波段，并删除以下的内容 $-$10db功率电平仅显示主频组件。要查看音调持续时间及其位置，请使用0%重叠。

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

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

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

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

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

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

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

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”， 10e-3) colorbar(ax3，“关闭”) linkaxes (ax₁,ax2 ax3],“x”）

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

xlim ([0.3 - 0.55])

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

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

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

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来重建啁啾。该过程对重构信号进行去噪处理。

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

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);[fridge, 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)

用于计算频谱图的参数给出了线性调频信号的清晰时频表示。pspectrum计算功率谱图，这意味着颜色值对应于真实功率水平，单位为dB。颜色条显示信号的功率水平是左右 $-$4 dB。

对数频率尺度可视化

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

Fs = 1e3;t = 0:1/Fs:10;Fo = 10;F1 = 400;Y = chirp(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 Hz增加到220 Hz。

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

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

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

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

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

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

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

结论

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

附录

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