当前位置:   article > 正文

小波分析: ​海温数据的时频域分解_海气大数据分析 小波分析

海气大数据分析 小波分析

模块安装

使用到的库是pycwt,安装非常简单,直接使用pip即可

pip install pycwt

复制

导入模块

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. import pycwt as wavelet
  4. from pycwt.helpers import find
  5. from matplotlib.pyplot import MultipleLocator

复制

下载数据

  1. # 从网页获取数据
  2. url = 'http://paos.colorado.edu/research/wavelets/wave_idl/nino3sst.txt'
  3. dat = np.genfromtxt(url, skip_header=19)
  4. print(dat.shape)

复制

(504,)

复制

设置参数

  1. title = 'NINO3 Sea Surface Temperature' # 标题
  2. label = 'NINO3 SST' # 标签
  3. units = 'degC' # 单位
  4. t0 = 1871.0 # 开始的时间,以年为单位
  5. dt = 0.25 # 采样间隔,以年为单位
  6. N = dat.size # 时间序列的长度
  7. t = np.arange(0, N) * dt + t0 # 构造时间序列数组

复制

去趋势&标准化

  1. p = np.polyfit(t - t0, dat, 1) # 线性拟合
  2. dat_notrend = dat - np.polyval(p, t - t0) # 去趋势
  3. std = dat_notrend.std() # 标准差
  4. var = std ** 2 # 方差
  5. dat_norm = dat_notrend / std # 标准化

复制

选择小波基函数

  1. mother = wavelet.Morlet(6) # Monther Wavelet: Morlet
  2. s0 = 2 * dt # Starting scale, in this case 2 * 0.25 years = 6 months
  3. dj = 0.25 # Twelve sub-octaves per octaves
  4. J = 7 / dj # Seven powers of two with dj sub-octaves
  5. alpha, _, _ = wavelet.ar1(dat) # Lag-1 autocorrelation for red noise

复制

计算小波变换(wave)和逆小波变换(iwave)

  1. wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J, mother)
  2. iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std

复制

计算小波能量谱

能量谱也叫能量谱密度,能量谱密度描述了信号或时间序列的能量如何随频率分布。能量谱是原信号傅立叶变换的平方。

  1. power = np.power(np.abs(wave),2)
  2. fft_power = np.power(np.abs(fft),2)
  3. period = 1 / freqs
  4. #power /= scales[:,None]

复制

计算小波功率谱

功率谱是功率谱密度函数,它定义为单位频带内的信号功率,是针对功率信号来说的。求功率谱就有了两种方法,分别叫做直接法和相关函数法:1、(傅立叶变换的平方)/(区间长度);2、自相关函数的傅里叶变换。

  1. signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha,
  2. significance_level=0.95,
  3. wavelet=mother)
  4. sig95 = np.ones([1, N]) * signif[:, None]
  5. sig95 = power / sig95

复制

  1. glbl_power = power.mean(axis=1)
  2. dof = N - scales # Correction for padding at edges
  3. glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha,
  4. significance_level=0.95, dof=dof,
  5. wavelet=mother)

复制

计算小波方差

  1. sel = find((period >= 2) & (period < 8))
  2. Cdelta = mother.cdelta
  3. scale_avg = (scales * np.ones((N, 1))).transpose()
  4. scale_avg = power / scale_avg # As in Torrence and Compo (1998) equation 24
  5. scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
  6. scale_avg_signif, tmp = wavelet.significance(
  7. var, dt, scales, 2, alpha,
  8. significance_level=0.95,
  9. dof=[scales[sel[0]],
  10. scales[sel[-1]]],
  11. wavelet=mother
  12. )

复制

