这个例子展示了如何执行和解释连续小波分析。这个例子可以帮助您回答一些常见的问题,例如:连续小波分析和离散小波分析的区别是什么?为什么连续小波分析中的频率或尺度是对数间隔的?在什么类型的信号分析问题中连续小波技术特别有用?
人们经常问的一个问题是:什么是连续关于连续小波分析?毕竟,你大概对在电脑上做小波分析感兴趣而不是用纸和笔,而在电脑上,没有什么是真正连续的。当这个词连续小波分析在科学计算环境中使用,这意味着每倍频程有一个以上小波或频率加倍的小波分析技术,其中小波之间的时间偏移为一个样本。由此产生的连续小波变换(CWT)具有两个在应用中非常有用的特性:
与离散小波技术相比,信号的频率内容被更精细地捕获。
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所有小波滤波器的性质,即它们在频率或带宽上的扩展与它们的中心频率成正比。换句话说,小波滤波器在高频时比在低频时更宽。由于时间和频率之间的互反关系,这意味着高频小波比低频小波更好地在时间上本地化(有更小的时间支持)。万博1manbetx为了看到这一点,从滤波器组中提取中心频率和时域小波。绘制两个小波的实部,一个高频和一个低频,以作比较。
psi=小波(fb);F=中心频率(fb);图形图(实(psi(9,:))ylabel)(“振幅”) yyaxis正确的绘图(真实(磅/平方英寸(16年底:))ylabel(“振幅”)轴紧S1 =“中心频率小波”;S2 = [num2str (F (9),' % 1.2 f ')”和“num2str(F(16号端),' % 1.2 f ')“周期/样本”];标题({S1;S2})包含(“样本”)传说(“0.25 c / s波”,“0.01 c / s波”)
在频率响应图中,中心频率不超过0.45 cycles/sample。找出滤波器组在一个八度或双倍频率中有多少小波滤波器。确认数字等于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)符号(F)格在ylabel (“周期/样本”)标题(“半对数的规模”)
图表明,小波中心频率不是线性间隔,这是一般情况下,与其他滤波器组。具体来说,中心频率呈指数递减,高中心频率之间的步长大于低中心频率之间的步长。为什么这对小波有意义?回想一下小波是常数q滤波器,这意味着它们的带宽与中心频率成比例。如果你有一个滤波器组每个滤波器都有相同的带宽,比如短时间傅里叶变换,或者光谱图,滤波器组,那么要保持滤波器之间的频率重叠你需要一个恒定的步长。然而,对于小波,步长应该与频率成正比,就像带宽一样。对于连续小波分析,最常见的间距是基数2^(1/)NV),在哪里NV是每倍频程的滤波器数,提升为整数次幂。比较而言,离散小波分析中专用的间距是提升为整数次幂的基数2。
参见[2.对离散小波分析进行了深入的研究。下表总结了离散小波和连续小波技术的主要异同点。对于离散技术,MATLAB中有代表性的算法的名称用括号表示。
由于小波同时在时间和频率上进行局部化,因此它们在许多应用中都很有用。对于连续小波分析,最常见的应用领域是时频分析。此外,连续小波变换的以下特性使它对某些类型的信号特别有用。
小波滤波器的常q特性,即高频小波持续时间短,低频小波持续时间长。
连续分析中小波滤波器间的单样本时移
为了理解哪些信号是连续小波分析的好候选信号,考虑一个双曲啁啾信号。
负载hyperbolicChirp图绘制(t, hyperbolchirp)轴紧包含(“时间(s)”) ylabel (“振幅”)标题(“双曲线啁啾”)
双曲啁啾是 和 。第一个组件在0.1至0.68秒期间处于激活状态,第二个组件在0.1至0.75秒期间处于激活状态。这些组件的每个瞬间的频率或瞬时频率如下所示: 和 循环/秒。这意味着瞬时频率在附近非常低 并随着时间接近0.8秒迅速增加。
获取啁啾信号的CWT,并绘制CWT以及绘制为虚线的真实瞬时频率。
Fs = 1 /意味着(diff (t));[cfs f] = cwt (hyperbolchirp Fs);helperHyperbolicChirpPlot (cfs f t)
我们看到CWT显示了一种时频表示,它准确地捕获了双曲啁啾的瞬时频率。现在用短时间傅里叶变换对相同的信号进行分析,短时间傅里叶变换采用恒定带宽滤波器。使用默认的窗口大小和重叠。
pspectrum (hyperbolchirp, 2048,的谱图)
请注意,当频率较低时,谱图无法区分两个不同的瞬时频率,或者当两个瞬时频率非常接近时,谱图无法区分两个不同的瞬时频率。根据输入长度2048,谱图默认使用128个样本的窗长进行计算。在2048hz的采样率下,这相当于62.5 msec的时间分辨率。基于默认的Kaiser窗口,这导致频率分辨率为41 Hz。这太大了,无法区分t=0附近的瞬时频率。如果我们为了更好地分辨低频而将频率分辨率降低大约1/2会发生什么?
pspectrum (hyperbolchirp, 2048,的谱图,“频率分辨率”, 20岁,...“重叠百分比”,90)
结果并不是更好,因为我们在较低频率下获得的任何频率分离都会导致较高的瞬时频率在时间上混在一起。频谱图是一种强大的时频分析技术,事实上在许多应用中可以说是最优的,但与所有技术一样,它也有局限性。这种特殊类型的信号对固定带宽滤波器提出了重大挑战。理想情况下,您希望在低频时具有窄频率支持的长时间响应,在高频时具有宽频率支持的短时间响应。小波变换正好提供了这一点,因此在这种情况下非常成功。说明这一点的等效方法是,小波变换在较高频率下具有更好的时间分辨率,在较低频率下具有更好的频率分辨率。万博1manbetx
在真实世界的信号中,高频事件通常持续时间较短,而低频事件持续时间较长。例如,加载并绘制人体心电图(ECG)信号。数据以180hz采样。
负载wecgtm=0:1/180:numel(wecg)*1/180-1/180;绘图(tm,wecg)网格在轴紧标题(“人体心电图”)包含(“时间(s)”) ylabel (“振幅”)
这些数据由缓慢变化的部分组成,这些部分被高频瞬变信号打断,这些高频瞬变信号代表着通过心脏传播的电脉冲。如果我们对这些脉冲事件在时间上的精确局部化感兴趣,我们希望分析窗口(过滤器)是短的。为了更精确的时间定位,我们愿意牺牲频率分辨率。另一方面,为了识别低频振荡,最好使用较长的分析窗口。
[cfs f] = cwt (wecg, 180);显示亮度图像(tm、f、abs (cfs))包含(“时间(s)”) ylabel (‘频率(Hz)’)轴xycaxis([0.025 - 0.25])标题(“心电图数据的CWT”)
CWT显示近30hz和37hz的稳态振荡以及表示心跳的瞬态事件(QRS复合体)。通过这种分析,您可以在相同的时频表示中同时分析这两种现象。
作为高频瞬变标点缓慢变化分量的另一个例子,考虑由尤利西斯飞船在1993年12月4日从21:00 UT到1994年5月24日的12:00 UT每小时记录太阳太阳磁场的时间系列。看[2.第218-220页,以获得该数据的完整描述。
负载太阳辐射强度绘图(sm_日期,sm)轴紧标题(《太阳磁场强度》)包含(“日期”) ylabel (“震级”)
原始时间数据显示了几个脉冲事件或冲击波的总体大致线性趋势。当快速太阳风或快速日冕物质抛射超过缓慢太阳风时,会产生冲击波。在时间序列中,我们看到重大冲击波结构大约发生在以下日期:2010年2月11日、2月26日h 10、4月3日和4月23日。
对数据进行连续小波分析,并将CWT幅度与时间序列数据绘制在同一图上。
sm_dates helperSolarMFDataPlot (sm)
CWT在脉冲事件在时间序列中发生的同时捕获脉冲事件。然而,CWT也揭示了隐藏在时间序列中的数据的低频特征。例如,从1994年2月11日的激波结构之前开始并延伸到激波结构之外,出现了一个低频稳态事件(接近0.04周期/天)。这是由于日冕物质抛射(CME)事件可能发生在不同的太阳区域,与更近的激波同时到达尤利西斯宇宙飞船。
在心电图和太阳磁场量级的例子中,我们有两个来自非常不同的机制的时间序列产生相似的数据:既有短时瞬变信号,也有长时低频振荡信号。用相同的时频表示来描述两者是有利的。
时间序列数据中的瞬态事件通常是信息量最大的事件之一。瞬态事件可能表示数据生成机制中的突然变化,您希望及时检测并定位该变化。例子包括机械故障、传感器问题、经济时间序列中的金融“冲击”等。作为一个例子,考虑下面的信号。
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毫秒时突然开启和关闭,并伴有额外的“缺陷”。获得信号的小波变换,并仅绘制最小尺度,或最高中心频率,小波系数。创建第二个轴并绘制原始时间数据。
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)信号和18赫兹正弦信号组成的信号的CWT。调幅-调频信号由公式定义
负载multicompsigsig = sig1 + sig2;类(团体、sampfreq“爱”)
虽然CWT显示了18-Hz的正弦波以及AM-FM分量,我们看到18-Hz的正弦波肯定出现在频率上的扩散。现在使用同步压缩来锐化频率估计。
海温,[f] =墓场(sig, sampfreq);图helperSSTPlot (sst、t、f)
同步压缩变换在很大程度上补偿了小波滤波器引入的时频扩展。同步压缩将相对较宽的时频峰值压缩为较窄的峰值脊.你可以提取这些脊然后从时频脊中重建单个的分量。
(冰箱,iridge) = wsstridge (sst 5 f,“NumRidges”,2); 图形轮廓(t、f、abs(sst));网格在;标题(“带脊线的同步压缩变换”);xlabel(的时间(秒));ylabel (“赫兹”);持有在;绘图(t、冰箱、,“k——”,“线宽”2);持有从;
最后根据时频脊线重构近似分量,并将结果与原始分量进行比较。
iridge xrec = iwsst (sst);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)“振幅”)包含(“时间(s)”)
两个信号都是由两个正弦波(10hz和50hz)组成的白噪声。正弦波的时间支持略有不同。万博1manbetx注意,10赫兹和50赫兹的分量Y
相对于中相应的组分滞后1/4个周期X
.绘制两个信号的小波相干性。
1000年图wcoherence (X, Y,,“阶段显示阈值”, 0.7)
小波相干估计捕捉到信号在50和10赫兹附近高度相关,并且只在正确的时间间隔内。图上的暗箭头表示信号之间的相位关系。这里的箭头指向上方,这表明了两个信号之间的90度相位关系。这正是由信号中的cos -sin项所决定的关系。
作为一个真实的例子,考虑厄尔尼诺地区3的数据和从1871年到2003年底的季节性全印度降雨量指数。这些数据每月抽样一次。尼诺3时间序列是在赤道太平洋西经90度至150度,北纬5度至南纬5度记录的每月摄氏度海面温度异常记录。全印度降雨量指数代表了印度的平均降雨量,以毫米为单位,剔除了季节性成分。
负载ninoairdata图subplot(2,1,1) plot(datyear,nino) title(“厄尔尼诺3区—海温异常”) ylabel (“度”)包含(“年”)轴紧次要情节(2,1,2)情节(datayear、空气)轴紧标题(“取消了印度所有降雨指数的季节性”) ylabel (“毫米”)包含(“年”)
计算两个时间序列之间的小波相干性。将相位显示阈值设置为0.7。要显示以年为单位的周期,指定采样间隔为一年的1/12。
图wcoherence(尼诺、空气、年(1/12),“阶段显示阈值”, 0.7);
这幅图显示了强相干的时间局部区域,其发生周期与典型的2至7年厄尔尼诺周期相对应。该图还显示,在这些时期,两个时间序列之间大约有3/8到1/2个周期延迟。这表明,与南美洲海岸记录的厄尔尼诺相一致的海洋变暖时期,与17000公里以外的印度的降雨量相关,但这种影响大约延迟了1/2个周期(1至3.5年)。
在这个例子中,你学到:
连续小波分析和离散小波分析的区别。
CWT的常数q和单样本时移特性如何让你同时分析非常不同的信号结构。
一种使用连续小波变换的时频重分配技术。
如何使用CWT比较两个信号的频率内容。
看见时频分析有关连续小波分析的更多突出主题和特色示例。
马拉特,s.g.。信号处理的小波变换:稀疏方法.第三。阿姆斯特丹 ;波士顿:爱思唯尔/学术出版社,2009。
Percival, Donald B.和Andrew T. Walden。时间序列分析的小波方法.剑桥:剑桥大学出版社,2000。https://doi.org/10.1017/CBO9780511841040。
类
|cwtfilterbank
|wcoherence
|墓场
|wsstridge
|pspectrum
(信号处理工具箱)