如何对时间序列进行小波分析，得出其周期？

阿里云国内75折回扣微信号：monov8

阿里云国际，腾讯云国际，低至75折。AWS 93折免费开户实名账号代冲值优惠多多微信号：monov8 飞机：@monov6

从信号处理角度进行分析简单的时间序列直接做各种谱分析频谱包络谱平方包络谱功率谱倒谱等等比如一些简单的旋转机械振动时间序列信号

如果频谱不好分析那可以分析如下图所示的时间序列的时频谱

给个简单的模拟信号的例子 t = 0:1/2000:1-1/2000; dt = 1/2000; x1 = sin(50*pi*t).*exp(-50*pi*(t-0.2).^2); x2 = sin(50*pi*t).*exp(-100*pi*(t-0.5).^2); x3 = 2*cos(140*pi*t).*exp(-50*pi*(t-0.2).^2); x4 = 2*sin(140*pi*t).*exp(-80*pi*(t-0.8).^2); x = x1+x2+x3+x4; figure; plot(t,x) title('Superimposed Signal')

其连续小波变换时频谱如下

一个模拟的轴承内圈故障振动信号带有明显的周期性

相应的频谱如下红色虚线代表故障特征频率及相应的倍频

包络谱如下

看一下相应的CWT时频谱很明显能看出冲击性

还可以试试小波相干与交叉小波分析小波相干、交叉小波可以很好地反映两个不同时间序列变化之间的“相关性”。小波相干分析一般反映序列间周期性“变化趋势”的一致性但不直接反映变化周期的强度关系交叉小波分析一般反映序列间“共有周期”的强度。

此外如果时频谱线能量发散时频脊线模糊还可以试试同步压缩之类的算法时间序列信号处理系列-基于Python的同步压缩变换 - 哥廷根数学学派的文章 - 知乎时间序列信号处理系列-基于Python的同步压缩变换 - 知乎当时间序列信号中噪声较大时为了有利于周期性分析不可避免地要进行降噪前处理比如K-SVD降噪

样条框架降噪

Morlet小波降噪

当待分析的时间序列过于复杂时那可能要引入模态分解多分辨分析比如小波分解经验模态分解及其变体变分模态分解经验小波变换局部均值分解辛几何模态分解各种各样的自适应分解算法基于小波脊线的时间序列分解

好多同学都对各种模态分解方法的时间序列处理感兴趣那就随便说一下实际上时间序列通常由多个具有物理意义的分量组成在很多时候为了更容易的研究信号我们希望在与原始数据相同的时间尺度上单独研究这些分量中的一个或多个理想情况下我们希望这些经MRA分解到的多个分量在物理上是有意义的可容易解释的。多分辨率分析MRA通常与小波或小波包相关联但诸如经验模态分解EMD变分模态分解VMD等模态分解方法也可以构成MRA。先给一个简单的合成信号信号以1000Hz的频率采样1秒钟。 Fs = 1e3; t = 0:1/Fs:1-1/Fs; comp1 = cos(2*pi*200*t).*(t>0.7); comp2 = cos(2*pi*60*t).*(t>=0.1 & t<0.3); trend = sin(2*pi*1/2*t); rng default wgnNoise = 0.4*randn(size(t)); x = comp1+comp2+trend+wgnNoise; plot(t,x) xlabel('Seconds') ylabel('Amplitude') title('Synthetic Signal')

该信号由3个主要分量组成频率为 60Hz的时间局部振荡分量、频率为 200 Hz的时间局部振荡分量和趋势项分量。趋势项分量为正弦曲线频率为0.5Hz。60Hz的振荡分量发生在 0.1到 0.3 秒之间而 200Hz的振荡分量发生在 0.7 到 1 秒之间。但这些分量从时域波形中无法分辨因此进行频域变换。 xdft = fft(x); N = numel(x); xdft = xdft(1:numel(xdft)/2+1); freq = 0:Fs/N:Fs/2; plot(freq,20*log10(abs(xdft))) xlabel('Cycles/second') ylabel('dB') grid on

从频率中可以更容易地辨别振荡分量的频率但时间局部性信号却丢失。为了同时定位时间和频率信息使用连续小波变换进行分析。

