Python中的功率谱 - 显著水平

Question

Python中的功率谱 - 显著水平

3

有很多关于如何使用Python计算功率谱的示例，例如使用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，但它测试的是分类数据具有给定频率的零假设。

- Emily

你可能需要自己编写代码，找到并使用方便的NumPy、SciPy或matplotlib函数。我猜这就是你所说的 - 来自www.atmos.washington.edu/~dennis/552_Notes_6b.pdf的6.2.6章节。 - wwii

你也可以在http://dsp.stackexchange.com/上提问。 - wwii

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Emily · Accepted Answer

我发现了以下公式，用于根据白噪声（或红噪声）的零假设计算功率谱中所有频谱峰值的显著性水平，参考文献[1]和[2]:

$\frac{\chi^{2}_{\phi,\alpha}{\phi}$ ,

使用白噪声（或红噪声）的理论功率谱 $Sp_{T}$ ，显著性水平 $\alpha$ 和自由度 $\phi$ 。功率谱的频率为 $\alpha$ ，其中h=0,1,...M， $n$ 是数据点数量。

在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., Practical Statistics, 1985, 公式（11-41）

[2]: Pankofsky, H.A.和Brier, G.W.，气象统计学的一些应用。宾夕法尼亚州立大学，1958年