python小波分析，频率普分析——代码修改

写在前面：

做毕业论文的时候想要使用小波分析，开始到处找小波分析的代码，matlab代码居多但没有显著性检验，唯一找到有检验的这篇文章基于python进行小波分析，频率谱分析（作者：撼沧）是python写的，但是只能画出简单图形的python菜菜，多次试图修改代码都没有成功，抱着侥幸心理拿给老师看，果然被否了。终于终于把它给搞出来了，感谢撼沧的文章（鞠躬.jpg）。

2022.3月更新：另外还要感谢Chris Torrence刚刚去他的GitHub发现小波分析9个月前被他修改过，有需要的可以直接去他GitHub找代码和示例数据进行修改。（但更新后的代码和下文代码标注的行数差了几行，但应该不影响使用）wavelets/wave_python at master · chris-torrence/wavelets · GitHubhttps://github.com/chris-torrence/wavelets/tree/master/wave_python

关于图形的分析请看撼沧的文章，本文仅介绍使用自己的数据出图的方法。

正文：

1、 先将waveletAnalysis.py与waveletFunctions.py以及自己的的数据放到项目文件夹下，在pycharm中显示如下图： 

2、撼沧文章中提到的示例数据：共504个，时间从1871年到1996年，每年4个数据

      笔者(我本人)的数据：共41个，时间从1980年到2020年，每年一个数据

3、  后续只对waveletAnalysis.py中的代码进行修改即可，需要修改部分我标红了，修改后为橙色

22行        sst = np.loadtxt('sst_nino3.dat') # input SST time series

                标红部分替换为自己数据的文件名（或路径），对于我的数据则修改为：

                sst = np.loadtxt('precip.dat')

36行        dt = 0.25

                表示时间间隔，我的数据是从1980-2020年每年一个数据，那么我修改为：

                dt = 1

37行        time = np.arange(len(sst)) * dt + 1871.0 # construct time array

                修改成自己数据的起始时间，我的数据从1980年开始则修改为：

                time = np.arange(len(sst)) * dt + 1980.0

38行        xlim = ([1870, 2000]) # plotting range

                该句为绘出图像的坐标范围，而非数据范围，示例数据是从1871-1996年的，因此示

                例图像前后都是留有一点空白的。

                我不打算留白，所以我将代码修改为：

                xlim = ([1980, 2020])  

 40行       dj = 0.25   # this will do 4 sub-octaves per octave

                不懂是什么含义，但是我修改成与dt相等的数值，图看起来比较对。

                （2022.3月更新，我懂了！这个参数是修改阴影等值线的密集程度的！）

                dj = 1

41行        s0 = 2 * dt # this says start at a scale of 6 months

                指起始的时间尺度，也就是图a、b纵轴的起点，可以按需修改，如果你研究的内容

                仅有较大的年代际变化，那么可以从较大的数值开始。此处我修改为：

                s0 = 1

                也就是说，图中2前面的刻度是1，只不过没有写出。

42行        j1 = 7 / dj # this says do 7 powers-of-two with dj sub-octaves each

                可以调控图b图c的纵轴时间尺度的最大值，同样可以按需调整（虚线下面的值受边界

                效应影响不太准确，所以保证取的时间长度能低于锥形虚线即可）。我修改为：

                j1 = 5 / dj 

58行          avg = np.logical_and(scale >= 2, scale < 8)

                 该语句是用来修改图c的，图c是小波分析的时间域，修改2和8是在修改窗口滑动平均的

                尺度，要根据自己数据的周期来进行修改，我个人理解就是一个滤波的作用，可以显示

                出周期信号强的年份。在我所使用的数据中是没有进行修改的，因为当时没用到图c。

102行      plt.plot(time, coi, 'k')

                该命令用于绘制影响锥曲线，我运行时该行报错，所以进行了修改

                plt.plot(time, coi[0:len(sst)], 'k',ls='--') # ls='--'是为了修改成虚线形式更美观

                （更新注：我现在使用python3.9该句没有报错，没报错时可以不用修改，仍然想用虚线

                就添加,ls='--'即可）

除了以上必要的修改，还有图表的标题以及横纵坐标，找到plt.xlabel、plt.ylabel、plt.title来修改相应的横坐标，纵坐标以及标题就可以啦。

我的数据出图是这样的（去掉了一张）：

以上就是修改的全部代码，可以按照自己的数据需求修改，我运行时由于matplotlib的版本问题还会有这样的提示：MatplotlibDeprecationWarning: The 'basey' parameter of __init__() has been renamed 'base' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.        plt3.set_yscale('log', basey=2, subsy=None)

可以将所有的basey和subsy分别修改成base和subs

中间可能有一些不太明晰之处，若有小伙伴懂得，还望不吝赐教。但是起码能改出自己的图嘿嘿嘿，改了短短几行代码，真是磨了我好久，我太菜了。此教程仅供比我还菜的同学参考使用哈哈哈，如果对你有用，请点个赞吧，嘿嘿。

我使用的precip.dat在这里，需要可自行下载precip.dat