How to use Ferret External Functions for computing FFT's

Below is an annotated version of the script ef_fft_demo.jnl

\cancel mode verify
! ef_fft_demo.jnl
! *am* 4/99
! Description: Demonstration of external functions for Fast Fourier Transforms:

! FFTA computes amplitude spectrum
! FFTP computes phase.

! NOTE: The argument must have an 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.
! - Set the units for the frequency axis to "1./(time units)" with the time units taken from the Ferret variable.

! Example:

! 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.

```

yes? USE monthly_navy_winds
yes? SHOW DATA
currently SET data sets:
1> /opt/local/ferret/fer_dsets/descr/monthly_navy_winds.des  (default)
name     title              I         J         K         L
UWND     ZONAL WIND        1:144     1:73      ...       1:132
VWND     MERIDIONAL WIND   1:144     1:73      ...       1:132
```

! Define the time series averaged over a region in space.
! Set the FFT, using an explicit time specification[l=1:132]
! 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 VIEW ul
SET VARIABLE/TITLE="Amplitude Spectrum" FFT_uwndfft
PLOT FFT_uwndff```

! 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 period vs FFT Amplitudes showing the first 24 months where most of the energy is.
! The PPL commands clean up the appearance of the plot.

``` SET VIEW ur
PLOT/VS/LINE/HLIMITS=0:30: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```

! Compute and plot the phase using fftp.

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

SET VIEW ll
PLOT FFT_uwndfftp

SET VIEW lr
PLOT/VS/LINE/HLIMITS=0:30: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```