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

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

考虑数字电话拨号的信号系统。这种系统产生的信号被称为双音多频(DTMF)信号。每个拨出的数字所产生的声音由两个正弦波的和组成 $-$或音调 $-$频率取自两个相互排斥的群体。每一对音调包含低组的一个频率(697 Hz, 770 Hz, 852 Hz，或941 Hz)和高组的一个频率(1209 Hz, 1336 Hz，或1477Hz)，代表一个独特的符号。以下是拨给电话簿按键的频率:

生成DTMF信号并收听它。

[tone, Fs] = helperDTMFToneGenerator();p = audioplayer(音调、Fs、16);玩(p)

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

N =元素个数(音调);t = (0: n - 1) / Fs;次要情节(2,1,1)阴谋(1 e3 * t,音调)包含(“时间(ms)”) ylabel (“振幅”)标题(DTMF信号的次要情节(2,1,2)pspectrum(音调,Fs,“漏”, 1“FrequencyLimits”(650、1500))

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

env =信封(音调,80,“rms”）;脉冲宽度(env, Fs)

ans =3×10.1041 0.1042 0.1047

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

这里你可以看到三个脉冲，每个大约100毫秒长。但是，您无法分辨拨了哪些号码。一个频域图可以帮助你弄清楚，因为它显示了信号中存在的频率。

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

f = [meanfreq(tone,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));

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

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

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

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

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

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

t =(0:长度(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”10 e - 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千赫，以获得关于每个脊的频率变化的更多细节。

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

注意，现在频率脊更好地定位在频率中。但是，由于频率分辨率和时间分辨率是成反比的，所以谱图的时间分辨率要小得多。设置99%的重叠以平滑时间窗口。使用一个“MinThreshold”的 $-$60db以删除不需要的背景内容。

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

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

时频再赋值

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

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

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

现在频率脊更清晰，在时间上也更精确。您还可以使用该函数本地化信号能量fsst，这将在下一节中讨论。

重建时频脊

考虑下面的录音，包括一个频率随时间下降的啁啾信号和一个最后的啪啪声。

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

让我们通过提取时频平面上的山脊来重建“splat”声音的一部分。我们使用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．沿着变换线画出山脊。

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

接下来，利用脊指数向量重建chirp信号iridge．在山脊两边各放一个箱子。绘制重构信号的谱图。

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

重建山脊消除了信号中的噪声。连续播放带噪信号和去噪信号，听其区别。

p = audioplayer ([yChirp; 0(大小(yChirp)); yrec], Fs, 16);玩(p);

测量能力

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

Fs = 1 e8;bw = 60 e6;t = 0:1 / Fs: 10 e-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 kHz采样的对数啁啾。啁啾的频率在10秒内从10hz增加到400hz。

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

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

甘氨胆酸ax =;斧子。YScale =“日志”；

三维瀑布可视化

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

Fs = 10 e3;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));

该信号还包含210hz的干扰，振幅为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．您学习了如何配置声谱图以获得对数频率刻度和三维可视化。最后，您学习了如何通过计算持久谱来找到干扰信号。

附录

本例中使用了以下辅助函数。