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