使用时频分析确定的数字在DTMF信号

你可以把几乎所有的时变信号分成足够短的时间间隔，使得信号在每个部分都是基本平稳的。时频分析最常用的方法是将信号分割成短周期，然后通过滑动窗口估计频谱。该pspectrum功能与使用“频谱”option计算每个滑动窗口的基于fft的光谱估计，并让您可视化信号的频率内容如何随时间变化。

考虑数字电话拨号的信号系统。这种系统产生的信号称为双音多频(DTMF)信号。每个拨号所产生的声音由两个正弦信号之和组成 $-$或音调 $-$与来自两个相互排斥的基团取频率。每对音调的含有低组（697赫兹，770赫兹，852赫兹，或941赫兹）和高组（1209赫兹，1336赫兹，或1477Hz）中的一个频率的一个频率，表示唯一的符号。下面是分配给一个电话垫的按钮频率：

产生DTMF信号和收听。

[音调，FS] = helperDTMFToneGenerator（）;P = audioplayer（色调，FS，16）;玩(p)

听着信号，你可以告诉大家，一个三位数字被拨打。但是，你不能告诉它是哪个号码。接着，在时间上和在650至1500赫兹频带频率域可视化的信号。设置'泄漏'的参数pspectrum函数以1来使用矩形窗口并提高频率分辨率。

N = numel（音调）;T =（0：N-1）/ FS;副区（2,1,1）图（1E3 * T，音）xlabel（'时间（ms）')ylabel ('振幅'）标题（DTMF信号的）副区（2,1,2）pspectrum（色调，FS，'泄漏'1，'FrequencyLimits'[650，1500]）

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

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

ans =3×10.1041 0.1042 0.1047

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

在这里你可以看到三个脉冲，每一个大约100毫秒长。但是，你不能告诉被拨哪些号码。频域图可以帮助你摸不着头脑，因为这表明存在于信号的频率。

通过在四个不同的频带推定平均频率找到的频率峰值。

F = [meanfreq（色调，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分贝功率电平以显现只有主频率分量。要一次查看使用0％重叠音的持续时间和它们的位置。

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

频谱编码频率功率电平的颜色。黄色的颜色表示频率内容以更高的功率;蓝色表示频率含量非常低的功率。一个强大的黄水平线表示一个音调的在特定频率的存在。该图清楚地显示了1336 Hz的音调在所有三个拨号数字的存在，告诉你，他们都是在键盘的第二列。从图中可以看出，最低频率，770赫兹，首先拨打。最高频率941赫兹，是未来。中间频率，852赫兹，排在最后。因此，所拨打的号码是508。

权衡时间和频率分辨率，以获取你的信号的最佳表示

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

请考虑以下信号，以4khz采样，由太平洋蓝鲸歌声的颤音部分组成:

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

该颤音信号由音调脉冲串的。看看所获得的时间信号和频谱pspectrum当没有指定分辨率和时间分辨率设置为10毫秒。设置'泄漏'参数为1以使用矩形窗口。既然我们要本地化的脉冲的时间位置，设置重叠％至0。最后，使用“MinThreshold”的 $-$60分贝从谱图视图去除背景噪声。

T =（0：长度（whaleTrill）-1）/ FS;图AX1 =副区（3,1,1）;情节（T，whaleTrill）AX2 =副区（3,1,2）;pspectrum（whaleTrill，FS，“频谱”，'OverlapPercent'0,...'泄漏'1，“MinThreshold”，-60）彩条（AX2，“关”) ax3 = subplot(3,1,3);pspectrum（whaleTrill，FS，“频谱”，'OverlapPercent'0,...'泄漏'1，“MinThreshold”,-60,“TimeResolution”，10E-3）彩条（AX3，“关”）linkaxes（[AX1，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”，3E3）

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

pspectrum (batsignal Fs,“频谱”，“FrequencyResolution”，3E3，...'OverlapPercent'99，“MinTHreshold”,-60)

新设置得到谱图清楚地示出了回声定位信号的四个频率脊。

时频再赋值

尽管我们已经能够确定四个频率的山脊，我们仍然可以看到每个脊分布在几个相邻频率点。这是由于在时间和频率两者中使用的窗口化方法的泄漏。

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

pspectrum (batsignal Fs,“频谱”，“FrequencyResolution”，3E3，...'OverlapPercent'99，“MinTHreshold”,-60,“再分配”，真正）

现在频率脊更加清晰和及时更好的本地化。您也可以使用功能定位信号能量fsst，这是在下一节中讨论。

重建时频岭

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

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

让我们通过提取在时频平面脊重建“图示”声的一部分。我们用fsst为了锐化带噪版本的splat信号的频谱，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。沿着变换绘制山脊。

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

接下来，使用脊索引向量重构所述线性调频脉冲信号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 = 1E8;BW = 60e6;t = 0时：1 / FS：10E-6;ICOMP =啁啾（T，-BW / 2，T（结束），（体重）/ 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分贝。

对数频标可视化

在某些应用中，在对数频率标度上可视化信号的光谱图可能是更好的选择。你可以通过改变YScale财产的ÿ-轴。例如，考虑在1千赫取样的对数线性调频脉冲。啁啾的频率从10Hz在10秒内增加到400赫兹。

FS = 1E3;t = 0时：1 / FS：10;FO = 10;F1 = 400;Y =啁啾（吨，FO，10中，f1，“对数”）;pspectrum（Y，FS，“频谱”，“FrequencyResolution”1，...'OverlapPercent'90，'泄漏',0.85,'FrequencyLimits',1 f / 2)

啁啾的频谱变成直线时在频率标度是对数的。

AX = GCA;ax.YScale =“日志”;

三维可视化瀑布

随着视图命令，您可以将信号的光谱图可视化为三维瀑布图。还可以使用。更改显示颜色色彩表函数。

FS = 10E3;t = 0时：1 / FS：2;X1 = VCO（锯齿（2 * PI *吨，0.5），[0.1 0.4] * Fs的，FS）;pspectrum（X1，FS，“频谱”，'泄漏',0.8)

视图（-45,65）的颜色表骨

使用持久性光谱寻找干扰

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

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

FS = 1000;T =（0：1 / FS：500）';X =啁啾（吨，180，T（端部），220）+ 0.15 * randn（大小（T））;

该信号还包含一个210赫兹的干扰，以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'，[100 290]）

计算信号的持久性光谱。现在，这两个信号分量都清晰可见。

图colormapparulapspectrum（X，FS，“坚持不懈”，'FrequencyLimits'，[100 290]，“TimeResolution”，1）

结论

在这个例子中，你学会了如何使用pspectrum函数执行时频分析，以及如何解释光谱数据和功率水平。你学会了如何更改时间和频率分辨率以改善信号的理解，以及如何使用锐化光谱和提取时频脊fsst，ifsst,tfridge。你学会了如何配置频谱情节获得对数频率刻度和三维可视化。最后，您了解了如何通过计算一个持久性的频谱找到干扰信号。

附录

以下的辅助函数被用于本实施例。