连续与离散小波分析

人们经常会问的一个问题是:什么是连续关于连续小波分析?毕竟，你大概对在计算机上做小波分析感兴趣，而不是用铅笔和纸，在计算机中，没有什么是真正连续的。当这个词连续小波分析在科学计算设置中使用，它是指每八度有多个小波或频率加倍的小波分析技术，其中小波之间的时间移动是一个样本。这使得得到的连续小波变换(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个小波滤波器VoicesPerOctave)．画出小波滤波器的频率响应。

freqz(神奇动物)

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

ψ=小波(神奇动物);F = centerFrequencies (fb);图绘制(真实(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

你可以在网站上读到关于连续小波分析和离散小波分析之间区别的更详细的解释连续和离散小波变换．

对数间隔中心频率

小波分析的一个方面，人们可能会发现有点令人困惑的是滤波器的对数间距。

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

次要情节(2,1,1)情节(F) ylabel (“周期/样本”)标题(小波中心频率的网格)在次要情节(2,1,2)semilogy (F)网格在ylabel (“周期/样本”)标题(“半对数的规模”）

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

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

时频分析

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

为了理解哪种信号适合连续小波分析，考虑一个双曲啁啾信号。

负载hyperbolicChirp图绘制(t, hyperbolchirp)轴紧包含(“时间(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 /意味着(diff (t));[cfs f] = cwt (hyperbolchirp 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)

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

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

负载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])标题(“心电数据的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(大小(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]) yyaxis正确的情节(t x,“——”) ylabel (“振幅”)包含(“秒”)举行从

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

提高时频分析

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

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

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

海温,[f] =墓场(sig, sampfreq);图helperSSTPlot (sst、t、f)

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

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

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

iridge xrec = iwsst (sst);helperPlotComponents (xrec sig1 sig2 t)

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

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

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

t = 0:0.001:2;X = cos(2* *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* *10*t).*(t>=0.6 & t<1.2)+.．.罪(2 *π* 50 * t)。*(t>= 0.4 & t<1.6)+ 0.35*randn(size(t));figure subplot(2,1,1) plot(t,X) ylabel(“振幅”subplot(2,1,2) plot(t,Y)“振幅”)包含(“时间(s)”）

这两个信号都是由两个正弦波(10hz和50hz)构成的白噪声。正弦波的时间支撑点略有不同。注意，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度的海域，以摄氏度为单位记录的每月海表温度异常。全印度降雨量指数表示印度平均降雨量，单位为毫米，剔除了季节性成分。

负载ninoairdataFigure subplot(2,1,1) plot(datayear,nino) title(“厄尔尼诺3区—海温异常”) ylabel (“度”)包含(“年”)轴紧次要情节(2,1,2)情节(datayear、空气)轴紧标题(“去季节性的全印度降雨指数”) ylabel (“毫米”)包含(“年”）

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

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

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

结论

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

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

参考文献

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

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