Viewing contents of file '../idllib/contrib/harris/filter.pro'
;-------------------------------------------------------------------
function filter,in,pl,ph,freq=freq,time=time,notch=notch,$
lowpass=lowpass,highpass=highpass,bandpass=bandpass,$
type=type, wts=freq_filter_wts
;+
; NAME: filter
;
; PURPOSE: perform digital filtering in the frequency domain
; via FFT
;
; CATEGORY: Time-series/Spectral analysis
;
; CALLING SEQUENCE:
; result = filter(time_series, /TIME, $
; lowperiod_cutoff, highperiod_cutoff, /BANDPASS)
;
; result = filter(time_series, /FREQ, $
; lowfreq_cutoff, highfreq_cutoff, /BANDPASS)
;
; result = filter(time_series, /FREQ, $
; lowfreq_cutoff, highfreq_cutoff, /NOTCH)
;
; result = filter(time_series, maxperiod, /LOWPASS)
;
; result = filter(time_series, $
; highfreq_cutoff, /FREQ, /LOWPASS)
;
; result = filter(time_series, $
; lowfreq_cutoff, /FREQ, /HIGHPASS)
;
; result = filter(time_series, /FREQ, WTS = wts)
;
; INPUTS: time_series = input time series to filter
; low_cutoff, high_cutoff = the limits (either in
; frequency or time, dependent on the keyword usage) for
; the filter. Only one paramter is required if the filter
; is one-sided (LOWPASS or HIGHPASS). These paramters
; are ignored if the WTS keyword is used
; KEYWORD PARAMETERS:
; FREQ, TIME = signals how the cutoff parameters will be
; interpretted. Either as PERIODS in the time domain or
; FREQUENCIES in the frequency domain. Default is /TIME
; These paramters are ignored if the WTS keyword is used
; LOWPASS = low-pass filter. Use the first cutoff parameter as
; the step function cutoff value
; HIGHPASS = high-pass filter. Use the first cutoff parameter as
; the step function cutoff value
; BANDPASS = band-pass filter. Use both cutoff parameters as
; the limits of the rectangular pass window.
; NOTCH = notch filter. Use both cutoff parameters as
; the limits of the rectangular rejection window.
; WTS = independently set the window weights in the
; frequency domain. This will override all other
; parameters.
; TYPE = can have 1 of four values (case insensitive).
; "LOWPASS", "HIGHPASS", "NOTCH", "FREQ_FILTER_WTS"
; OUTPUTS:
; result = filtered time series
;
; COMMON BLOCKS:
; none.
; SIDE EFFECTS:
; none.
; MODIFICATION HISTORY:
; Written by: R.A.Vincent, Physics Dept., University of Adelaide,
; Early 1989.
; Enhanced: T.J.Harris, Physics Dept., University of Adelaide,
; Late 1989.
;
;-
;in = input time series (may be complex)
;pl and ph are the appropriate cutoff periods in time units
; or freq units if keyword 'freq' set.
;pl,ph included in the pass band
;freq_filter_wts will override all other parameters and be used as
;the filter weights
if (keyword_set(lowpass)) then type = 'lowpass'
if (keyword_set(highpass)) then type = 'highpass'
if (keyword_set(bandpass)) then type = 'bandpass'
if (keyword_set(notch)) then type = 'notch'
if (keyword_set(freq_filter_wts)) then begin
type = 'freq_filter_wts'
pl = 0 & ph = 0
endif
sz = size(in)
nmax = sz(sz(0)+2)
datatype = sz(sz(0)+1)
ntop = nmax-1
temp = complexarr(nmax)
if (datatype ne 6) then temp = complex(in) else temp = in
temp = fft(temp,-1)
if (n_elements(ph) le 0) then ph = pl
temp2 = ph > pl
pl = pl < ph
ph = temp2
if (keyword_set(freq)) then begin
;the bandpass limits are in frequency units 0 --> ntop/2
npl = (fix(pl) > 0) < ntop
nph = (fix(ph) > 0) < ntop
endif else begin
;the bandpass limits are in time units 0 --> nmax
npl = (fix(float(nmax)/ph) > 0) < ntop
nph = (fix(float(nmax)/pl) > 0) < ntop
endelse
temp2 = nph > npl
npl = npl < nph
nph = temp2
CASE (strupcase(type)) OF
'LOWPASS': begin
if (npl lt ntop) then temp(npl+1:ntop-npl) = 0.0
end
'HIGHPASS': begin
if (npl gt 0) then temp(0:npl-1) = 0.0
if (npl gt 2) then temp(ntop-npl+2:ntop) = 0.0
end
'BANDPASS': begin
;if (npl gt 0) then temp(0:npl-1) =0.0
;if (nph lt ntop) then temp(nph+1:ntop-nph) =0.0
;if (npl gt 2) then temp(ntop-npl+2:ntop) =0.0
temp2 = temp
temp = temp2 * 0.0
temp(npl:nph) = temp2(npl:nph)
temp(ntop-nph+1:ntop-npl+1)=temp2(ntop-nph+1:ntop-npl+1)
end
'NOTCH': begin
temp(npl:nph) =0.0
temp(ntop-nph+1:ntop-npl+1)=0.0
end
'FREQ_FILTER_WTS': temp = temp*freq_filter_wts
ENDCASE
return,fft(temp,1)
end