Re: Harmonic analysis

[Billy Kessler <William.S.Kessler@xxxxxxxx>]

The following script will find any (single) harmonic of a time series
that consists of complete years. I also have a (more complicated)
script that does time series with blanks or partial years.

Billy K
-------------------------------------------------------------------
\canc mode verify
! Description: finds any annual harmonic. Does not plot, just sets up.

! arg1=variable name
! arg2=harmonic number wanted
! arg3=number of time gridpoints in a complete year
! arg4=number of complete years
! arg5=number of days in a year (default=365.2425)
!    -> deals with various length years (may be 360 or 365.25 days, etc)
!    but MUST BE COMPLETE YEARS. 

! need this so ferret knows what "t" is
set grid $1

! number of timesteps to be used
define symbol nts$1 `$3*$4`

! set frequency
let/q tunits$1=t
! tunits[l=1@ddf] = 1/number of sec in time axis unit (ferret does deriv in sec)
! therefore tunits[l=1@ddf]*86400 = 1 if time axis in days, = 86400 if in sec
let/q om1$1_$2= $2 / ($5"365.2425" * tunits$1[l=1@ddf]*86400)

let/q pi=4.*atan(1.)
let/q om1rad$1_$2=2.*pi*om1$1_$2

! terms which are products of trig terms and data
let/q dat$1=$1
let/q ysin$1_$2 =dat$1*(sin(om1rad$1_$2*t))
let/q sumys$1_$2=ysin$1_$2[l=1:($nts$1)@sum]
let/q ycos$1_$2=dat$1*(cos(om1rad$1_$2*t))
let/q sumyc$1_$2=ycos$1_$2[l=1:($nts$1)@sum]

! get harmonic components
! for complete years these are very simple. Note factor 2/(# timesteps)
let/q aa$1_$2 = (2/($nts$1)) * sumys$1_$2 
let/q bb$1_$2 = (2/($nts$1)) * sumyc$1_$2

! now get amp and phase (also make phase 0:360
let/q phi$1_$2=(360./(2.*pi))*atan2(aa$1_$2,bb$1_$2)
let/q phi$1_$2b=if (phi$1_$2 GE 0) then phi$1_$2 else phi$1_$2+360
let/q amp$1_$2=(aa$1_$2^2+bb$1_$2^2)^.5

! reconstruct a time series from the harmonic component
let/q recon$1_$2=aa$1_$2*sin(om1rad$1_$2*t)+bb$1_$2*cos(om1rad$1_$2*t)

! add definitions of mean and residual after mean and harmonic
let/q mean$1 = dat$1[l=1:($nts$1)@ave]
let/q resid$1_$2 = dat$1 - mean$1 - recon$1_$2

! add calculation of variance fraction
let/q var$1 = dat$1[l=1:($nts$1)@var]
let/q vfrac$1_$2 = 100*((amp$1_$2^2)/2)/var$1

! make good names for definition info below
define symbol dol1 $1
define symbol dol2 $2

say 
say Have defined:
say mean($dol1) = mean (fn of x,y,z)
say aa($dol1)_($dol2),bb($dol1)_($dol2) = harmonic coefficients (fn of x,y,z)
say amp($dol1)_$2 = harmonic amplitude (fn of x,y,z)
say phi($dol1)_$2,phi($dol1)_$2b = phase (-180:180 and 0:360 (b)) (fn of x,y,z)
say recon($dol1)_$2 = reconstructed time series (fn of x,y,z,t)
say resid($dol1)_$2 = residual after mean and harmonic (fn of x,y,z,t) 
say vfrac($dol1)_($dol2) = variance fraction (percent) (fn of x,y,z)
say
say Note: phase is relative to T0 of the original data's time axis!!
say --->>> sho axis [name of time axis]. Check T0.
say

set mode/last verify