Date: Fri, 05 Jun 2009 15:59:57 -0600
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
;-----------------------------------------------------------------------------
;
; 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
