连续与离散小波分析

人们经常有一个问题:什么是连续关于连续小波分析?毕竟，你可能对在计算机上做小波分析感兴趣，而不是用铅笔和纸，在计算机中，没有什么是真正连续的。当术语连续小波分析在科学计算设置中使用，它是指每八度以上有一个小波的小波分析技术，或频率的两倍，其中小波之间的时间漂移是一个样本。这使得得到的连续小波变换(CWT)具有两个在应用中非常有用的特性:

此外，在CWT的大多数应用中，与实值小波相对的复值小波是有用的。主要原因是复值小波包含相位信息。看到比较两种信号的时变频率内容这是一个相位信息有用的例子。本教程中的例子只使用复值小波。

参见[1]用于详细处理小波信号处理，包括用复值小波进行连续小波分析。

过滤器或每八度的声音

一个常用来表示每八度的小波滤波器的数量的术语是每八度的声音．为了说明这一点，构造一个默认连续小波滤波器组。

Fb = cwtfilterbank

fb = cwtfilterbank with properties: VoicesPerOctave: 10 Wavelet: 'Morse' SamplingFrequency: 1 SamplingPeriod: [] PeriodLimits: [] SignalLength: 1024 FrequencyLimits: [] TimeBandwidth: 60 Wavelet parameters: [] Boundary: 'reflection'

如果检查滤波器组的属性，默认情况下，每个八度(或10个)有10个小波滤波器VoicesPerOctave)．画出小波滤波器的频率响应。

freqz(神奇动物)

小波滤波器在连续分析中占有重要地位常q所有小波滤波器的特性，即它们在频率或带宽上的扩展与它们的中心频率成正比。换句话说，小波滤波器在高频率时比在低频率时更宽。由于时间和频率之间的倒数关系，这意味着高频小波比低频小波在时间上有更好的局部化(有更小的时间支持)。万博1manbetx为了看到这一点，从滤波器组中提取中心频率和时域小波。画出两个小波的实部，一个高频，一个低频，以便比较。

