The 33 hour low-pass filter
A 33 hour low-pass filter is applied to the hourly time series. The plot below was created with these Ferret commands (33hrfilt is a local file containing the filter weights, the weight values are shown in the output above)
DEFINE AXIS/T=-33:33:1/UNITS=hour TAX1 DEFINE GRID/t=TAX1 g2 FILE/VAR=WEIGHTS33T/GRID=G2 "./33hrfilt" SET VAR/TITLE="33 Hour Filter Weights" WEIGHTS33T PLOT WEIGHTS33T FRAME/FILE=weights33.gif
The plot below shows an example of the high diurnal and tidal frequency energies contained in mooring time series data. The red line is the fourier transform of the above filter. Its sharp cutoff at 33 hours helps preserve low frequency energies and virtually eliminate the high frequency energies which could be aliased to lower frequencies if we did not apply this filter before sub-sampling to one day intervals. Below are the Ferret commands which produced this plot
USE "http://dods.shore.mbari.org/cgi-bin/nph-nc/data/OASISdata/netcdf/hourlyM4.nc" ! Define the time series at a point in space. ! Set the FFT, using explicit time specification. ! Plot the amplitude spectrum vs frequency. LET FFT_uwndtmp = U_UNCORRECTED_HR[Z=30] ! 30 meters depth LET FFT_uwndtim = MISSING(FFT_uwndtmp, 0.0) ! ffta needs data, assign to 0 LET FFT_uwndfft = ffta(FFT_uwndtim) SET VARIABLE/TITLE="U Velocity Energy Spectrum"/units=Amplitude FFT_uwndfft ! For easier interpretation, invert the frequency axis and plot the ! spectrum vs period: hours/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 72 hours where the most energy is. ! The PPL ccommands clean up the plot appearance. SET WIN/SIZE=0.5 PLOT/VS/LINE/HLIMITS=0:72:3/VLIMITS=0:2/TITLE="Amplitude Spectrum"/SET_UP per[l=1:`FFT_nf`], FFT_uwndfft PPL XFOR (I2) PPL XLAB Period, hours/cycle PPL YLAB Energy of M4 ADCP U Velocity Component at 30m go unlabel 4 go unlabel 5 go unlabel 6 go unlabel 8 PPL PLOT ! Plot the filter's transform DEFINE AXIS/T=-33:33:1/UNITS=hour TAX1 DEFINE GRID/t=TAX1 g2 FILE/VAR=WEIGHTS33T/GRID=G2 "./33hrfilt" SET VAR/TITLE="33 Hour Filter Weights" WEIGHTS33T LET FFT_filt = ffta(WEIGHTS33T) LET FFT_filt_nf = `FFT_filt,return=lend` LET FFT_filt_freq1 = FFT_nyquist/ FFT_filt_nf ! Define a frequency axis. DEFINE AXIS/T=`FFT_filt_freq1`:`FFT_nyquist`:`FFT_filt_freq1` filtaxis DEFINE GRID/T=filtaxis gfiltfftfreq LET af = T[g=gfiltfftfreq] ! Define the period from the frequency axis. LET perfilt = 1./af ! Plot the filter over the transform of the data LET W33NORM = FFT_filt*31.43 SET VAR/TITLE="33 Hour Filter Transform" W33NORM PLOT/OVER/LINE/VS perfilt[l=1:`FFT_filt_nf`], W33NORM FRAME/FILE=w33_freq.gif
To use FFTA on data from other moorings which have multiple deployments, and hence gaps in the time axis you need to first regrid the data onto a regular time axis. Here are the Ferret command to do that:
set mem/size=20 def sym tst="`U_UNCORRECTED_HR,return=tstart`" def sym ten="`U_UNCORRECTED_HR,return=tend`" def axis/t=():():1/units=hour treg let u_reg = U_UNCORRECTED_HR[gt=treg] ! Can now pass u_reg to FFTA, or save it to a file...
Last updated: Sat Nov 17 20:31:58 2001 by $Id: hourly2daily.pl,v 1.4 2001/09/20 22:28:15 mccann Exp $