[ferret_users] Fitting a red noise spectrum to a power spectrum

[Alejandro Ludert <aludert@xxxxxxxxxx>]

Hello Ferret Community. I require some help Fitting a red noise spectrum to my calculated power spectrum of a time series.

I have calculated the power spectrum of a time series using the method described in the .jnl file ef_fft_demo.jnl. 

Now I need to fit a red noise spectrum to my plot and I am having a bit of trouble doing it. Here is a copy of the .jnl file I'm using:

\cancel mode verify

! ef_fft_demo.jnl

! *am* 4/99

!set mode verify

set mode metafile

! Define viewports----------------------------------------------

!define the a col

define view/xlim=0.,.49/ylim=.67,1/text=0.1 a1

define view/xlim=0.,.49/ylim=.34,.66/text=0.1 a2

define view/xlim=0.,.49/ylim=0.,.33/text=0.1 a3

!define b col

define view/xlim=.50,1/ylim=.67,1/text=0.1 b1

define view/xlim=.50,1/ylim=.34,.66/text=0.1 b2

define view/xlim=.50,1/ylim=0.,.33/text=0.1 b3

!----------------------------------------------------------------

MESSAGE

use/regulart seof.ts.nc

cancel mode logo

! S-EOF 1 obs----------------------------------------------------

! Define the time series at a point in space.

! Set the FFT, using explicit time specification.

! Plot the amplitude spectrum vs frequency.

LET FFT_uwndtim = ts[i=1, j=1, k=1]

LET FFT_uwndfft = ffta(FFT_uwndtim[l=1:24])

!  Get the frequency increment used in the FFT.

LET FFT_nf = `FFT_uwndfft,return=lend`

LET FFT_nyquist = 0.5

LET FFT_freq1 = FFT_nyquist/ FFT_nf

!  Define a frequency axis.

DEFINE AXIS/T=`FFT_freq1`:`FFT_nyquist`:`FFT_freq1` faxis

DEFINE GRID/T=faxis gfftfreq

LET a = T[g=gfftfreq]

!  Define the period from the frequency axis.

LET per = 1./a

!  Plot as a "Y VS X" plot, showing the first 24 months where the most energy is.

!  The PPL ccommands clean up the plot appearance.

SET VIEW ul

PLOT/VS/LINE/HLIMITS=0:24:2/xlimits=0:2/ylimits=0:13/nolab/SET_UP per[l=1:`FFT_nf`], FFT_uwndfft

PPL XFOR (I2)

PPL XLAB Period, years/cycle

PPL YLAB

ppl axtype,2,1

ppl xlab Period(years); ppl ylab Power; ppl title "S-EOF 1 OBS"

PPL PLOT

! Calculate the red noise spectrum

let var = ts[i=1, j=1, k=1, l=@var]

let p = ts[i=1, j=1, k=1, l=1:24]

let q = ts[i=1, j=1, k=1, l=@shf:-1]

go variance

let var_white = ((1- correl^2.)*var^2.)^0.5

let s1 = (4*var_white^2.)/24

let s2 = 1+correl^2.- 2*correl*cos(2*3.14159*per/24)

let s = s1/s2

!set view ur

PLOT/VS/LINE/HLIMITS=0:24:2/l=1:`FFT_nf`/nolab/SET_UP per[l=1:`FFT_nf`], s

PPL PLOT

! clean up

! restore plot state 

CANCEL SYMBOL FFT_*

cancel mode metafile

sp gksm2ps -p portrait -l cps -o t.eps metafile.plt

The error I get is: 

 *** NOTE: Ambiguous coordinates on T axis: 1+CORREL^2.- 2*CORREL*COS(2*3.14159*PER/24)

 **ERROR: inconsistent sizes of data regions: T axis

          Q has 12 points (L=01:12)

          _expression_ has 24 points (L=01:24)

PLOT/VS/LINE/HLIMITS=0:24:2/l=1:12/nolab/SET_UP per[l=1:12], s

Command file, command group, or REPEAT execution aborted

I can't get the power spectrum density function s to be properly evaluated at each point and then plotted over my previously calculated power spectrum. It seems that the problem is related to the correlation, when i replace the numerical values in the equations for s1 and s1 it works a bit better, but still not well. When I do that The spectrum has an improper magnitude and I'm not sure how to re scale it.

The formulas for the white noise variance and the spectral density function are equations 8.46 and 8.77 from Wilks second edition.

As a very new user to Ferret, your help will be greatly appreciated.

Thank you in advance for your time, Alejandro Ludert.

Attachment: seof.ts.nc
Description: Cdf file