Re: [ferret_users] frequency_response_curve

[Edward Cokelet - NOAA Federal <edward.d.cokelet@xxxxxxxx>]

Chinnu, I attach a .jnl file that computes and plots the frequency response function of the lsl_lowpass filter.  The idea is to define a Dirac delta function centered in the center of a time axis (1024 hours long in this example).  Then compute the least squares Lanczos lowpass filter of this function over the time interval at a fixed cutoff period (35 hours in this example) with a series of different filter span widths (ranging from 3 to 257 hours in this example).  The delta function and the sequence of low-passed filtered values are plotted in the time domain in window 1.  Then compute the fast fourier transform amplitude of each of the sequence of low pass filtered values and plot those in the frequency domain in window 2.

P.S. When I programmed this in 2011, there was what I determined to be a ferret bug that claimed that the index l (“el”) had missing values.  I got around that by writing out the values of the sequence of low pass filtered values as netcdf files and reading them back in.  That code is in this example.  I have not tried a later ferret version without that step.  I interpreted it as a ferret bug a the time, but it might be a ferret “feature” as interpreted by other people.    

If you have a problem with this, e-mail me privately, and we can try to sort it out.  FYI, I am not a ferret developer, but just a user like yourself.  
Here’s the ferret code:  




! Computes and plots frequency response function for the Least Squares Lanczos filter

! Programmed by E.D. Cokelet, NOAA/PMEL, Seattle
! Last modified 5 April 2011

! Define delta function

define axis/t=1:1024:1/unit=hour t_ax			! An hourly time axis that is 1024 hours long
let delta = if (l eq 512) then 1 else 0*t[gt=t_ax]	! Define a delta function at the central time

let cutoff_period = 35					! Low frequency cutoff (=half-amplitude point) at 35 hours

let lp3  = lsl_lowpass(delta[l=1:1024], cutoff_period,  3)
let lp5  = lsl_lowpass(delta[l=1:1024], cutoff_period,  5)
let lp9  = lsl_lowpass(delta[l=1:1024], cutoff_period,  9)
let lp17 = lsl_lowpass(delta[l=1:1024], cutoff_period, 17)
let lp33 = lsl_lowpass(delta[l=1:1024], cutoff_period, 33)
let lp65 = lsl_lowpass(delta[l=1:1024], cutoff_period, 65)
let lp129 = lsl_lowpass(delta[l=1:1024], cutoff_period, 129)
let lp257 = lsl_lowpass(delta[l=1:1024], cutoff_period, 257)

! Write out results to get around a ferret bug that keeps resetting the range of l and then claiming missing values
list/format=cdf/clobber/file=lp3.nc lp3
list/format=cdf/clobber/file=lp5.nc lp5
list/format=cdf/clobber/file=lp9.nc lp9
list/format=cdf/clobber/file=lp17.nc lp17
list/format=cdf/clobber/file=lp33.nc lp33
list/format=cdf/clobber/file=lp65.nc lp65
list/format=cdf/clobber/file=lp129.nc lp129
list/format=cdf/clobber/file=lp257.nc lp257

set window 1
plot delta,lp3,lp5,lp17,lp33,lp65,lp129,lp257
frame/format=gif/file=freq_response_delta_function_lsl_lowpass.gif

! Cancel variables and read back in from files to get around the ferret bug
can var lp3,lp5,lp9,lp17,lp33,lp65,lp129,lp257

set window 2

use lp3
let fft3  = ffta( lp3[d=lp3,l=2:1023] )
plot/vlimits=0:1.05  fft3/fft3[l=@max]

use lp5
let fft5  = ffta( lp5[d=lp5,l=3:1022] )
plot/over/line fft5/fft5[l=@max]

use lp9
let fft9  = ffta( lp9[d=lp9,l=5:1020] )
plot/over/line fft9/fft9[l=@max]

use lp17
let fft17  = ffta( lp17[d=lp17,l=9:1016] )
plot/over/line fft17/fft17[l=@max]

use lp33
let fft33  = ffta( lp33[d=lp33,l=17:1008] )
plot/over/line fft33/fft33[l=@max]

use lp65
let fft65  = ffta( lp65[d=lp65,l=33:992] )
plot/over/line fft65/fft65[l=@max]

use lp129
let fft129  = ffta( lp129[d=lp129,l=65:960] )
plot/over/line fft129/fft129[l=@max]

use lp257
let fft257  = ffta( lp257[d=lp257,l=129:896] )
plot/over/line fft257/fft257[l=@max]

plot/over/vs/line/dash/color=black/thick=2/nolabel {`1/cutoff_period`, `1/cutoff_period`},{0., 1.0}

frame/format=gif/file=freq_response_fft_lsl_lowpass.gif





------------------------------------------------------------------------------------------
Edward D. (Ned) Cokelet, Ph.D. Oceanographer
NOAA/PMEL                         off:  (206) 526-6820
7600 Sand Point Way NE     fax: (206) 526-6485
Seattle, WA 98115-6349

Edward.D.Cokelet@xxxxxxxx

The contents of this message are mine personally and do not necessarily reflect any position of the Government or the National Oceanic and Atmospheric Administration.


BEGIN:VCARD
VERSION:3.0
N:Cokelet;E. D. (Ned);;;
FN:E. D. (Ned) Cokelet
ORG:NOAA/PMEL;
TITLE:Oceanographer\, PhD
EMAIL;type=INTERNET;type=WORK:Edward.D.Cokelet@xxxxxxxx
TEL;type=WORK;type=pref:(206) 526-6820
TEL;type=WORK;type=FAX:(206) 526-6485
item1.ADR;type=WORK:;;7600 Sand Point Way N.E.;Seattle;WA;98115-6439;
item1.X-ABADR:us
NOTE:The contents of this message are mine personally and do not necessarily reflect any position of the Government or the National Oceanic and Atmospheric Administration.
CATEGORIES:PMEL Cyclists
X-ABUID:B7E1C068-924B-11D8-B26B-000A9595D2D2\:ABPerson
END:VCARD