Re: low-pass filter problem

Guangshan,

See answers below to your additional questions.

guangshan chen wrote:
> Hi Dave,
>
> Thanks for the help and the scripts.
> Please see my question below on response curve
>
> Thank you very much!
>
> Guangshan
>
>
> On Jun 5, 2009, at 4:59 PM, Dave Allured wrote:
>
>> Guangshan,
>>
>> The number of filter weights is 33. That is not adequate for a good
>> quality Lanczos digital filter with fca=0.06.
>>
>> 1. Informal rule: Always use a number of weights that is at least 5
>> times 1/fca. By this rule you should have at least 5 * 1/0.06 = 83
>> weights.
>>
>> 2. Always check the filter response curves when making digital filters!
> I know more the number or weights more better. I will remember these two
> rules.
>
>>
>> Please see the attached filter response curve for 33 weights. You can
>> clearly see that at the annual frequency, f=0.83, the 33
> the annual frequency, f= 0.083.
>
>> weight filter transmits about 8% of the original signal.

> How to tell 8% of the original signal transmitted? Because I know little
> about the response curves, could you tell me how to check it? what is
> the curve meaning?

f=0.83 was a mistake in my original reply. I should have said only
f=0.083.

The response curve is the amount of signal transmitted by the filter
at different frequencies. Zero = no signal, 1.0 = 100 percent
signal. The filter can be represented as a multiplying function
a(f) over the frequency domain:

   y(f) = a(f) * x(f)

Where x = input signal, y = output signal, a = filter response, and
f = frequency, which is the independent variable.

Now look at the response curve for your original filter with 33
weights (see below). At f = 0.083 (horizontal axis), the function
value is approximately a(0.83) = 0.08 = 8 percent (vertical axis).
So at 0.083 cycles per month, and frequencies very close to that,
about 8 percent of the original signal gets through the filter. In
your case this is signal that you did not want.

My informal rule

>> This explains the unexpected peak at f=0.83 in your output plot.
>>
>> Now compare that point to the same frequency on the 83 weight filter
>> response curve.

> Why is this filter better?

The 83 weight Lanczos filter is better than 33 because the
transition band between 1 and 0 is much narrower in the frequency
domain. This is called "sharp cutoff". This includes more of the
frequencies that you want, and omits more that you don't want.

The number of weights is a trade-off. More weights means sharper
cutoff, but also more data points lost at the beginning and end of
the time series and adjacent to missing values. More weights also
means more computer time, but for many simple applications the
computer time is not significant.

My informal rule nweights = 5 * 1/fca is a first guess for this
trade-off. It can be adjusted more carefully for specific
applications. Some applications need much sharper cutoff between
two frequencies. Others can tolerate a more gradual cutoff.

It all comes back to the same thing. When making a new filter, it
pays to evaluate the response curve in the context of the specific
application.

--Dave

