Example 小波分析

简介

一般情况下，小波周期图水平轴对应原来时间序列的时间，垂直轴代表变化的周期，颜色代表变化周期的强度。该图里，黄色代表变化周期的高强度。

在最上面的图中表示 2 个变量的时间序列（Alongstream Velocity、Cross-Stream Velocity）,下面的 2 张图中为该 2 个变量的小波分析，随着时间序列的变化，变化周期的强度变化。

代码

from __future__ import division
import numpy
from matplotlib import pyplotimport pycwt as wavelet
from pycwt.helpers import findurl = 'http://paos.colorado.edu/research/wavelets/wave_idl/nino3sst.txt'
dat = numpy.genfromtxt(url, skip_header=19)
title = 'NINO3 Sea Surface Temperature'
label = 'NINO3 SST'
units = 'degC'
t0 = 1871.0
dt = 0.25  # 获取NINO数据# We also create a time array in years.
N = dat.size
t = numpy.arange(0, N) * dt + t0p = numpy.polyfit(t - t0, dat, 1)
dat_notrend = dat - numpy.polyval(p, t - t0)
std = dat_notrend.std()  # Standard deviation
var = std ** 2  # Variance
dat_norm = dat_notrend / std  # Normalized datasetmother = wavelet.Morlet(6)
s0 = 2 * dt  # Starting scale, in this case 2 * 0.25 years = 6 months
dj = 1 / 12  # Twelve sub-octaves per octaves
J = 7 / dj  # Seven powers of two with dj sub-octaves
alpha, _, _ = wavelet.ar1(dat)  # Lag-1 autocorrelation for red noisewave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J,mother)
iwave = wavelet.icwt(wave, scales, dt, dj, mother) * stdpower = (numpy.abs(wave)) ** 2
fft_power = numpy.abs(fft) ** 2
period = 1 / freqssignif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha,significance_level=0.95,wavelet=mother)
sig95 = numpy.ones([1, N]) * signif[:, None]
sig95 = power / sig95glbl_power = power.mean(axis=1)
dof = N - scales  # Correction for padding at edges
glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha,significance_level=0.95, dof=dof,wavelet=mother)sel = find((period >= 2) & (period < 8))
Cdelta = mother.cdelta
scale_avg = (scales * numpy.ones((N, 1))).transpose()
scale_avg = power / scale_avg  # As in Torrence and Compo (1998) equation 24
scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha,significance_level=0.95,dof=[scales[sel[0]],scales[sel[-1]]],wavelet=mother)# Prepare the figure
pyplot.close('all')
pyplot.ioff()
figprops = dict(figsize=(11, 8), dpi=72)
fig = pyplot.figure(**figprops)# First sub-plot, the original time series anomaly and inverse wavelet
# transform.
ax = pyplot.axes([0.1, 0.75, 0.65, 0.2])
ax.plot(t, iwave, '-', linewidth=1, color=[0.5, 0.5, 0.5])
ax.plot(t, dat, 'k', linewidth=1.5)
ax.set_title('a) {}'.format(title))
ax.set_ylabel(r'{} [{}]'.format(label, units))# Second sub-plot, the normalized wavelet power spectrum and significance
# level contour lines and cone of influece hatched area. Note that period
# scale is logarithmic.
bx = pyplot.axes([0.1, 0.37, 0.65, 0.28], sharex=ax)
levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16]
bx.contourf(t, numpy.log2(period), numpy.log2(power), numpy.log2(levels),extend='both', cmap=pyplot.cm.viridis)
extent = [t.min(), t.max(), 0, max(period)]
bx.contour(t, numpy.log2(period), sig95, [-99, 1], colors='k', linewidths=2,extent=extent)
bx.fill(numpy.concatenate([t, t[-1:] + dt, t[-1:] + dt,t[:1] - dt, t[:1] - dt]),numpy.concatenate([numpy.log2(coi), [1e-9], numpy.log2(period[-1:]),numpy.log2(period[-1:]), [1e-9]]),'k', alpha=0.3, hatch='x')
bx.set_title('b) {} Wavelet Power Spectrum ({})'.format(label, mother.name))
bx.set_ylabel('Period (years)')
#
Yticks = 2 ** numpy.arange(numpy.ceil(numpy.log2(period.min())),numpy.ceil(numpy.log2(period.max())))
bx.set_yticks(numpy.log2(Yticks))
bx.set_yticklabels(Yticks)# Third sub-plot, the global wavelet and Fourier power spectra and theoretical
# noise spectra. Note that period scale is logarithmic.
cx = pyplot.axes([0.77, 0.37, 0.2, 0.28], sharey=bx)
cx.plot(glbl_signif, numpy.log2(period), 'k--')
cx.plot(var * fft_theor, numpy.log2(period), '--', color='#cccccc')
cx.plot(var * fft_power, numpy.log2(1./fftfreqs), '-', color='#cccccc',linewidth=1.)
cx.plot(var * glbl_power, numpy.log2(period), 'k-', linewidth=1.5)
cx.set_title('c) Global Wavelet Spectrum')
cx.set_xlabel(r'Power [({})^2]'.format(units))
cx.set_xlim([0, glbl_power.max() + var])
cx.set_ylim(numpy.log2([period.min(), period.max()]))
cx.set_yticks(numpy.log2(Yticks))
cx.set_yticklabels(Yticks)
pyplot.setp(cx.get_yticklabels(), visible=False)# Fourth sub-plot, the scale averaged wavelet spectrum.
dx = pyplot.axes([0.1, 0.07, 0.65, 0.2], sharex=ax)
dx.axhline(scale_avg_signif, color='k', linestyle='--', linewidth=1.)
dx.plot(t, scale_avg, 'k-', linewidth=1.5)
dx.set_title('d) {}--{} year scale-averaged power'.format(2, 8))
dx.set_xlabel('Time (year)')
dx.set_ylabel(r'Average variance [{}]'.format(units))
ax.set_xlim([t.min(), t.max()])pyplot.savefig(fname="./wavelets.png", dpi=300)
pyplot.show()

