一、Python时间序列小波分析——实例分析

news2024/10/6 20:33:55

小波分析是在Fourier分析基础上发展起来的一种新的时频局部化分析方法。小波分析的基本思想是用一簇小波函数系来表示或逼近某一信号或函数。
小波分析原理涉及到傅里叶变换,并有多种小波变换,有点点小复杂。但是不会原理没关系,只要会应用并解释就可以。
在时间序列分析中,小波分析主要用于时间序列的消噪和滤波,信息量系数和分形维数的计算,突变点的检测和周期成分的识别以及多时间尺度的分析等。

小波分析通常以matlab软件为主,具体可以参考 小波系数等值线图和小波方差图绘制教学,视频教学可以参考 全网最简单的小波系数等值线图和小波方差图绘制小白教学。
matlab大多数人是不会用滴,而且安装、上手慢,不如试试Python实现小波分析。
最近在做项目中使用小波分析做时间序列分析,从零开始,学习过程还是有点痛苦,经过学习,已经初步能够解释说明输出图,特此记录一下,供大家参考。
在学习中曾经想在某宝找人做,结果要价600元,一下激发了我的征服欲,赚了自己600元,哈哈哈。
Python上手简单,实现方便。PyWavelets 是一个免费的开源库,用于 Python 中的小波变换。全网对于小波分析使用python进行实例分析解释说明几乎没有,下面一步一步的介绍,给大家省600元。

pip install pycwt   #安装

PyWavelets 的主要特点是:

  • 1D、2D 和 nD 正向和反向离散小波变换(DWT 和 IDWT)
  • 一维、二维和 nD 多级 DWT 和 IDWT
  • 一维和二维稳态小波变换(Undecimated Wavelet Transform)
  • 一维和二维小波包分解和重建
  • 一维连续小波变换
  • 计算小波和缩放函数的近似值
  • 超过 100 个内置小波滤波器并支持自定义小波
  • 单精度和双精度计算
  • 真实而复杂的计算
  • 与 Matlab Wavelet Toolbox ™ 兼容的结果
    PyWavelets给出了一个例子,我感觉解释说明不清晰。
from __future__ import division
import numpy
from matplotlib import pyplot

import pycwt as wavelet
from pycwt.helpers import find

url = '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  # In years

N = dat.size
t = numpy.arange(0, N) * dt + t0

p = 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 dataset

mother = 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 noise

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

power = (numpy.abs(wave)) ** 2
fft_power = numpy.abs(fft) ** 2
period = 1 / freqs

power /= scales[:, None]
signif, 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 / sig95

glbl_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.show()

在这里插入图片描述

后续将以自己的数据集对输出图进行解释说明,其中输出图包括小波系数实部等值线图,小波系数方差图,小波系数模的平方图,以及主周期趋势图,分析不同时间尺度下的时间周期。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/350451.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

单臂路由配置

单臂路由实现不同VLAN间通信 链路类型 交换机连接主机的端口为access链路 交换机连接路由器的端口为Trunk链路 子接口 路由器的物理接口可以被划分成多个逻辑接口 每个子接口对应一个VLAN网段的网关 我们需要知道的是单臂路由和虚拟机软件原理比较相似。他们都是依托于物理设备…

按键70秒,Root轻松得:Linux惊现高危漏洞

导读E安全11月16日讯 Linux被发现高危漏洞(CVE-2016-4484)&#xff0c;攻击者可以通过持续按下Enter键70秒钟来获取 root initramfs shell&#xff0c;进而破坏Linux boxes。漏洞存在于Linux流行变体中的统一密钥设置(LUKS)。通过访问shell&#xff0c;攻击者可以解密Linux机…

标签管理系统

电子墨水屏最大的不便利之处在于它是单色屏&#xff0c;显示速度不够快&#xff0c;因为显示内容的单色性&#xff0c;要求需要显示的数据处理好后&#xff0c;以单色位图的方式把数据信息发送到屏幕上才可以正常显示。 我们把这个制作显示内容的方式叫做模版&#xff0c;有了…

Zynq 裸机 PS + PL 双网口实现之 lwip 库文件修改

基于 xilinx vivado 2017.4 库文件 lwip141_v2_0 的修改&#xff1a; 添加对 PHY 芯片 ksz9031 的支持&#xff1b; 添加 SDK 中 LWIP 参数设置对话框 emio_options 选项&#xff1b; 添加 XPAR_GMII2RGMIICON_0N_ETH0_ADDR 和 XPAR_GMII2RGMIICON_0N_ETH1_ADDR 宏配置&#…

布隆过滤器和布谷鸟过滤器详解

今天和大家分享下布隆过滤器和布谷鸟过滤器 一.布隆过滤器 1.简单介绍 布隆过滤器是用于检索一个元素是否在一个集合中的算法&#xff0c;是一种用空间换时间的查询算法。 2.实现原理 布隆过滤器的存储结构是一个bitmap结构&#xff0c;初始值都是0&#xff0c;如下图所示&am…

内核模块(传参和依赖)

目录 一、模块传参 二、模块依赖 三、内核空间和用户空间 四、执行流 五、模块编程与应用编程的比较 六、内核接口头文件查询 七、小作业 一、模块传参 module_param(name,type,perm);//将指定的全局变量设置成模块参数 name:全局变量名 type&#xff1a; 使用符号 …

MFC消息机制

