蟒蛇的功率谱-显著水平

问题描述

用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年