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


  


oar.pmel.contact_ferret@noaa.gov 30 Apr, 1999 404 Not Found

Not Found

The requested URL "/footer.txt" was not found on this server.