Re: low-pass filter problem

Dave's rules should be adhered to by all when developing filters.

See: http://www.ncl.ucar.edu/Applications/filter.shtml
Example 2 is appropriate to your original question.

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!
>
> Please see the attached filter response curve for 33 weights. You can
> clearly see that at the annual frequency, f=0.83, the 33 weight filter
> transmits about 8% of the original signal. 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.
>
> 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
>
> ------------------------------------------------------------------------
>
>
> ------------------------------------------------------------------------
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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 Sat Jun 06 2009 - 08:52:41 MDT