时频分析实践介绍

您所在的位置:网站首页 matlab如何对数据进行频率分析 时频分析实践介绍

时频分析实践介绍

2023-07-21 04:19| 来源: 网络整理| 查看: 265

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

您可以将几乎任何时变信号划分为足够短的时间区间,这样,信号在每个区间内基本上是平稳的。通常,时频分析是通过将一个信号分割为若干短周期并在滑动窗内估计频谱来进行的。与 'spectrogram' 选项结合使用的 pspectrum 函数计算每个滑动窗内基于 FFT 的频谱估计值,让您直观地看到信号的频率成分如何随时间变化。

以某数字电话拨号的信号系统为例。这种系统产生的信号称为双音多频 (DTMF) 信号。每个拨号号码生成的声音由两个正弦波 - 即音调 - 的总和组成,频率取自两个互斥的组。每对音调包含一个低组频率(697 Hz、770 Hz、852 Hz 或 941 Hz)和一个高组频率(1209 Hz、1336 Hz 或 1477 Hz),并表示一个唯一的符号。以下是分配给电话键盘上各按键的频率:

生成一个 DTMF 信号并聆听该信号。

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

聆听该信号,可分辨出拨了一个三位号码。然而,您无法分辨出拨打的具体号码。接下来,在时域和 650 至 1500 Hz 频带的频域上可视化信号。将 pspectrum 函数的 'Leakage' 参数设置为 1,以使用矩形窗并提高频率分辨率。

N = numel(tones); t = (0:N-1)/Fs; subplot(2,1,1) plot(1e3*t,tones) xlabel('Time (ms)') ylabel('Amplitude') title('DTMF Signal') subplot(2,1,2) pspectrum(tones,Fs,'Leakage',1,'FrequencyLimits',[650, 1500])

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

env = envelope(tones,80,'rms'); pulsewidth(env,Fs)ans = 3×1 0.1041 0.1042 0.1047 title('Pulse Width of RMS Envelope')

此处您可以看到三个脉冲,每个脉冲长度大约为 100 毫秒。然而,您无法分辨出拨打的是哪个数字。频域图有助于您解决此问题,因为它显示信号中存在的各频率。

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

f = [meanfreq(tones,Fs,[700 800]), ... meanfreq(tones,Fs,[800 900]), ... meanfreq(tones,Fs,[900 1000]), ... meanfreq(tones,Fs,[1300 1400])]; round(f)ans = 1×4 770 852 941 1336

通过将估计的频率与电话键盘的图相匹配,您可以识别出拨打的按键是“5”、“8”和“0”。然而,频域图没有提供任何类型的时间信息来让您知道这些数字的拨打顺序。这些数字的组合可以是“580”、“508”、“805”、“850”、“085”或“058”。为了解决此难题,请使用 pspectrum 函数计算频谱图,并观察信号的频率成分如何随时间变化。

计算 650 至 1500 Hz 频带的频谱图,并去除低于 - 10 dB 功率水平的内容,以仅可视化主频率分量。要查看音调持续时间及其在时间上的位置,请使用 0% 重叠。

pspectrum(tones,Fs,'spectrogram','Leakage',1,'OverlapPercent',0, ... 'MinThreshold',-10,'FrequencyLimits',[650, 1500]);

频谱图的颜色对频率功率水平进行编码。黄色表示功率较高的频率成分;蓝色表示功率非常低的频率成分。一条亮黄色水平线表示存在某特定频率的音调。该图清楚地显示拨打的所有三个数字中都存在一个 1336 Hz 音调,告诉您它们都在键盘的第二列中。从图中可以看出,先拨打的是最低频率 770 Hz。然后拨打的是最高频率 941 Hz。最后拨打的是中间频率 852 Hz。因此,拨打的号码是 508。

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