>> Also I have attached a generic NCL script so that you can easily check
>> response curves for different specifications. It uses command line
>> parameters, so you don't need to change the script. There are usage
>> instructions in the header.
>>
>> 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/psd1/
>>
>> Guangshan Chen wrote:
>>> Hi all,
>>> I would like to use a low-pass filter on one monthly data to remove
>>> the signal with period less than 1 year.
>>> When I check the result, I find there still a signal in high
>>> frequency band.
>>> My script as follow:
>>> print("generates the weights for low pass filter")
>>> ihp = 0
>>> nsigma = 1.
>>> befaft = 17 - 1 ; number of points to use before and after
>>> current point
>>> nwt = 2*befaft + 1 ; total number of weights
>>> fca = 0.06
>>> fcb = -999.
>>> wts_t = filwgts_lancos (nwt, ihp, fca, fcb, nsigma)
>>> ; opt=0 is chosen for wgt_runave when applying the weights to, say, a
>>> time series
>>> opt = 0
>>> ; low pass data
>>> x_pass = wgt_runave_Wrap(x,wts_t,opt)
>>> x_new = x_pass(lat|:,lon|:,time|24:tmax-24)
>>> x_org_ts = x({lat|0},{lon|180},time|:)
>>> printVarSummary(x_org_ts)
>>> x_new_ts = x_new({lat|0},{lon|180},time|:)
>>> printVarSummary(x_new_ts)
>>> ;------------------------------------------------------------------
>>> ; set parameter for spectrum
>>> ;------------------------------------------------------------------
>>> print("calculating spectrum of eof time series....")
>>> d = 1 ;Remove the series mean and least squares linear trend.
>>> sm = 7
>>> pct = 0.10
>>> spec_org = specx_anal(x_org_ts,d,sm,pct) splt_org =
>>> specx_ci(spec_org, 0.1, 0.95) spec_new =
>>> specx_anal(x_new_ts,d,sm,pct) splt_new = specx_ci(spec_new,
>>> 0.1, 0.95) Any suggestions?? Also I have another question:
>>> how to set the precision of the Y axis?
>>> I used these two sources:
>>> * res_spec_at_tmYLAutoPrecision = False*
>>> * res_spec_at_tmYLPrecision = 4*
>>> Anything else??
>>> Thanks in advance.
>>> Guangshan<lanczos.response.0.06-33pt.png><lanczos.response.0.06-83pt.png>
>>>
>> ;-----------------------------------------------------------------------------
>>
>> ;
>> ; Plot filter responses.
>> ;
>> ; 2007-may-24 plot.lanczos:
>> ; Original version, by Dave Allured.
>> ; Lanczos, Fca=1/120, nwgt=61, high pass.
>> ; 2007-may-25 Add command line option to print weights.
>> ; 2007-jun-15 Add plot for time domain response.
>> ; 2007-jun-18 30 day, 119 weights.
>> ; 2008-mar-04 Add ps and pdf command line options.
>> ;
>> ; 2008-mar-18 plot-filter:
>> ; Plot Lanczos filter response, all three filter types.
>> ; Generalize. Add command line options for filter specs.
>> ; Add auto centering for plots.
>> ; 2009-jun-04 Add PNG output type.
>> ; 2009-jun-04 Generalize, specify fca/fcb rather than t1/t2.
>> ; Fix bottom alignment of weights curve.
>> ;
>> ;-----------------------------------------------------------------------------
>>
>> ;
>> ; Command line example:
>> ;
>> ; ncl plot-filter.lanczos.ncl lp=1 fca=.06 nwgt=83
>> ;
>> ; Command line parameters.
>> ;
>> ; hp=1 High pass Filter types, specify only one at a time.
>> ; lp=1 Low pass
>> ; bp=1 Band pass
>> ;
>> ; fca=n.n Low cutoff frequency (time steps, 0.0 < fca < 0.5)
>> ; fcb=n.n High cutoff frequency, for band pass only (0.0 < fca <
>> fcb)
>> ;
>> ; nwgt=nnn Number of filter weights
>> ;
>> ; (default) Output plot to X11 window
>> ; ps=1 Output plot to postscript
>> ; eps=1 Output plot to encapsulated postscript
>> ; pdf=1 Output plot to PDF
>> ; png=1 Output plot to PNG
>> ;
>> ; table=1 Also print table of filter weights
>> ;
>> ;-----------------------------------------------------------------------------
>>
>>
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>
>> begin
>>
>> ; Get command line parameters.
>>
>> formats = (/ "ps", "eps", "pdf", "png", "x11" /) ; output selector
>> table
>> x11 = 1 ; default = X11 window
>> ifmt = ind (isdefined (formats(:))) ; find selector in
>> table
>> output = formats(ifmt(0))
>>
>> ; Get filter specs.
>>
>> filter_table = (/ "lp", "hp", "bp" /)
>> filter_names = (/ "lowpass", "highpass", "bandpass" /)
>> alg_name = "lanczos"
>>
>> flags = isdefined (filter_table(:)) ; look up filter cmd in
>> table
>>
>> if (.not. any (flags)) then
>> print ("*** Abort: No filter type specified (hp=1, lp=1, bp=1).")
>> exit
>> end if
>>
>> if (num (flags) .gt. 1) then
>> print("*** Abort: More than one filter type was specified (hp,
>> lp, bp).")
>> exit
>> end if
>>
>> itype = ind (flags)
>> filter_name = filter_names(itype)
>> print ("Filter type = " + filter_name)
>>
>> if (filter_name .ne. "bandpass") then ; this makes fca, fcb etc.
>> work
>> fcb = fca ; out correctly in all cases
>> end if
>>
>> t1 = 1.0 / fcb ; also t2 will always be highest
>> t2 = 1.0 / fca
>>
>> if (filter_name .eq. "bandpass") then
>> freq_str = fca + "-" + fcb
>> ref_line_freq = (/ fca, fcb /)
>> ref_line_time = (/ t1, t2 /)
>> else
>> freq_str = fca + ""
>> ref_line_freq = fca
>> ref_line_time = t1
>> end if
>>
>> ; Generate filter weights using NCL built-in Lanczos function.
>>
>> nsigma = 1. ; fixed filter param.
>>
>> wgt = filwgts_lancos (nwgt, itype, fca, fcb, nsigma)
>>
>> xlim1 = 0 ; set X limits of frequency domain plot
>> xlim2 = max ((/ fca + fcb, fcb * 1.5 /)) ; try to center
>> response curve
>> xlim2 = min ((/ xlim2, 0.5 /)) ; limit to Nyquist freq. (0.5)
>>
>> xrev1 = 0 ; set X limits of time domain plot
>> xrev2 = max ((/ t1 + t2, t2 * 1.6 /)) ; try to center response
>> curve
>>
>> ; Plot filter response.
>>
>> x = wgt_at_freq
>> y = wgt_at_resp
>>
>> spec_string = freq_str + " " + filter_name +", " + nwgt + " weights"
>>
>> plotfile = alg_name + ".response." + freq_str + "-" + nwgt + "pt"
>>
>> wks = gsn_open_wks (output, plotfile)
>>
>> res = True
>> res_at_tiMainString = "Lanczos filter, Response, " + spec_string
>> res_at_tiXAxisString = "Frequency (cycles per time step)"
>>
>> res_at_tiMainFontHeightF = 0.018
>> res_at_tiXAxisFontHeightF = 0.018
>>
>> res_at_gsnXRefLine = ref_line_freq
>> res_at_gsnXRefLineDashPattern = 1
>>
>> res_at_tmXMajorGrid = True ; implement x grid
>> res_at_tmXMajorGridThicknessF = 1.0 ; 2.0 is default
>> res_at_tmXMajorGridLineDashPattern = 2 ; select short dash lines
>>
>> res_at_tmYMajorGrid = True ; implement y grid
>> res_at_tmYMajorGridThicknessF = 1.0 ; 2.0 is default
>> res_at_tmYMajorGridLineDashPattern = 2 ; select short dash lines
>>
>> res_at_gsnMaximize = True ; make plot fill the page
>>
>> res_at_trXMinF = xlim1 ; constrain plot horizontally
>> res_at_trXMaxF = xlim2
>>
>> res_at_trYMinF = -0.1 ; constrain plot vertically
>> res_at_trYMaxF = 1.1
>>
>> plot = gsn_csm_xy (wks, x, y, res)
>>
>> ; Plot filter weights.
>>
>> half = (nwgt-1) / 2.
>> xi = fspan (-half, half, nwgt)
>>
>> res_at_tiMainString = "Lanczos filter, Weights, " + spec_string
>> res_at_tiXAxisString = "Weight Position"
>>
>> res_at_trYMinF = min (wgt) - 0.05 ; constrain plot vertically
>>
>> delete (res_at_gsnXRefLine)
>> delete (res_at_gsnXRefLineDashPattern)
>> delete (res_at_trXMinF)
>> delete (res_at_trXMaxF)
>>
>> plot = gsn_csm_xy (wks, xi, wgt, res)
>>
>> ; Plot response in time domain.
>>
>> xrev = 1 / x(1:) ; transform plot to time domain
>> yrev = y(1:) ; omit divide by zero at first point
>>
>> res_at_gsnXRefLine = ref_line_time ; draw ref. line at cutoff
>> cycle timet
>> res_at_gsnXRefLineDashPattern = 1
>>
>> res_at_tiMainString = "Lanczos filter, Response, " + spec_string
>> res_at_tiXAxisString = "Cycle Time (inverse frequency)"
>>
>> res_at_trXMinF = xrev1 ; constrain plot horizontally
>> res_at_trXMaxF = xrev2
>>
>> plot = gsn_csm_xy (wks, xrev, yrev, res)
>>
>> ; Print filter weights, if selected.
>>
>> if (isdefined ("table")) then
>> print (wgt+"")
>> end if
>>
>> end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk