NCL Home>
Application examples>
Data Analysis ||
Data files for some examples
Example pages containing:
tips |
resources |
functions/procedures
NCL: Filters
Lanczos Filter Weights
Filters require that a set of weights be applied to data.
The weights may be applied in the spatial
(eg, smth9) or time domains. The focus
of the following examples will be on application to the temporal domain.
The filwgts_lanczos function may be used to
create a set of weights that have characteristics specified
by the user. The wgt_runave or
wgt_runave_Wrap functions
can be used to apply the weights. The use of wgt_runave
or wgt_runave_Wrap
will result in the loss of data at the beginning and end of
the series being filtered. There is no way to avoid this data loss.
Each Lanczos filter will have an associated response function.
This shows how weights which are applied in the (say)
time domain will affect amplitudes in the frequency domain.
Often, the phrase half-power point is used to describe the
overall effectivesness of a filter. This is the
frequency (period) at which the power or variance
is 50% of maximum. This is the frequency (period)
at which the amplitude response is 0.707 [0.707^2=0.5].
Note: Be aware that the frequency at which the
amplitude response is 0.5 is sometimes referred to as
the half-power point. This is not correct.
In response to a question posted to ncl-talk.ucar.edu, Dave Allured (CU/CIRES),
posted the following:
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*.
For 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)
See Example 7.
Butterworth Bandpass Filter
The Butterworth filter is a signal processing filter designed to have
as flat a frequency response as possible in the passband. NCL has a
function bw_bandpass_filter which is optimized for
narrow band applications.
Reference:
Electronic Supplement to Development of a Time-Domain, Variable-Period
Surface Wave Magnitude Procedure for Application at Regional and
Teleseismic Distances, Part I: Theory; David R. Russell
Bulletin of the Seismological Society of America
http://www.seismosoc.org/publications/BSSA_html/bssa_96-2/05055-esupp/
Lanczos and Butterworth Response Functions
Figures which show the response functions of Lanczos or
Butterworth (V6.3.0) filters also show the ideal filter response
in the frequency domain. For band-pass filters, this ideal response
is sometimes referred to as the 'boxcar' response.
Lanczos Response Function:
The filters_8.ncl
script will generate one or more response functions for different
Lanczos filter lengths. See Example 8.
Butterworth Response Function (6.3.0):
The filters_9.ncl
script will generate one or more response functions for different
Butterworth pole values.
filter_1.ncl: Reads in a variable
from a netCDF file, averages all latitudes between 15S and 15N,
applies a user specified filter on the "time" dimension,
and creates a color Hovmueller diagram.
gsn_csm_hov is the graphical
interface that creates Hovmueller
diagrams.
filters_2.ncl:
Low pass filters via filwgts_lanczos:
- Read the monthly Southern Oscillation Index (Darwin only)
- Using different filter lengths (49 [solid], 121 [dashed]),
generate two low pass filters that (ideally) will remove fluctuations
less than 24 months (2 years).
- Using different filter lengths (121 [solid], 241 [dashed]),
generate two low pass filters that (ideally) will remove fluctuations
less than 120 months (10 years).
The blue (49 weights) and red (121 weights) curves are for the
24-month filter. In this case, not much was gained by using much
longer filter (121 vs 49). Also, using the 49-point filter returned
more data at the ends.
The green (121 weights) and black (241 weights) curves are for the
120-month filter. The response functions are not much different and,
as a result, the decadal curves are not much different.
filters_3.ncl:
Band pass filters via filwgts_lanczos:
Generate and apply a 20-100 day band pass filter using 201 weights.
The response function for the filter is also shown.
The Lanczos filter weights are applied to the entire
daily series spanning 1980-2005.
The plotted values do not have any missing points
'end-points' because the plotted time period (1995-1998)
is a subset of the entire series.
filters_4.ncl:
Comparison: band pass filters via filwgts_lanczos and via FFT:
This example compares results returned using
filwgts_lanczos and FFTs [ezfftf and ezfftb].
The band pass period for this example is 30-to-60 days.
The FFT filtering is performed by applying a 'boxcar' cutoff
in frequency space (1/30 and 1/60) while the Lanczos
weights are applied in time. The Lanczos band pass filter in this example is 731 points.
This means that exactly one year of data will be
'lost' at the beginning and end of the series to which the Lanczos filter is applied.
(a) Read entire 1980-2005 daily anomaly data [here, anomolous zonal wind]
for a specific location [0,100E]
(b) Apply the Lanczos filter weights to the entire 1980-2005 series. The loss
of data at the end points will affect 1980 and 2005. [BPF, black line]
(c) Select a four year [1995-1998] temporal window to demonstrate
how using a segment of an infinitely long [well, ok, 1980-2005 :-) ]
would compare.
(d) Apply the Lanczos filter to the 4-year temporal window only. The filtered
series for 1995 and 1998 will be set to missing. [bpf, red line]
(e) Perform a Fourier analysis (ezfftf) on the input anomalies for the
4-year temporal window; set the Fourier coefficents outside the
desired frequency band to zero (boxcar); perform a Fourier synthesis (ezfftb)
to return to temporal window space. [FFT, blue line]
For simplicity, no tapering (taper) or detrending (dtrend) was performed.
If the
BPF line is considered "truth", then,
in this case (a fairly broad band pass), the
FFT approach yields similar results to the
BPF over the temporal window subperiod. There are some minor phase differences
outside at the 'ends' and one issue in 1997 but, overall, pretty good agreement.
In addition, there is apparently no loss of data at the beginning and end of the
sub-period. Why not use the FFT approach all the time?
If the number of weights had been kept the same, but the band pass period
was much shorter, then the fourier synthesis would likely show more
distinct differences.
filters_5.ncl:
Band pass filters via FFTs: End point problems:
This demonstrates one 'problem' with using boxcar type truncation
of Fourier coefficients in the frequency domain. Consider a series of
length N (here, N=1000 days) which will have N/2 frequencies (0 to 0.5 cycles/day).
In this example, the 20-100 day band pass period suggested by MJO Clivar will be used.
(a) Create Fourier coefficients (real and imaginary) which are 1.0 in the
frequency band of interest and 0.0 elsewhere. This is the boxcar response in
frequency space (top figure). In temporal space, it would take an infinite
number of weights to yield the boxcar frequency response.
(b) Perform a Fourier synthesis (ezfftb) on (a) to return to temporal
space (middle figure).
(c) Generate Lanczos weighting using 201 weights; 'map' the Lanczos response function into FFT
frequency space via linear interpolation; apply the mapped weights to the
boxcar Fourier coefficients in the frequency domain; perform a Fourier
synthesis (bottom figure).
The middle figure shows that a Fourier synthesis of the boxcar Fourier coefficients yields
spurious and significant 'end-effects' and, in this idealized case, minor ripples at
all times. In the real world, the real and imaginary coefficients
are not all equal [here, 1.0], the synthesis can have additional unwanted effects.
The reconstructed series can reflect local and, possibly, significant (in time) occurring features.
Is there a way to use an FFT and get the desirable and predictable characteristics
Lanczos weighting? Weighting the boxcar Fourier coefficients by the Lanczos responses
yields slightly smaller ripples and end-effects (bottom figure).
Punch line: The Fourier approach yields values but you are not
quite sure of what you are getting. At least with the
Lanczos filter approach, you know what you are getting.
filters_6.ncl:
FFT speed and Lanczos clarity:
(a) Develop a Lanczos filter that has good characteristics
(b) Map the Lanczos response function into FFT frequency space
via linear interpolation
(c) Apply the mapped weights to the Fourier coefficients.
(d) Perform the Fourier synthesis.
(e) In the time domain, manually set the end-points to missing.
There is no way around this!!!!
filters_7.ncl:
Matching data with different sampling intervals.
A user had daily data which was treated with a 10-50 day
band pass filter. A different series was availble but
values were available every 3rd day. How to create a filter
that will operate on the latter series which will 'match' that used
on the daily data. (See the above comments by Dave Allured.)
The figure shows a comparison. Do they match exactly? Of course not!
(a) Information was lost by sampling every 3rd day. There is no way
to recover this lost information; (b) the sereies are not plotted
at the exact same time points.
Does the new filtered series capture the essentials of the
fully samples series? Yes!
filters_8.ncl:
Lanczos Response Functions
This script will generate one or more response functions for different
Lanczos filter lengths. This allows users to assess which filter length
is acceptable for their application. There is a trade-off
amongst bandwidth, filter length and response.
This example was run twice to get the 2 figures.
filters_9.ncl:
Butterworth Response Functions (V6.3.0)
This script will generate one or more response functions for different
Butterworth pole values. The left figure shows the bandpass relative to the full
frequency span. The right figure shows a blowup of the region of interest.
bfband_1.ncl:
This script reads a time series at a user specified grid point. It applies
a Butterworth band Pass filter (
bw_bandpass_filter) to the series.
The raw (top) and filtered/envelope (bottom) series illustrates usage and results.
bfband_2.ncl:
Use Lanczos filter (
filwgts_lanczos) with 1461 weights (=4*365+1) to filter a
long time series. Basically 2 years of daily data are lost at
each end of the series.
Also, apply a Butterworth filter (
bw_bandpass_filter). No data are lost.
The rightmost figure shows the beginning of the full time series.
bfband_3.ncl:
Plot spatial patterns of 40-50 day Butterworth band passed data every 5 days.