Re: Power Spectrum

From: Dennis Shea (shea AT cgd.ucar.edu)
Date: Fri May 06 2005 - 19:04:19 MDT


> Could anyone tell me how the power spectrum is calculated in the
>function specx_anal€¯?
>

---
NCL usees FFTPACK to compute fourier info. The documentation
for the forward fft is at:
     ftp://ftp.ucar.edu/dsl/lib/fftpack/ezfftf.f
---

[1] detrend the series [optional] [2] taper the series [optional] [3] calculate the variance of the detrended/tapered series [4] forward fft on series [5] square the coef [periodogram ~2 dof] [6] smooth the periodogram estimates [7] normalize [6] so that the area under the curve is equal to the variance calculated in [3]

area_under_curve = SUM { S(f)*df*frac } where frac=1.0 except at the beginning and end values where frac=0.5 ====================

> If my understanding based on the manual is right, the power spectrum >should be return to the array > My concern is if spcx=|F(w)| or spcx= |F(w)|^2. or >spcx=(|F(w)|^2) / (N^2) ? Here F(w) is the Fourier Transform of x(t), | | means absolute value, N is the length of the F(w). > Neither ... The returned apc estimates have been normailized to get the area under the curve. FYI: The documentation for the forward fft is at: ftp://ftp.ucar.edu/dsl/lib/fftpack/ezfftf.f ---

You could do the above in NCL in a step-by-step process:

======================================================== In the last few days, I have been twice asked how NCL's specx_anal "works". It is pretty basic and for the most part can be well simulated by use of several of NCL's built-in functions.

[a] input a series x[*] with no missing values. [b] remove the mean an, optionally, detrend the series. http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/dtrend.html [c] taper the result of [b] http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/taper.html [d] Calculate the variance of [c] [d] perform an FFT [ezfftf] on [c] http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/fft.html [e] square the real (rc) and imaginary (ic) coef from [d]. This gives the periodogram. periodogram = rc^2 + ic^2 ; array syntax The periodgram has a bandwidth of 1/N (=df [delta_frequency]) but no statistical significance. [f] perform a running average over the periodogram values using either runave or wrunave. This creates a spectrum. http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/runave.html http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/wrunave.html A runave of say 5 pts creates a bandwidth of 5*df [approx] and approx 10 degrees of freedom. [g] normalize the area under the curve [f] to the input variance. Remember the first and last bandwidths are only 0.5*df. The numbers will not match 'spec_anal' exactly but should be very close. ++++++++++++++++++++++++++++++++++++++++++++++++++ good luck D

_______________________________________________ ncl-talk mailing list ncl-talk@ucar.edu http://mailman.ucar.edu/mailman/listinfo/ncl-talk



This archive was generated by hypermail 2b29 : Sun May 08 2005 - 22:05:34 MDT