绘制小波分析结果

  1. # Prepare the figure
  2. fig = plt.figure(figsize=(11, 8))
  3. # First sub-plot, the original time series anomaly and inverse wavelet transform.
  4. ax = plt.axes([0.1, 0.75, 0.65, 0.2])
  5. ax.plot(t, iwave, '-', linewidth=1, color=[0.5, 0.5, 0.5])
  6. ax.plot(t, dat, 'k', linewidth=1.5)
  7. ax.set_title('a) {}'.format(title))
  8. ax.set_ylabel(r'{} [{}]'.format(label, units))
  9. # Second sub-plot, the normalized wavelet power spectrum and significance
  10. # level contour lines and cone of influece hatched area. Note that period
  11. # scale is logarithmic.
  12. bx = plt.axes([0.1, 0.37, 0.65, 0.28], sharex=ax)
  13. levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16]
  14. bx.contourf(t, np.log2(period), np.log2(power), np.log2(levels), extend='both', cmap=plt.cm.viridis)
  15. extent = [t.min(), t.max(), 0, max(period)]
  16. bx.contour(t, np.log2(period), sig95, [-99, 1], colors='k', linewidths=1, extent=extent)
  17. bx.fill(np.concatenate([t, t[-1:] + dt, t[-1:] + dt,t[:1] - dt, t[:1] - dt]),
  18. np.concatenate([np.log2(coi), [1e-9], np.log2(period[-1:]),np.log2(period[-1:]), [1e-9]]),
  19. 'k', alpha=0.3, hatch='x')
  20. bx.set_title('b) {} Wavelet Power Spectrum ({})'.format(label, mother.name))
  21. bx.set_ylabel('Period (years)')
  22. Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),np.ceil(np.log2(period.max())))
  23. bx.set_yticks(np.log2(Yticks))
  24. bx.set_yticklabels(Yticks)
  25. # Third sub-plot, the global wavelet and Fourier power spectra and theoretical
  26. # noise spectra. Note that period scale is logarithmic.
  27. cx = plt.axes([0.77, 0.37, 0.2, 0.28], sharey=bx)
  28. cx.plot(glbl_signif, np.log2(period), 'k--')
  29. cx.plot(var * fft_theor, np.log2(period), '--', color='#cccccc')
  30. cx.plot(var * fft_power, np.log2(1/fftfreqs), '-', color='#cccccc',linewidth=1)
  31. cx.plot(var * glbl_power, np.log2(period), 'k-', linewidth=1.5)
  32. cx.set_title('c) Global Wavelet Spectrum')
  33. cx.set_xlabel(r'Power [({})^2]'.format(units))
  34. cx.set_xlim([0, 6])
  35. cx.set_ylim(np.log2([period.min(), period.max()]))
  36. cx.set_yticks(np.log2(Yticks))
  37. cx.set_yticklabels(Yticks)
  38. plt.setp(cx.get_yticklabels(), visible=False)
  39. # Fourth sub-plot, the scale averaged wavelet spectrum.
  40. dx = plt.axes([0.1, 0.07, 0.65, 0.2], sharex=ax)
  41. dx.axhline(scale_avg_signif, color='k', linestyle='--', linewidth=1.)
  42. dx.plot(t, scale_avg, 'k-', linewidth=1.5)
  43. dx.set_title('d) {}--{} year scale-averaged power'.format(2, 8))
  44. dx.set_xlabel('Time (year)')
  45. dx.set_ylabel(r'Average variance [{}]'.format(units))
  46. ax.set_xlim([t.min(), t.max()])
  47. plt.savefig('wavelet_analysis.jpg',dpi=600)
  48. plt.show()

复制

绘制小波方差

  1. Cdelta = mother.cdelta
  2. scale_avg = (scales * np.ones((N, 1))).transpose()
  3. scale_avg = power / scale_avg # As in Torrence and Compo (1998) equation 24
  4. scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
  5. scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha,
  6. significance_level=0.95,
  7. dof=[scales[0],scales[-1]],
  8. wavelet=mother)
  9. fig = plt.figure(figsize=(14, 6))
  10. ax = fig.add_subplot(111)
  11. ax.plot(t-t0, scale_avg, 'k-')
  12. ax.axhline(scale_avg_signif, color='k', linestyle='--', linewidth=1.)
  13. ax.set_xlim(0,126)
  14. ax.grid()
  15. x_major_locator=MultipleLocator(5)
  16. ax.xaxis.set_major_locator(x_major_locator)

复制

绘制小波系数

  1. fig = plt.figure(figsize=(14, 6))
  2. ax = fig.add_subplot(111)
  3. cf = plt.contourf(t, np.log2(period), np.log2(wave.real),
  4. np.arange(-4,4,0.5), extend='both', cmap=plt.cm.jet)
  5. #plt.clabel(cf,colors='k')
  6. Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),np.ceil(np.log2(period.max())))
  7. ax.set_yticks(np.log2(Yticks))
  8. ax.set_yticklabels(Yticks)
  9. ax.set_ylabel('Period (years)')
  10. plt.colorbar()

复制

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/码创造者/article/detail/856679
推荐阅读
相关标签
  

闽ICP备14008679号