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