1.消息映射消息映射是一个将消息和成员函数相互关联的表。比如&#xff0c;框架窗口接收到一个鼠标左击消息&#xff0c;MFC将搜索该窗口的消息映射&#xff0c;如果存在一个处理WM_LBUTTTONDOWN消息的处理程序&#xff0c;然后就调用OnButtonDown。2.消息映射机制2.1 声明宏 写…

教务查询系统简介

教务查询系统简介 项目核心代码展示 service层如下&#xff1a; Teacher老师Service层&#xff1a; public interface TeacherService {//根据id更新老师信息void updateById(Integer id, TeacherCustom teacherCustom) throws Exception;//根据id删除老师信息void removeB…

SheetJS的通用电子表格对象简介和使用

简言 “通用电子表格格式”&#xff08;CSF&#xff09;是SheetJS使用的对象模型。 例如使用xlsx插件时&#xff0c;获得的excel文件数据对象就是依据这个模型设计的。 SheetJs通用电子表格对象 介绍 cdn导入&#xff1a; <script lang"javascript" src"ht…

uprobe 实战

观测数据源 目前按照我的理解&#xff0c;和trace相关的常用数据源–探针 大致分为四类。 内核 Trace point kprobe 用户程序 USDT uprobe 在用户程序中&#xff0c;USDT是所谓的静态Tracepoint。和内核代码中的Trace point类似。实现方式是在代码开发时&#xff0c;使用USDT…

Visual Studio开启clang-tidy代码检查

在CLion中有针对C的静态代码检查工具clang-tidy&#xff0c;感觉非常好用&#xff0c;能养成好的编码习惯&#xff0c;后来写Qt转入了VS平台&#xff0c;想要继续使用功能一致的clang-tidy体验&#xff0c;所以研究出来在VS中开启clang-tidy的方法。 版本&#xff1a;Visual S…

XML 基础知识 XXE 漏洞原理解析及实验(基础篇)

XML 介绍 XXE全称XML外部实体注入&#xff0c;所以在介绍XXE漏洞之前&#xff0c;先来说一说什么是XML以及为什么使用XML进而再介绍一下XML的结构。 XML全称 可拓展标记语言&#xff0c;与HTML相互配合后&#xff1a; HMTL用来显示数据 HTML的焦点在于数据的外观&#xff08;…

双网卡(有线和wifi)同时连接内网和外网

双网卡&#xff08;有线和wifi&#xff09;同时连接内网和外网 Win10技巧&#xff1a;如何修改有线/WiFi网络优先级&#xff1a;https://www.ithome.com/html/win10/253612.htm双网卡实现两个网络的自由访问&#xff1a;https://blog.51cto.com/ghostlan/1299090Linux服务器安…

【Linux】网络编程 - 基础概念

目录 一.OSI七层模型vsTCP/IP五层模型 1.一些周边概念 2.OSI七层模型 3.TCP/IP五层模型 4.网络传输流程图 二.什么是MAC地址 三.什么是IP/IP地址 1.什么是IP 2.什么是IP地址 四.什么是端口号 一.OSI七层模型vsTCP/IP五层模型 1.一些周边概念 局域网vs广域网 网络互…

LeetCode——2341. 数组能形成多少数对

一、题目 给你一个下标从 0 开始的整数数组 nums 。在一步操作中&#xff0c;你可以执行以下步骤&#xff1a; 从 nums 选出 两个 相等的 整数 从 nums 中移除这两个整数&#xff0c;形成一个 数对 请你在 nums 上多次执行此操作直到无法继续执行。 返回一个下标从 0 开始、…

亚马逊云科技重磅发布《亚马逊云科技汽车行业解决方案》

当今&#xff0c;随着万物智联、云计算等领域的高速发展&#xff0c;创新智能网联汽车和车路协同技术正在成为车企加速发展的关键途径&#xff0c;推动着汽车产品从出行代步工具向着“超级智能移动终端”快速转变。挑战无处不在&#xff0c;如何抢先预判&#xff1f;随着近年来…

31-Golang中的二维数组

二维数组的使用方式 使用方式一&#xff1a;先声明/定义再赋值 1.语法&#xff1a;var数组名 [大小] [大小]类型2.比如&#xff1a;var arr [2] [3]int,再赋值 package main import ("fmt" )func main() {//定义/声明数组var arr [4][6]int//赋初值arr[1][2] 1ar…

volatile,内存屏障

volatile的特性可见性: 对于其他线程是可见,假设线程1修改了volatile修饰的变量,那么线程2是可见的,并且是线程安全的重排序: 由于CPU执行的时候,指令在后面的会先执行,在指令层级的时候我们晓得volatile的特性后,我们就要去volatile是如何实现的,这个很重要&#xff01;&#…

金三银四面试必备的软件测试八股文,看完拿捏面试官

1、问&#xff1a;你在测试中发现了一个 bug&#xff0c;但是开发经理认为这不是一个 bug&#xff0c;你应该怎样解决&#xff1f;首先&#xff0c;将问题提交到缺陷管理库里面进行备案。然后&#xff0c;要获取判断的依据和标准&#xff1a; 根据需求说明书、产品说明、设计文…

蓝屏怎么办电脑蓝屏怎么办?蓝屏问题详细分析

蓝屏怎么办电脑蓝屏怎么办&#xff1f;最近很多小伙伴在咨询这个问题&#xff0c;其实电脑蓝屏了进不去&#xff0c;我们可以重新启动电脑&#xff0c;如果进入系统后还是直接蓝屏&#xff0c;那么你可以尝试一下&#xff0c;关机重启&#xff0c;然后在进入系统的时候&#xf…