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

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

考虑数字电话拨号的信号系统。这种系统产生的信号被称为双音多频(DTMF)信号。每个拨出的数字所产生的声音由两个正弦波的和组成 $-$或音调 $-$频率从两个相互排斥的群体中取出。每对音调含有低组（697Hz，770 Hz，852 Hz或941Hz）的一个频率，以及高组的一个频率（1209Hz，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]）

信号的时域图确认存在三个能量突发，对应于三个推动的按钮。为了测量突发的长度，您可以采用RMS包络的脉冲宽度。

ent =信封（音调，80，'rms'）;脉冲宽（ENV，FS）

ans =3×10.1041 0.1042 0.1047

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

在这里，您可以看到三个脉冲，每个脉冲大约100毫秒长。但是，您无法拨打拨打的数字。频域绘图可帮助您解决此问题，因为它显示了信号中存在的频率。

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

f = [算法（音调，FS，[700 800]），......meanfreq(音调,Fs 900 [800]),......Meanfreq（音调，FS，[900 1000]），......Meanfreq（音调，FS，[1300 1400]）]];圆形（f）

ans =1×4770 852 941 1336

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

将频谱图计算在650到1500 Hz频段上并删除下面的内容 $-$10 DB功率电平仅可视化主频率分量。及时查看音调持续时间和它们的位置使用0％重叠。

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

谱图的颜色编码频率电量。黄色表示具有更高功率的频率内容;蓝色表示具有非常低的功率的频率内容。强烈的黄色水平线表示特定频率处的音调存在。该曲线清楚地显示了所有三位拨号的数字中的1336 Hz音调，告诉您它们都在键盘的第二列上。从情节中，您可以先看到最低频率，770 Hz首先拨打。最高频率，941 Hz是接下来的。中频，852赫兹持续了。因此，拨号号码为508。

交易时间和频率分辨率以获得信号的最佳表示

这pspectrum.功能将信号划分为段。更长的段提供更好的频率分辨率;较短的段提供更好的时间分辨率。可以使用该段长度使用“FrequencyResolution”和“TimeResolution”参数。当没有指定频率分辨率或时间分辨率值时，pspectrum.尝试根据输入信号的长度在时间和频率分辨率之间找到一个良好的平衡。万博 尤文图斯

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

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

轨道信号由一系列色调脉冲组成。查看时间信号和频谱图pspectrum.未指定分辨率时，将时间分辨率设置为10毫秒时。设定“漏”参数为1使用矩形窗口。由于我们想要本地化脉冲的时间位置，因此将重叠百分率设置为0。最后，使用a“MinThreshold”的 $-$60 dB以从频谱图视图中删除背景噪声。

t =（0：长度（鲸状）-1）/ fs;图AX1 =子图（3,1,1）;图（T，鲸状物）AX2 =子图（3,1,2）;PSPectrum（鲸鱼，FS，'谱图'那'重叠的'0,......“漏”，1，“MinThreshold”，-60）颜色栏（AX2，“关闭”) ax3 = subplot(3,1,3);PSPectrum（鲸鱼，FS，'谱图'那'重叠的'0,......“漏”，1，“MinThreshold”, -60,“TimeResolution”，10e-3）颜色栏（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”，3e3）

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

pspectrum (batsignal Fs,'谱图'那“FrequencyResolution”，3e3，......'重叠的', 99,“MinTHreshold”, -60)

新设置产生了一种频谱图，该频谱图清楚地显示了回声机器的四个频率脊。

时频再赋值

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

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

pspectrum (batsignal Fs,'谱图'那“FrequencyResolution”，3e3，......'重叠的', 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信号ir．在山脊的每一侧包括一个垃圾箱。绘制重建信号的谱图。

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'重叠的'，90）

用来计算谱图的参数给出了线性调频信号清晰的时频表示。pspectrum.计算功率谱图，这意味着颜色值对应于DB中的真正功率电平。彩色条显示信号的功率电平周围 $-$4 dB。

对数频率尺度可视化

在某些应用中，最好是在对数频率尺度上可视化信号的谱图。你可以通过改变YScale财产的y设在。例如，考虑一个在1 kHz采样的对数啁啾。啁啾的频率在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，......'重叠的'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)

查看（-45,65）ColorMap骨

使用持久谱寻找干扰

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

考虑嵌入在宽带信号中的干扰窄带信号。产生一个啁啾采样在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．您学习了如何配置声谱图以获得对数频率刻度和三维可视化。最后，您学习了如何通过计算持久谱来找到干扰信号。

附录

在此示例中使用以下辅助功能。