时间序列如何分析周期性?（时间序列如何分析周期性问题）

从信号处理角度进行分析

简单的时间序列直接做各种谱分析（频谱，包络谱，平方包络谱，功率谱，倒谱等等）

比如一些简单的旋转机械振动时间序列信号

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

给个简单的模拟信号的例子

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的同步压缩变换 - 哥廷根数学学派的文章 - 知乎 https://zhuanlan.zhihu.com/p/554189692

当时间序列信号中噪声较大时，为了有利于周期性分析，不可避免地要进行降噪前处理

比如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 年代中期开始，美国宏观经济波动性减弱的时期，很难从原始数据中辨别出来。