Re: time step when using filwgts_lanczos

Dave, Nice explanation. THX.

I took a bandpass example [4] from

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

This had values every day. The daily values were then sampled
every 3rd day [decimated or subsampled ] to simulate the
problem described by Cecile.

; ***********************************************
; create the filter weights and apply
; ***********************************************
   ihp = 2 ; band pass
   sigma = 1.0 ; Lanczos sigma

   nWgt = 201 ; loose 100 each end
   fca = 1./50. ; start freq
   fcb = 1./10. ; last freq
   wgt = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )
   xbpf = wgt_runave_Wrap ( X, wgt, 0 ) ; entire series

   nWgt3 = nWgt/3 ; 83 [DA] ; loose mWgt/2 each end
   fca = 3./50. ; start freq
   fcb = 3./10. ; last freq
   WGT3 = filwgts_lanczos (nWgt3, ihp, fca, fcb, sigma )
   XBPF3 = wgt_runave_Wrap ( X(::3), WGT3, 0 ) ; decimate entire series

A sample of the daily filtered series and the filtered
decimated series is at:

      http://www.cgd.ucar.edu/~shea/filter.decimate3.png

A perfect match? ... of course not ... for several reasons:

[a] The decimated series is has 'cruder' spacing.
[b] Information was lost by the subsampling. There is *no way*
to recover this lost information.

Useful? YES!

I would say that the filtered subsampled data captures
all the essential features of the filtered daily
sampled series.

D

Dave Allured wrote:
> 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
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Oct 29 2008 - 09:10:09 MDT