Power spectrum in python - significance levels(蟒蛇的功率谱-显著水平)
本文介绍了蟒蛇的功率谱-显著水平的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!
问题描述
用python计算功率谱的例子有很多,例如Plotting power spectrum in python:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
data = np.random.rand(301) - 0.5
ps = np.abs(np.fft.fft(data))**2
time_step = 1 / 30
freqs = np.fft.fftfreq(data.size, time_step)
idx = np.argsort(freqs)
plt.plot(freqs[idx], ps[idx])
但如何计算功率谱的95%或99%显著水平(零假设:白噪声)?我找到了scipy.stats.chisquare,但这测试了分类数据具有给定频率的零假设。
推荐答案
我在[1]和[2]中根据功率谱的所有频谱峰值的白(或红)噪声的零假设,找到了以下公式来计算显著程度:
,
与白(或红)噪声的理论功率谱、显著水平和自由度。当h=0,1,...M时,功率谱的频率是数据点的数量。
在Python中:
import pylab as plt
import numpy as np
from scipy.stats import chi2
###
fft=np.fft.fft(data) ; n=len(fft)
abs=np.absolute(fft)**2
## frequencies (30min resolution)
f_u01=np.zeros(n/2+1,float)
f_u01=np.linspace(0,1,num=(n/2.+1))/(30*60*2)
### Variance of data as power spectrum of white noise
var=N.var(data)
### degrees of freedom
M=n/2
phi=(2*(n-1)-M/2.)/M
###values of chi-squared
chi_val_99 = chi2.isf(q=0.01/2, df=phi) #/2 for two-sided test
chi_val_95 = chi2.isf(q=0.05/2, df=phi)
### normalization of power spectrum with 1/n
plt.figure(figsize=(5,5))
plt.plot(fft[0:n/2],abs[0:n/2]/n, color='k')
plt.axhline(y=(var/n)*(chi_val_99/phi),color='0.4',linestyle='--')
plt.axhline(y=(var/n)*(chi_val_95/phi),color='0.4',linestyle='--')
[1]:Schönwiese,C.-D.,Praktische Statistik,1985,公式(11-41)
[2]潘科夫斯基,H.A.和布里尔,G.W.,统计学在气象学中的一些应用。宾夕法尼亚州立大学,1958年这篇关于蟒蛇的功率谱-显著水平的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持编程学习网!
沃梦达教程
本文标题为:蟒蛇的功率谱-显著水平
基础教程推荐
猜你喜欢
- 使 Python 脚本在 Windows 上运行而不指定“.py";延期 2022-01-01
- 使用Python匹配Stata加权xtil命令的确定方法? 2022-01-01
- 如何在 Python 中检测文件是否为二进制(非文本)文 2022-01-01
- 如何在Python中绘制多元函数? 2022-01-01
- 症状类型错误:无法确定关系的真值 2022-01-01
- Python 的 List 是如何实现的? 2022-01-01
- 合并具有多索引的两个数据帧 2022-01-01
- 使用 Google App Engine (Python) 将文件上传到 Google Cloud Storage 2022-01-01
- 将 YAML 文件转换为 python dict 2022-01-01
- 哪些 Python 包提供独立的事件系统? 2022-01-01