Psi = wavelet (fb);F = centerFrequencies(fb);Figure plot(real(psi(9，:))) ylabel(“振幅”) yyaxis正确的情节(真实(ψ(end-16,:))) ylabel (“振幅”)轴紧S1 =“具有中心频率的小波”；S2 = [num2str(F(9))，' % 1.2 f '）和"num2str (F (end-16),' % 1.2 f '）“周期/样本”];标题({S1;S2})包含(“样本”)传说(“0.25 c/s小波”，“0.01 c/s小波”）

在频率响应图中，中心频率不超过0.45个周期/样本。求滤波器组在一个八度或频率的两倍内有多少个小波滤波器。确认数字等于VoicesPerOctave．

numFilters10 = nnz(F >= 0.2 & F<= 0.4)

numFilters10 = 10

在创建筛选器组时，可以指定不同数量的筛选器。创建一个滤波器组，每个八度有8个声音，并绘制频率响应。

Fb = cwtfilterbank(“VoicesPerOctave”8);图freqz (fb)

确认滤波器组中每个八度有八个滤波器。

F = centerFrequencies(fb);numFilters8 = nnz(F >= 0.2 & F<= 0.4)

numFilters8 = 8

您可以阅读连续和离散小波分析之间的差异的更详细的解释连续和离散小波变换．

对数间隔中心频率

小波分析中人们会感到困惑的一个方面是滤波器的对数间距。

绘制小波滤波器组的中心频率。

subplot(2,1,1) plot(F) ylabel(“周期/样本”)标题(“小波中心频率”网格)在子图(2,1,2)符号学(F)网格在ylabel (“周期/样本”)标题(“半对数的规模”）

图显示，小波中心频率不是线性间隔，通常情况下与其他滤波器组。具体来说，中心频率呈指数递减，因此高中心频率之间的步长大于低中心频率之间的步长。为什么这对小波有意义?回想一下，小波是常数q滤波器，这意味着它们的带宽与它们的中心频率成正比。如果你有一个滤波器组，其中每个滤波器都有相同的带宽，就像在短时傅里叶变换或频谱图中，滤波器组，那么要保持滤波器之间的恒定频率重叠，你需要一个恒定的步长。然而，对于小波，步长应该像带宽一样与频率成比例。对于连续小波分析，最常见的间隔是2^(1/NV),NV是每个八度频域的过滤器数量，提高到整数次幂。作为比较，离散小波分析中专门使用的间隔是以2为底的整数次幂。

参见[2]对离散小波分析进行了深入的处理。下表总结了离散小波技术和连续小波技术之间的主要异同。对于离散技术，括号中提供了MATLAB中代表性算法的名称。

时频分析

由于小波同时在时间和频率上进行局部化，因此它们在许多应用中都很有用。对于连续小波分析，最常见的应用领域是时频分析。此外，CWT的下列特性使它对某些类型的信号特别有用。

为了理解哪些信号是连续小波分析的良好候选者，考虑一个双曲啁啾信号。

负载hyperbolicChirp图图(t，双曲chirp)轴紧包含(“时间(s)”) ylabel (“振幅”)标题(双曲唧唧喳喳的）

双曲啁啾是 $$\mathrm{罪}（\frac{15\pi }{0．8-t}）$$和 $$\mathrm{罪}（\frac{5\pi }{0．8-t}）$$．第一个组件在0.1秒到0.68秒之间激活，第二个组件在0.1秒到0.75秒之间激活。每个时刻的频率，或者这些分量的瞬时频率是 $$\frac{7．5}{（0．8-t{）}^2}$$和 $$\frac{2．5}{（0．8-t{）}^2}$$周期分别/秒。这意味着瞬时频率在附近非常低 $$t＝0$$并随着时间接近0.8秒而迅速增加。

获得啁啾信号的CWT，并将CWT与真实瞬时频率绘制为虚线。

Fs = 1/mean(diff(t));[cfs,f] = cwt(双波啁啾，Fs);helperHyperbolicChirpPlot (cfs f t)

我们看到，CWT显示了一个时频表示，准确地捕捉了双曲啁啾的瞬时频率。现在用使用恒定带宽滤波器的短时傅里叶变换分析相同的信号。使用默认窗口大小和重叠。

pspectrum (hyperbolchirp, 2048,的谱图）

请注意，当频率较低时，频谱图无法区分两个不同的瞬时频率，或者等效地，当两个瞬时频率接近时。基于输入长度2048，频谱图默认使用128个样本的窗口长度计算。在2048 Hz的采样率下，这相当于62.5 msec的时间分辨率。基于默认的Kaiser窗口，这将导致频率分辨率为41 Hz。这太大了，无法区分t=0附近的瞬时频率。如果我们将频率分辨率降低大约1/2，希望更好地分辨较低的频率，会发生什么?

pspectrum (hyperbolchirp, 2048,的谱图，“FrequencyResolution”, 20岁,.．.“OverlapPercent”, 90)

结果并不好，因为我们在较低频率上获得的任何频率分离都会导致较高的瞬时频率在时间上涂抹在一起。频谱图是一种强大的时频分析技术，事实上在许多应用中可以说是最优的，但像所有技术一样，它也有局限性。这种特殊类型的信号对固定带宽滤波器提出了重大挑战。理想情况下，您希望在低频处有较窄的频率支持的长时间响应，在高频处有较宽的频率支持的短时间响应。万博1manbetx小波变换恰恰提供了这一点，因此在这种情况下非常成功。另一种等效的说法是小波变换在较高频率有更好的时间分辨率，在较低频率有更好的频率分辨率。

在现实信号中，高频事件通常持续时间较短，而低频事件持续时间较长。例如，加载并绘制一个人体心电图(ECG)信号。数据采样频率为180hz。

负载wecgTm = 0:1/180:数字(wecg)*1/180-1/180;wecg情节(tm)网格在轴紧标题(“人体心电图”)包含(“时间(s)”) ylabel (“振幅”）

这些数据由缓慢变化的高频瞬态组成，这些高频瞬态代表了电脉冲在心脏中传播的过程。如果我们对及时准确地定位这些脉冲事件感兴趣，我们希望分析窗口(过滤器)很短。在这种情况下，我们愿意牺牲频率分辨率以获得更准确的时间定位。另一方面，为了识别低频振荡，最好有一个较长的分析窗口。

[cfs,f] = cwt(wecg,180);显示亮度图像(tm、f、abs (cfs))包含(“时间(s)”) ylabel (的频率(赫兹))轴xyCaxis ([0.025 0.25])“ECG数据的CWT”）

CWT显示30hz和37hz附近的稳态振荡以及表示心跳的瞬态事件(QRS复合体)。通过这种分析，您能够以相同的时频表示同时分析这两种现象。

另一个高频瞬变的例子是缓慢变化的分量，考虑尤利西斯号宇宙飞船从1993年12月4日21:00到1994年5月24日12:00在太阳南极上空每小时记录的太阳磁场大小的时间序列。参见[2第218-220页对该数据的完整描述。

负载solarMFmagnitudes情节(sm_dates, sm)轴紧标题(“太阳磁场大小”)包含(“日期”) ylabel (“级”）

原始数据显示了几个脉冲事件或冲击波的总体大致线性趋势。当快速太阳风或快速日冕物质抛射超过缓慢太阳风时，就会产生激波。在时间序列中，我们看到重大的冲击波结构大约在以下日期发生:2月11日，2月26日，3月10日，4月3日和4月23日。

对数据进行连续小波分析，并将CWT大小与时间序列数据绘制在同一图上。

sm_dates helperSolarMFDataPlot (sm)

CWT在时间序列中同时捕捉脉冲事件。然而，CWT也揭示了隐藏在时间序列中的数据的低频特征。例如，从1994年2月11日激波结构之前开始并延伸到激波结构之外，有一个低频稳态事件(接近0.04个周期/天)。这是由于日冕物质抛射(CME)事件可能发生在不同的太阳区域，与更近的激波同时到达尤利西斯飞船。

对于ECG和太阳磁场幅度的例子，我们有两个时间序列，它们来自非常不同的机制，产生类似的数据:具有短时间瞬态和较长时间低频振荡的信号。用相同的时频表示来描述两者是有利的。

本地化瞬变

时间序列数据中的瞬态事件通常是信息量最大的事件之一。瞬态事件可以指示数据生成机制中的突然变化，您希望及时检测和本地化这种变化。例子包括机械故障、传感器问题、经济时间序列中的金融“冲击”等等。作为一个例子，考虑下面的信号。

rng默认的；Dt = 0.001;T = 0:dt:1.5-dt;addNoise = 0.025*randn(size(t));x = cos(2 *π* 150 * t) * (t > = 0.1 & t < 0.5) +罪(2 *π* 200 * t) * 0.7 (t > & t < = 1.2);x = x+addNoise;X ([222 800]) = X ([222 800]) +[-2 2];图;情节(t。* 1000 x);包含(的毫秒）;ylabel (“振幅”）;

该信号由两个正弦分量组成，在222毫秒和800毫秒时突然打开和关闭，并伴有额外的“缺陷”。获得信号的CWT，并只绘制最佳尺度，或最高中心频率，小波系数。创建第二个轴并绘制原始时间数据。

CFS = cwt(x,1000，“爱”）;情节(t、abs (cfs (1:)),“线宽”(2) ylabel“级”甘氨胆酸)组(,“XTick”，[0.1 0.222 0.5 0.7 0.8 1.2]正确的情节(t x,“——”) ylabel (“振幅”)包含(“秒”)举行从

最优尺度的CWT系数可以定位数据中的所有突变。CWT系数显示正弦分量在222毫秒和800毫秒时打开和关闭以及定位正弦中的缺陷的峰值。由于CWT系数与数据具有相同的时间分辨率，为了检测数据中的瞬态变化，通常使用CWT精细尺度的震级。要最准确地定位这些瞬态，请使用最佳尺度(最高频率)系数。

锐化时频分析

任何使用滤波器的时频变换，如CWT中的小波，或短时傅里叶变换中的调制窗口，都必然会在时间和频率上抹黑信号的图像。信号能量在时间和频率上定位的不确定性来自于滤波器在时间和频率上的扩展。同步压缩是一种试图通过沿频率轴“压缩”变换来补偿这种涂抹的技术。为了用小波来说明这一点，获得由幅频调制(AM-FM)信号和18hz正弦波组成的信号的CWT。调幅-调频信号由方程定义

负载multicompsigSig = sig1+sig2;类(团体、sampfreq“爱”）

虽然CWT显示了18hz正弦波和AM-FM分量，但我们看到18hz正弦波在频率上确实是分散的。现在使用同步压缩来锐化频率估计。

[sst,f] = wsst(sig,sampfreq);图helperSSTPlot (sst、t、f)

同步压缩变换在很大程度上补偿了小波滤波器带来的时频扩展。同步压缩将时频相对较宽的峰值压缩为较窄的峰值脊．你实际上可以提取这些脊并从时频脊中重建各个分量。

(冰箱,iridge) = wsstridge (sst 5 f,“NumRidges”2);图;轮廓(t、f、abs (sst));网格在；标题(“山脊同步压缩变换”）;包含(的时间(秒)）;ylabel (“赫兹”）;持有在；情节(t,冰箱,“k——”，“线宽”2);持有从；

最后从时频脊重建分量的近似，并将结果与原始分量进行比较。

Xrec = iwsst(sst, ridge);helperPlotComponents (xrec sig1 sig2 t)

比较两种信号的时变频率内容

通常你有两个信号，它们可能在某种程度上是相关的。一个信号可以决定另一个信号的行为，或者两个信号可能只是因为两个信号的一些外部影响而相互关联。如果你得到了两个信号的CWT，你可以使用这些变换来得到小波相干性,时变相关性的度量。

每个时间瞬间和中心频率的小波相干性产生一个在0到1之间的相干值，它量化了两个信号之间的相关性强度。由于使用了复值小波，因此还可以使用相位信息来推断信号之间的超前-滞后关系。例如，考虑以下两个信号。

T = 0:0.001:2;X = cos(2*pi*10*t).*(t>=0.5 & t<1.1)+.．.因为(2 *π* 50 * t)。*(t>= 0.2 & t< 1.4)+0.25*randn(size(t));Y = sin(2*pi*10*t).*(t>=0.6 & t<1.2)+.．.罪(2 *π* 50 * t)。*(t>= 0.4 & t<1.6)+ 0.35*randn(size(t));图subplot(2,1,1) plot(t,X) ylabel(“振幅”) subplot(2,1,2) plot(t,Y)“振幅”)包含(“时间(s)”）

两个信号都由两个正弦波(10hz和50hz)构成白噪声。正弦波的时间支持略有不同。万博1manbetx注意，10hz和50hz组件Y相对于相应的组件滞后了1/4个周期X．绘制两个信号的小波相干性。

1000年图wcoherence (X, Y,,“PhaseDisplayThreshold”, 0.7)

小波相干性估计捕捉到信号在50和10 Hz附近高度相关，并且仅在正确的时间间隔内。图中黑色箭头表示信号之间的相位关系。这里箭头指向直线向上，表示两个信号之间的90度相位关系。这正是由信号中的余弦项决定的关系。

作为一个真实世界的例子，考虑一下厄尔尼诺区域3的数据和1871年至2003年底的去季节化全印度降雨指数。数据每月抽样一次。尼诺3时间序列记录了赤道太平洋西经90度至150度、北纬5度至南纬5度的每月海面温度异常情况。全印度降雨指数代表印度平均降雨量，单位为毫米，剔除了季节性成分。

负载ninoairdata图subplot(2,1,1) plot(datayear,nino)标题(“厄尔尼诺3区—海温异常”) ylabel (“度”)包含(“年”)轴紧子图(2,1,2)图(数据年，空气)轴紧标题(“去季节化全印度降雨指数”) ylabel (“毫米”)包含(“年”）

计算两个时间序列之间的小波相干性。设置相位显示阈值为0.7。如果需要以年为单位显示周期，请指定采样间隔为1/12年。

图wcoherence(尼诺、空气、年(1/12),“PhaseDisplayThreshold”, 0.7);

该图显示了与典型的2至7年厄尔尼诺周期相对应的时间范围内的强一致性区域。该图还显示，在这些时期，两个时间序列之间大约有3/8到1/2个周期的延迟。这表明，与南美洲海岸记录的厄尔尼诺现象相一致的海洋变暖时期与大约17000公里外印度的降雨量相关，但这种影响延迟了大约1/2个周期(1至3.5年)。

结论

在这个例子中，你学到了:

看到时频分析有关连续小波分析的更多突出主题和特色示例。

参考文献

[1]马拉特s.g.信号处理的小波之旅:稀疏方式．阿姆斯特丹第三版 ;波士顿:爱思唯尔/学术出版社，2009。

[2]珀西瓦尔，唐纳德·B和安德鲁·t·沃顿。时间序列分析的小波方法．剑桥:剑桥大学出版社，2000。https://doi.org/10.1017/CBO9780511841040。