\cancel mode verify
! ef_fft_demo.jnl
! *am* 4/99
set mode verify
! Description: Demonstration of external functions for FFT:  
! FFTA computes amplitude spectrum
! FFTP computes phase.

! The variable must have the explicit time specification.
! Planned upgrades to FFTA and FFTP:
!  - Allow the function to operate on the full time axis from the 
!    input grid without requiring the explicit time specification.
!  - Give units for the frequency axis of "1./(time units)" , where the
!    time units are taken from the Ferret variable.

! We will plot an FFT for the monthly_navy_winds and look for the
! annual cycle. Say "show data" to see the length of the time axis. 

MESSAGE
USE monthly_navy_winds
SHOW DATA

! 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 = uwnd[x=150e:130w@ave,y=20n:40n@ave]
LET FFT_uwndfft = ffta(FFT_uwndtim[l=1:132])
SET VARIABLE/TITLE="Amplitude Spectrum" FFT_uwndfft

SET VIEW ul
PLOT FFT_uwndfft

!  For easier interpretation, invert the frequency axis and plot the 
!  spectrum vs period: months/cycle

!  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 ur
PLOT/VS/LINE/HLIMITS=0:24:2/TITLE="Amplitude Spectrum"/SET_UP per[l=1:`FFT_nf`], FFT_uwndfft
PPL XFOR (I2)
PPL XLAB Period, months/cycle
PPL YLAB 
PPL PLOT


!  Next we will plot the phase
MESSAGE

LET FFT_uwndfftp = fftp(FFT_uwndtim[l=1:132])
SET VARIABLE/TITLE="FFT Phase"/UNITS="deg" FFT_uwndfftp

SET VIEW ll
PLOT FFT_uwndfftp

!  And the phase in months/cycle

SET VIEW lr
PLOT/VS/LINE/HLIMITS=0:24:2/TITLE="FFT Phase"/SET_UP per[l=1:`FFT_nf`],FFT_uwndfftp
PPL XFOR (I2)
PPL XLAB Period, months/cycle
PPL YLAB Deg
PPL PLOT

! clean up
! restore plot state 
CANCEL SYMBOL FFT_*
CANCEL VARIABLE FFT_*
SET MODE/LAST VERIFY	! restore echo behavior