Re: time step when using filwgts_lanczos

From: Dave Allured <dave.allured_at_nyahnyahspammersnyahnyah>
Date: Tue, 28 Oct 2008 18:39:34 -0600

Cecile,

I think example 3 on this page is closest to what you want:

   http://www.ncl.ucar.edu/Applications/filter.shtml

The band pass Lanczos filter is computed and applied in a very few
lines near the top of the program. The two essential function calls
that must be used together for this filter are:

   wgt = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )
   xBPF = wgt_runave ( x, wgt, 0 )

Note: NCL's filter functions operate over discrete time steps in a
time series. Their time metric is time steps, *not* calendar time
or real time. This is a common point of confusion about digital
filters.

Therefore the filters *do not care* about your time units or time
coordinate variable. To get correct filter parameters, you must
manually (or with clever programming) express your time domain
parameters in terms of *time steps*.

 From your example, if the desired filter is 10 to 50 days, and the
time series is on 3-day time steps, then:

   dt = 3 days per time step
   t1 = 50 days (low frequency cutoff, expressed in time domain)
   t2 = 10 days (high frequency cutoff, expressed in time domain)

   fca = dt/t1 = 3./50. = 0.06 time steps (low frequency cutoff)
   fcb = dt/t2 = 3./10. = 0.30 time steps (high frequency cutoff)

The number of filter weights is also important for sharp filter
cutoffs with minimal Gibbs oscillation. My informal rule of thumb
is to use at least 5 times 1/fca. So for this example:

   nwgt >= 5 * 1/fca = 5 * 50./3. = 83 weights

See the reference "Lanczos Filtering in One and Two Dimensions" by
C. Duchon on the filwgts_lanczos NCL page, for more details about
selecting the number of weights.

I recommend that you examine closely the filter response curves
returned by filwgts_lanczos, wgt_at_freq and wgt_at_resp, before using
this filter in a live application.

Dave Allured
CU/CIRES Climate Diagnostics Center (CDC)
http://cires.colorado.edu/science/centers/cdc/
NOAA/ESRL/PSD, Climate Analysis Branch (CAB)
http://www.cdc.noaa.gov/

Cecile Dardel wrote:
> Hi everyone,
>
> I'm trying to filter some satellite data, which have different time
> units and time steps.
> My original NCL script was to band-pass filter daily data, that did have
> a data for each day.
> (it works well).
> Now, I have to deal with weekly data, but the time unit remains being
> the day.
> Another data set I have to filter has a 3-day time step, the unit still
> being the day.
>
> So I don't know how to "tell" the filtering function all these different
> time information.
> I need to filter between 10 and 50 days.
> When I put liminf = 10 and limp = 50, the filter doesn't seem to
> understand I want to band-pass filter between 10 and 50 days.
>
> Do you know how I can do that ?
>
> Hereafter you can find the beginning of my script.
>
> Thank you very much for any help !
>
> Cecile
>
>
> ;=============================================================
> type_filt = 2 liminf = 10 limsup = 50 date1 =
> 20000101
> date2 = 20051231
> ;=============================================================
> if ((data.eq."SSH").and.(expid.eq."TP_ERS"))then
> new_units = "days since 1992-10-14" ; because ut_calendar is
> sensitive to the unit case !!!
> end if
>
> if ((data.eq."SSH").and.(expid.eq."TPJ_JPL"))then
> new_units = "days since 1992-01-01 12:00:00"
> date1 = 20000102
> date2 = 20051231
> end if
>
> print("Reading field")
> r=addfile(indir+infile,"r")
> time = r->time
>
> time_at_units = new_units
> printVarSummary(time)
> date = ut_calendar(time,2) print(date)
> i1 = ind(date.eq.date1)
> i2 = ind(date.eq.date2)
>
> fld1= r->$fldname1$(i1:i2,{south:north},{west:east})
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Oct 28 2008 - 18:39:34 MDT

This archive was generated by hypermail 2.2.0 : Wed Oct 29 2008 - 09:22:41 MDT