从CWT时频谱图中可以看出60Hz和200Hz分量的时间范围但没有发现趋势项分量。为了分离出信号的分量并单独进行分析接下来使用多分辨分析直接在时域中进行相关操作。多分辨分析通过将信号分成不同分辨率的分量进而缩小分析范围而提取不同分辨率的信号分量相当于分解数据在不同时间尺度上的变化或等效地在不同频带上进行分析。首先采用离散小波变换的变体最大重叠离散小波变换对信号进行多分辨分析分解层数为8。关于最大重叠离散小波变换的相关内容请查看如下文献。

最大重叠离散小波变换的8层多分辨分析分解如下

如果从上向下看会看到所分解的分量变得越来越平滑即分量频率越来越低。回想一下原始信号包含3个主要分量一个 200 Hz 的高频振荡成分、一个 60 Hz 的低频振荡成分和一个趋势成分它们都被加性噪声破坏了。从D2 图中可以看出时间局部化的高频分量被分解出来而下面的两个图包含较低频率的振荡分量这是多分辨率分析的一个重要方面最后S8子图中包含了趋势项分量。除了小波多分辨分析经验模态分解 (EMD) 是一种所谓的数据自适应多分辨技术。 EMD 在不使用固定基函数的情况下递归地从数据中提取不同的分辨率成分关于EMD相关文献浩如烟海不做赘述了。EMD的多分辨分析分解如下所示

虽然MRA分解分量的数目不同但 EMD MRA和小波 MRA会产生相似的信号波形在 EMD MRA分解中高频振荡成分位于第1个本征模态函数中 (IMF1)低频振荡成分主要位于IMF2和IMF3中IMF6 中的趋势项分量与小波技术提取的趋势分量非常相似。自适应多分辨分析的另一种技术是变分模态分解 (VMD)VMD 从信号中提取固有模式函数或振荡模式并不使用固定基函数进行分析。EMD在时域上递归以逐步提取低频IMF分量而VMD 首先识别频域中的信号峰值并同时提取所有模式相关文献如下 Dragomiretskiy, Konstantin, and Dominique Zosso. “Variational Mode Decomposition.” IEEE Transactions on Signal Processing 62, no. 3 (February 2014): 531–44. https://doi.org/10.1109/TSP.2013.2288675. VMD的多分辨分析分解如下所示

由上图可知与小波和EMD类似VMD将3个分量基本分离了出来。还有一种数据自适应多分辨分析技术经验小波变换 (EWT) EWT根据分析信号的频率构造 Meye小波进而进行自适应小波之前写过EWT相关的内容经验小波变换在信号处理及轴承故障诊断中的应用 - 哥廷根数学学派的文章 - 知乎https://zhuanlan.zhihu.com/p/53 EWT的多分辨分析分解如下所示

与之前的EMD和小波MRA类似EWT分解出了相关的振荡分量用于执行分析的滤波器及其通带信息如下

下面考虑一段神户地震信号源于1995 年 1 月 16 日在澳大利亚霍巴特的塔斯马尼亚大学记录从 20:56:51 (GMT) 开始以 1 秒的间隔持续 51分钟。 figure plot(T,kobe) title('Kobe Earthquake Seismograph') ylabel('Vertical Acceleration (nm/s^2)') xlabel('Time') axis tight grid on

以最大重叠离散小波变换为例其8层MRA分解如下

从D4和D5子图中可以看出初级与延迟次级波分量地震波中的分量以不同的速度传播初级波比次级剪切波传播的更快。将信号分解为若干分量的目的通常是去除某些分量以减轻对信号分析的影响MRA技术的关键是重建原始信号的能力如下

每种方法的最大重建误差约为10^(-12) 或更小表明它们可以完美地对信号进行重建。在很多研究中我们对趋势项不感兴趣由于趋势项一般位于最后一个 MRA 分解分量中因此只需将该分量去除然后进行重建。

此外再删除第1个MRA分解分量看起来主要是噪声

在前面我们将趋势项删除然而在许多应用中趋势项可能是我们的主要研究部分因此可视化几种MRA方法所提取的趋势项分量。

根据以上的分析小波MRA技术可以更平滑且最准确地提取趋势项EMD提取了一个平滑的趋势项但它相对于真实趋势幅度发生了偏移而 VMD似乎比小波和EMD更偏向于提取振荡分量。在前面的示例中强调了多分辨分析在检测数据中的振荡分量和总体趋势中的作用然而MRA还可以定位和检测信号中的瞬态成分。为了说明这一点以1947年第一季度至 2011 年第四季度美国实际国内生产总值 (GDP) 数据垂直的黑线标志着“大缓和”的开始标志着从 1980 年代中期开始美国宏观经济波动性减弱的时期很难从原始数据中辨别出来。