pspectrum 函数将一个信号划分为若干个段。较长的段提供更好的频率分辨率;较短的段提供更好的时间分辨率。可以使用 'FrequencyResolution' 和 'TimeResolution' 参数控制段的长度。当未指定频率分辨率或时间分辨率值时,pspectrum 会尝试根据输入信号长度在时间分辨率和频率分辨率之间找到良好的平衡。

假设有以下信号,其采样率为 4 kHz,包含太平洋蓝鲸鲸吟的颤音部分:

load whaleTrill p = audioplayer(whaleTrill,Fs,16); play(p)

颤音信号由一系列音调脉冲组成。查看在未指定分辨率和时间分辨率设置为 10 毫秒时的时间信号和由 pspectrum 获得的频谱图。将 'Leakage' 参数设置为 1 以使用矩形窗。由于我们要定位脉冲的时间位置,因此将重叠百分比设置为 0。最后,使用 -60 dB 的 'MinThreshold' 去除频谱图视图中的背景噪声。

t = (0:length(whaleTrill)-1)/Fs; figure ax1 = subplot(3,1,1); plot(t,whaleTrill) ax2 = subplot(3,1,2); pspectrum(whaleTrill,Fs,'spectrogram','OverlapPercent',0, ... 'Leakage',1,'MinThreshold',-60) colorbar(ax2,'off') ax3 = subplot(3,1,3); pspectrum(whaleTrill,Fs,'spectrogram','OverlapPercent',0, ... 'Leakage',1,'MinThreshold',-60,'TimeResolution', 10e-3) colorbar(ax3,'off') linkaxes([ax1,ax2,ax3],'x')

pspectrum 选择的 47 毫秒时间分辨率不够小,不足以定位频谱图中的所有颤音脉冲。另一方面,10 毫秒的时间分辨率足以在时间上定位每个颤音脉冲。如果我们放大几个脉冲,这会变得更加清晰:

xlim([0.3 0.55])

现在加载一个信号,其中包含由一只大棕色蝙蝠发出的回声定位脉冲。信号以 7 微秒的采样间隔测量。分析信号的频谱图。

load batsignal Fs = 1/DT; figure pspectrum(batsignal,Fs,'spectrogram')

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

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3)

注意现在频率脊能够更好地在频率上定位。然而,由于频率分辨率和时间分辨率成反比,频谱图的时间分辨率明显更小。设置 99% 的重叠来对时间窗进行平滑处理。使用 -60 dB 的 'MinThreshold' 来去除不需要的背景内容。

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3, ... 'OverlapPercent',99,'MinTHreshold',-60)

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

时频重排

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

pspectrum 函数能够在时间和频率上估计每个频谱估计值的能量中心。如果您将每个估计值的能量重新赋予最接近新时间和频率中心的 bin,您可以校正该窗的一些泄漏。您可以通过使用 'Reassign' 参数来完成此操作。将此参数设置为 true 会计算重排后的信号频谱图。

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3, ... 'OverlapPercent',99,'MinTHreshold',-60,'Reassign',true)

现在频率脊变得更尖锐,在时间上可更好地定位。您也可以使用 fsst 函数来定位信号能量,这将在下一节中讨论。

重新构造时频脊线

假设有以下音频录制,其中包含频率随时间降低的 chirp 信号,且最后有一个啪嗒声。

load splat p = audioplayer(y,Fs,16); play(p) pspectrum(y,Fs,'spectrogram')

让我们通过在时频平面中提取脊来重新构造“啪嗒”声的一部分。我们使用 fsst 来锐化含噪版本啪嗒信号的频谱,使用 tfridge 来识别 chirp 声的脊,并使用 ifsst 来重新构造 chirp。该过程对重新构造的信号进行去噪。

将高斯噪声添加到“啪嗒”声的 chirp 部分。添加的噪声对用廉价麦克风录制的音频进行仿真。检查时频频谱内容。

rng('default') t = (0:length(y)-1)/Fs; yNoise = y + 0.1*randn(size(y)); yChirp = yNoise(t


【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3