Re: low-pass filter problem

Hi Dave,

Thanks for your kindness.
Now I am clear.

Bests

Guangshan
On Jun 8, 2009, at 3:44 PM, Dave Allured wrote:

> 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<lanczos.response.0.06-33pt.png><lanczos.response.0.06-83pt.png>
> _______________________________________________
> 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 Mon Jun 08 2009 - 15:17:34 MDT