实验结果

Example 等值线图

代码

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as npnp.random.seed(19680801)
npts = 200
ngridx = 100
ngridy = 200
x = np.random.uniform(-2, 2, npts)
y = np.random.uniform(-2, 2, npts)
z = x * np.exp(-x**2 - y**2)fig, (ax1, ax2) = plt.subplots(nrows=2)# Create grid values first.
xi = np.linspace(-2.1, 2.1, ngridx)
yi = np.linspace(-2.1, 2.1, ngridy)# Linearly interpolate the data (x, y) on a grid defined by (xi, yi).
triang = tri.Triangulation(x, y)
interpolator = tri.LinearTriInterpolator(triang, z)
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)# Note that scipy.interpolate provides means to interpolate data on a grid
# as well. The following would be an alternative to the four lines above:
# from scipy.interpolate import griddata
# zi = griddata((x, y), z, (xi[None, :], yi[:, None]), method='linear')ax1.contour(xi, yi, zi, levels=14, linewidths=0.5, colors='k')
cntr1 = ax1.contourf(xi, yi, zi, levels=14, cmap="RdBu_r")fig.colorbar(cntr1, ax=ax1)
ax1.plot(x, y, 'ko', ms=3)
ax1.set(xlim=(-2, 2), ylim=(-2, 2))
ax1.set_title('grid and contour (%d points, %d grid points)' %(npts, ngridx * ngridy))# Tricontourax2.tricontour(x, y, z, levels=14, linewidths=0.5, colors='k')
cntr2 = ax2.tricontourf(x, y, z, levels=14, cmap="RdBu_r")fig.colorbar(cntr2, ax=ax2)
ax2.plot(x, y, 'ko', ms=3)
ax2.set(xlim=(-2, 2), ylim=(-2, 2))
ax2.set_title('tricontour (%d points)' % npts)plt.subplots_adjust(hspace=0.5)plt.savefig(fname="./test4.png", dpi=300)
plt.show()

