写在前面
如题,基于自相关函数的标准化功率谱估计方法,通过python来实现,发现关于功率谱的教程貌似还挺少。
以下是理论依据:
以下直接来介绍具体的方法,这里使用的数据是Niño 3.4 SST Index,计算的数据来源为 HadISST1 dataset.
https://psl.noaa.gov/gcos_wgsp/Timeseries/Nino34/
主要是SST数据在 5S-5N and 170-120W的区域平均结果
Time Interval: Monthly
Time Coverage: 1870 to Sep 2021
Update Status: Periodically updated
References:
Rayner N. A., D. E. Parker, E. B. Horton, C. K. Folland, L. V. Alexander, D. P. Rowell, E. C. Kent, A. Kaplan, Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century, J. Geophys. Res., 108 (D14), 4407, doi:10.1029/2002JD002670, 2003.
以上网址还有可以下载其他的一些指数
功率谱简单介绍
功率谱分析是以傅里叶变换为基础的频域分析方法,其意义为将时间序列的总能量分解到不同频率上的分量,根据不同频率的波的方差贡献诊断出序列的主要周期,从而确定周期的主要频率,或者说这一段时间序列隐含的显著周期。
对于一个样本量为n的离散时间序列,样本时间序列为:
可以通过两种完全等价的方法进行功率谱估计
方法1、直接使用傅里叶变换
序列可以展开为傅里叶级数,
其中,为傅里叶系数,他们可以通过以下求得
其中,k为波数,不同波数k的功率谱值为:
方法2、通过自相关函数间接作出连续功率谱估计
可以通过谱密度与自相关函数互为傅里叶变换的重要性质来进行功率谱估计,下面我将基于这种方法来实现。
方法和步骤
1. 计算样本时间序列的自相关系数
其自相关系数为:
其中,m为最大滞后数,一般取n/3 ~ n/10
以下是具体的计算公式:
2. 计算粗谱
在得到自相关系数后,根据自相关系数计算原始谱,以下是计算公式:
其中,k为波数,取值从0到m,
k=0,1,...,m
在实际计算中,需要考虑端点特性,采用以下计算公式:
3. 计算平滑功率谱
这里使用hanning 平滑方法,许多包例如numpy和scipy也包含了hanning平滑方法,这里还是手算来实现,当作训练了
4.确定频率和周期
这没啥好说的
5. 对谱进行显著性检验
这里通过计算红噪声标准谱的方法进行检验,以下是红/白噪声标准谱的计算公式:
给定显著性水平,计算红/白噪声检验置信上限谱值:
自由度计算公式为:
在显著水平,自由度为下,比较和,如果≥,则认为k波数对应的周期波动是显著的
6. 制作并分析谱图
绘制功率谱密度估计值的曲线,横轴为频率,绘制置信检验曲线,寻找显著周期
以下是最终的检验结果
从图中我们可以看出,超过95%置信区间的最显著周期大概在40个月左右
总结
以上基于自相关函数的方法实现了时间序列的功率谱估计,相关测试数据和代码放到了GitHub网址上,感兴趣的可以去玩玩:
https://github.com/Blissful-Jasper/jianpu_record
★部分材料来自于书本以及课堂课件:现代气候统计诊断与预测技术
给我点在看的人
越来越好看