CWT-Based时频分析
这个例子展示了如何使用连续小波变换(CWT)对信号进行时间和频率的联合分析。这个例子讨论了瞬态的局部化，其中CWT优于短时傅里叶变换(STFT)。算例还说明了如何利用逆小波变换合成时频局部化的信号近似。
在一些例子中对CWT和STFT进行了比较。您必须有信号处理工具箱™来运行光谱图的例子。
调制信号的时频分析
加载一个二次啁啾信号并显示它的声谱图。在t= 0时，信号的频率大约从500hz开始，在t=2时降低到100hz，在t=4时又增加到500hz。采样频率为1khz。
负载quadchirp；fs = 1000;[S、F T] =光谱图(quadchirp, 100、98128 fs);helperCWTTimeFreqPlot (S T F,“冲浪”，二次啁啾的短时傅里叶变换，“秒”，“赫兹”）
使用CWT获得信号的时频图。
(慢性疲劳综合症,f) = cwt (quadchirp fs,“WaveletParameters”[14200]);helperCWTTimeFreqPlot (cfs tquad f,“冲浪”，“二次啁啾的CWT”，“秒”，“赫兹”）
带有凹凸小波的连续小波变换产生的时频分析与短时傅立叶变换非常相似。
频率和幅度调制经常发生在自然信号中。使用CWT获得一个大棕色蝙蝠(Eptesicus Fuscus)发出的回声定位脉冲的时频分析。采样间隔为7微秒。使用bump小波，每八度有32个声音。感谢伊利诺伊大学贝克曼中心的Curtis Condon、Ken White和Al Feng提供的蝙蝠数据，以及在本例中使用蝙蝠的许可。
负载batsignalt = 0: DT(元素个数(batsignal) * DT) DT;[cfs f] = cwt (batsignal“撞”1 / DT,“VoicesPerOctave”、32);helperCWTTimeFreqPlot (cfs t * 1 e6, f / 1000“冲浪”，“蝙蝠回声定位的CWT”，.．.微秒的，“赫兹”）
获取并绘制蝙蝠数据的STFT。
F (S, T) =光谱图(batsignal、50、48128、1 / DT);helperCWTTimeFreqPlot (S, t . * 1 e6, f / 1 e3,“冲浪”，“蝙蝠回声定位(STFT)”，.．.微秒的，“赫兹”）
对于模拟信号和自然调制信号，CWT提供了类似于STFT的结果。
利用连续小波变换检测振荡中的瞬态
在时频分析的某些情况下，连续小波变换比短时间傅里叶变换提供更多信息的时频变换。当信号被瞬变损坏时，就会发生这样的情况。这些瞬变现象的出现和消失通常具有物理意义。因此，除了表征信号中的振荡成分外，重要的是能够本地化这些瞬态。为了模拟这一点，创建一个由两个频率为150和200赫兹的正弦波组成的信号。采样频率为1khz。正弦波有不相交的时间支撑点。金宝app150赫兹的正弦波出现在100到300毫秒之间。200赫兹的正弦波从700毫秒到1秒。另外，在222毫秒和800毫秒有两个瞬态。 The signal is corrupted by noise.
rng默认的；dt = 0.001;t = 0: dt: 1 dt;addNoise = 0.025 * randn(大小(t));x = cos(2 *π* 150 * t) * (t > = 0.1 & t < 0.3) +罪(2 *π* 200 * t) * (t > 0.7);x = x + addNoise;X ([222 800]) = X ([222 800]) +[-2 2];图;情节(t。* 1000 x);包含(的毫秒）;ylabel (“振幅”）;
放大这两个瞬态，可以看到它们代表150和200赫兹振荡中的扰动。
次要情节(2,1,1)情节(t(184:264)。* 1000 x (184:264));ylabel (“振幅”)标题(瞬态的)轴紧；次要情节(2,1,2)情节(t(760:840)。* 1000 x (760:840));ylabel (“振幅”)轴紧；包含(的毫秒）
利用解析Morlet小波获得并绘制CWT。
图;[cfs f] = cwt (x, 1 / dt,“爱”）;轮廓(t * 1000 f, abs (cfs))网格在c = colorbar;包含(的毫秒) ylabel (“频率”) c.Label.String =“级”；
分析Morlet小波比bump小波具有较差的频率局部化，但具有较好的时间局部化。这使得Morlet小波成为瞬态局部化的更好选择。绘制幅值平方精细尺度系数来演示瞬态的局部化。
图;wt = cwt (x, 1 / dt,“爱”）;情节(t。* 1000、abs (wt (1:)) ^ 2);包含(的毫秒）;ylabel (“级”）;网格在；标题(解析Morlet小波-精细尺度系数）;持有从
小波收缩使瞬变的时间定位具有很高的精度，而拉伸使振荡在150和200赫兹的频率定位成为可能。
短时傅里叶变换只能将瞬态局部化到窗口的宽度。你必须以分贝(dB)为单位绘制短时傅立叶变换图，以便能够可视化瞬变。
[S、F T] =光谱图(x, 50、48128、1000);冲浪(t . * 1000 F, 20 * log10 (abs (S)),“edgecolor”，“没有”）;视图(0,90);轴紧；阴影插值函数；colorbar包含(“时间”), ylabel (“赫兹”）;标题(“STFT”）
瞬变现象在短时傅立叶变换中只出现在功率的宽带增加。比较STFT在第一次瞬态出现之前(居中于183毫秒)和之后(居中于223毫秒)获得的短时功率估价值。
图;情节(20 * log10 (abs (S (:, 80))));持有在；情节(20 * log10 (abs (S (:, 100))),“r”）;传奇('T = 183 msec'，'T = 223毫秒')包含(“赫兹”）;ylabel (“数据库”）;持有从；
短时傅里叶变换不能提供与CWT相同程度的瞬态定位能力。
用逆CWT去除时间局域频率分量
创建一个由指数加权正弦波组成的信号。有两个25-Hz的分量——一个以0.2秒为中心，另一个以0.5秒为中心。有两个70赫兹的组件——一个以0.2秒为中心，另一个以0.8秒为中心。前25-Hz和70-Hz分量在时间上同时发生。
t = 0:1/2000:1-1/2000;dt = 1/2000;x1 =罪(50 *π* t) * exp(-50 *π* (t - 0.2) ^ 2);x2 =罪(50 *π* t) * exp(-100 *π* (t - 0.5) ^ 2);x3 = 2 * cos(140 *π* t) * exp(-50 *π* (t - 0.2) ^ 2);x4 = 2 * sin(140 *π* t) * exp(-80 *π* (t - 0.8) ^ 2);x = x1 + x2 + x3 + x4;图;情节(t, x)标题(叠加信号的）
获取并显示CWT。
类(x, 2000);标题(“使用默认莫尔斯小波的解析CWT”）;
通过将CWT系数归零去除大约0.07到0.3秒的25hz分量。使用逆CWT (icwt)来重建信号的近似值。
(慢性疲劳综合症,f) = cwt (x, 2000);T1 = 07;T2 = .33;F1 = 19;F2 = 34;cfs(f > F1 & f < F2, t> T1 & t < T2) = 0;xrec = icwt (cfs);
显示重构信号的CWT。初始的25hz组件被移除。
类(xrec, 2000)
绘制原始信号和重建信号。
次要情节(2,1,1);情节(t, x);标题(原始信号的）;次要情节(2,1,2);情节(t, xrec)标题(“去掉前25赫兹分量的信号”）;
最后，将重构信号与没有以0.2秒为中心的25-Hz分量的原始信号进行比较。
y = x2 + x3 + x4;图;xrec情节(t)在情节(t y“r——”)传说(“逆CWT近似”，“没有25hz的原始信号”）;持有从
通过解析CWT确定精确频率
当你使用解析小波获得正弦波的小波变换时，解析CWT系数实际上编码了频率。
为了说明这一点，考虑从人耳获得的耳声发射。耳声发射(OAEs)是由耳蜗(内耳)发出的，它的存在表明听力正常。加载和绘制OAE数据。数据以20千赫采样。
负载dpoae情节(t。* 1000,dpoaets)包含(的毫秒) ylabel (“振幅”）
这种发射是由一个从25毫秒开始到175毫秒结束的刺激引起的。根据实验参数，发射频率为1230 Hz。获取并绘制CWT。
类(dpoaets“撞”, 20000);包含(的毫秒）;
您可以通过找到频率最接近1230 Hz的CWT系数并检查它们的幅度作为时间的函数来研究OAE的时间演化。用时间标记标记唤起刺激的开始和结束。
[dpoaeCWT f] = cwt (dpoaets“撞”, 20000);[~, idx1230] = min (abs (f - 1230));cfsOAE = dpoaeCWT (idx1230:);情节(t。* 1000、abs (cfsOAE));持有在甘氨胆酸AX =;情节(25 [25],[AX.YLim (1) AX.YLim (2)),“r”) plot([175 175]，[AX.YLim(1) AX.YLim(2)]，“r”)包含(的毫秒), ylabel (“级”）;标题(rocky系数大小的）
在唤醒刺激开始和耳声发射之间存在一定的延迟。一旦唤起刺激终止，耳声发射立即开始衰减。
另一种分离发射的方法是使用反CWT在时域重建频率局部逼近。
在[11501350]Hz的频率范围内，通过反转CWT重构频率局域发射近似。将原始数据连同重建和指示唤起刺激开始和结束的标记绘制出来。
xrec = icwt(dpoaeCWT,f，[1150 1350]);图绘制(t。* 1000,dpoaets)在；xrec情节(t。* 1000年,“r”) AX = gca;ylim = AX.YLim;ylim情节(25 [25],“k”)情节(ylim (175 175),“k”)包含(的毫秒) ylabel (“振幅”)标题(“发射的频率局部重建”）
在时域数据中，您可以清楚地看到在触发刺激的应用和终止时，发射是如何启动和关闭的。
重要的是要注意，即使为重建选择了一个频率范围，解析小波变换实际上是对发射的确切频率进行编码。为了证明这一点，对解析CWT重建的发射近似进行傅里叶变换。
xdft = fft (xrec);频率= 0:2e4 /元素个数(xrec): 1 e4;xdft = xdft(1:元素个数(xrec) / 2 + 1);图绘制(频率、abs (xdft));包含(“赫兹”）;ylabel (“级”)标题(基于cwt的信号近似的傅里叶变换）;[~, maxidx] = max (abs (xdft));流('频率为%4.2f Hz\n'频率(maxidx));
频率为1230.00 Hz
结论与进一步阅读
在这个例子中，您学习了如何使用CWT来使用解析小波对一维信号进行时频分析类．您看到了一些信号的例子，其中CWT提供了与短时傅里叶变换相似的结果，还有一个例子，CWT可以提供比短时傅里叶变换更可解释的结果。最后，您学习了如何使用以下方法重构信号的时间尺度(频率)局部逼近icwt．
有关CWT的更多信息，请参见小波工具箱文档。
参考:Mallat, S。《信号处理的小波之旅:稀疏方法》，学术出版社，2009年。
附录
本例中使用了以下辅助函数。
你也可以从以下列表中选择一个网站:
