3

There are a lot of examples how to calculate a power spectrum with python, e.g. 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])

But how can you calculate the 95% or 99% significance level of the power spectrum (null hypothesis: white noise)? I found scipy.stats.chisquare, but that tests the null hypothesis that the categorical data has the given frequencies.

Emily
  • 163
  • 1
  • 7
  • You might have to roll your own, finding and using convenient NumPy, SciPy or matplotlib *functions*. I assume this is what you are talking about - 6.2.6 of www.atmos.washington.edu/~dennis/552_Notes_6b.pdf . – wwii Oct 19 '16 at 19:01
  • You might also ask over at http://dsp.stackexchange.com/ – wwii Oct 19 '16 at 19:03

1 Answers1

1

I found following formula to calculate the significance level according to the null-hypothesis of white (or red) noise for all spectral peaks of the power spectrum in [1] and [2]:

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

with the theoretical power spectrum of white (or red) noise Sp_{T}, the significance level \alpha and the degrees of freedom \phi. The frequencies of the power spectrum are \alpha, for h=0,1,...M and n is the number of data points.

In 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, formula (11-41)

[2]: Pankofsky, H.A. and Brier, G.W., Some Applications of Statistics to Meteorology. Pennsylvania State Univ., 1958

Community
  • 1
  • 1
Emily
  • 163
  • 1
  • 7
  • could you help clarify: 1) `var=N.var(data)` should be `var=np.var(data)`? and 2), why `phi=(2*(n-1)-M2/.)/M` in the code, and in the equation it is shown as `phi=(2*n-M2/.)/M`? Thanks – Jason Aug 05 '21 at 14:24
  • I also would like to know that, I can't find the book in ref [1] anywhere. – paulo Jan 25 '22 at 21:43