实验结果

参考资料

[1] Python 画图，点线图；
[2] Python 画图，柱状图；
[3] python matplotlib 画图（柱状图）总结；
[4] 数据可视化之美-动态图绘制（以Python为工具）；
[5] 图形美不胜收，25 个可视化案例，Matplotlib 始终都是数据可视化绕不开的 Python 库；
[6] Python matplotlib官方链接；
[7] 数据可视化之美+点、线、面组合（以Python为工具）；

python小波变换绘制频谱图和等线图相关推荐

用Python pyecharts v1.x 绘制图形（二）:折线图、折线面积图、散点图、雷达图、箱线图、词云图
文章目录关于pyecharts 折线图折线面积图散点图雷达图箱线图词云图其他关于pyecharts pyecharts是一个用于生成echart(百度开源的数据可视化javascrip ...
【python】seaborn绘制小提琴图和箱线图
文章目录箱线图增强箱线图小提琴图如果看过这篇用python给女朋友挑钻石,那么可以直接从箱线图开始读. seaborn是matplotlib的补充包,提供了一系列高颜值的figure,并且集成 ...
python 获取股市数据 baostock + 画K线图 mpl_finance
python 获取股市数据 baostock + 画K线图 mpl_finance 获取股票数据安装baostock库 baostock库的特性获取A股K线数据核心代码如下完整代码如下画K线 ...
数据可视化实验：python数据可视化-柱状图，条形图，直方图，饼图，棒图，散点图，气泡图，雷达图，箱线图，折线图
数据可视化实验:python数据可视化实验8-12:大数据可视化工具-python 目录 1柱状图 2条形图 3直方图 4饼图 5棒图 6散点图 7气泡图 8雷达图 9箱线图 10折线图 1柱状图 ...
R语言使用ggpubr包的ggline函数绘制各种漂亮形式的线图实战
R语言使用ggpubr包的ggline函数绘制各种漂亮形式的线图实战目录 R语言使用ggpubr包的ggline函数绘制各种漂亮形式的线图实战
使用python中的库matplotlib绘制箱线图(boxplot)
文章目录数据准备相关概念:Q1,Q2,Q3,IQR,QQ1,QQ3 绘制箱线图分析箱线图数据准备假设有如下数据: 我们要分别绘制出age这列数据的箱线图和%fat这列数据的箱线图. 相关概念 ...
【数据处理】python matplotlib 画箱线图；箱线图介绍；如何画箱线图
一.箱线图介绍假设一组数据有n个数,将它们从小到大排列,分为四等分.位于第25%(n+1)位置的数字是第一四分位数Q1.位于第50%(n+1)位置的数字是第二四分位数Q2,也是中位数.位于第75%( ...
数据可视化——R语言ggplot2包绘制精美的小提琴图(并箱线图或误差条图组合)
数据可视化--R语言ggplot2包绘制精美的小提琴图(并箱线图或误差条图组合) 概述:R语言使用ggplot2工具包绘制小提琴图.为了使数据表达更加丰富,同时将小提琴图与箱线图和误差条图相结合.另外 ...
软件设计模式（观察者模式）——模拟股票价格变动并绘制分时图和K线图
一.观察者模式简介: 首先看百度百科上对观察者模式的简介:观察者模式(Observer)完美的将观察者和被观察的对象分离开.举个例子,用户界面可以作为一个观察者,业务数据是被观察者,用户界面观察业务数 ...

python小波变换绘制频谱图和等线图

目录

Example 小波分析

简介

代码

实验结果

Example 等值线图

代码

实验结果

参考资料

python小波变换绘制频谱图和等线图相关推荐

最新文章

热门文章