;+
;*****************************************************************************************
;
; FUNCTION : vector_bandpass.pro
; PURPOSE : This program does a bandpass filter on the input data using IDL's
; built-in FFT.PRO routine. The data is first padded with zeroes
; to ensure the number of elements remains an integer power of 2.
; The user defines the input vector array of data, the sample rate,
; and frequency range(s) before running the program, then tells the
; program whether a low-pass (i.e. only return low frequency
; signals), high-pass, or middle frequency bandpass filter. The
; program eliminates the postitive AND negative frequency bins in
; frequency space to ensure symmetry before performing the inverse
; FFT on the data.
;
; CALLED BY: NA
;
; CALLS:
; power_of_2.pro
;
; REQUIRES: NA
;
; INPUT:
; DAT : [N,3]-Array of magnetic or electric field data
; SR : Scalar defining the sample rate (mHz, Hz, kHz, etc. doesn't
; matter as long as everything is consistent)
; LF : Scalar defining the low frequency cutoff (Default = 0)
; [Note: MUST be same units as SR]
; HF : Scalar defining the high frequency cutoff (Default = Nyquist)
; [Note: MUST have same units as SR AND LF]
;
; EXAMPLES:
; htr_mfi2tplot,DATE=date
; get_data,'WIND_B3_HTR(GSE,nT)',DATA=mag
; magf = mag.Y
; tt = mag.X ; -Unix time (s since 01/01/1970)
; nt = N_ELEMENTS(tx)
; evl = MAX(tt,/NAN) - MIN(tt,/NAN) ; -Event length (s)
; nsps = ((nt - 1L)/evl) ; -Approx Sample Rate (Hz)
; lfmf1 = vector_bandpass(magf,nsps,15d-2,15d-1,/LOWF)
;
; KEYWORDS:
; LOWF : If set, program returns low-pass filtered data with freqs
; below LF
; MIDF : If set, program returns bandpass filtered data with freqs
; between LF and HF **[Default]**
; HIGHF : If set, program returns high-pass filtered data with freqs
; above HF
;
; CHANGED: 1) Fixed Low Freq. bandpass to get rid of artificial
; zero frequency bin created by FFT calc. [01/14/2009 v1.0.1]
; 2) Fixed case where NaN's are in data [01/18/2009 v1.0.2]
; 3) Changed program my_power_of_2.pro to power_of_2.pro
; and renamed [08/10/2009 v2.0.0]
;
; CREATED: 12/30/2009
; CREATED BY: Lynn B. Wilson III
; LAST MODIFIED: 08/10/2009 v2.0.0
; MODIFIED BY: Lynn B. Wilson III
;
;*****************************************************************************************
;-
FUNCTION rbsp_vector_bandpass,dat,srt,lf,hf,LOWF=lowf,MIDF=midf,HIGHF=highf
d = REFORM(dat)
no = N_ELEMENTS(d[*,0])
;-----------------------------------------------------------------------------------------
; -Get rid of any NaN's so FFT can actually work
;-----------------------------------------------------------------------------------------
bbx = WHERE(FINITE(d[*,0]) EQ 0,bx,COMPLEMENT=ggx)
bby = WHERE(FINITE(d[*,1]) EQ 0,by,COMPLEMENT=ggy)
bbz = WHERE(FINITE(d[*,2]) EQ 0,bz,COMPLEMENT=ggz)
IF (bx GT 0L) THEN d[bbx,0] = 0d0
IF (by GT 0L) THEN d[bby,1] = 0d0
IF (bz GT 0L) THEN d[bbz,2] = 0d0
;-----------------------------------------------------------------------------------------
; -Pad data with zeros to change no => 2^m {m = integer}
;-----------------------------------------------------------------------------------------
d2x = power_of_2(d[*,0])
d2y = power_of_2(d[*,1])
d2z = power_of_2(d[*,2])
d2 = [[d2x],[d2y],[d2z]]
nd = N_ELEMENTS(d2x)
n_m = nd/2L + 1L ; -mid point element
;-----------------------------------------------------------------------------------------
; -Calc FFT frequencies
;-----------------------------------------------------------------------------------------
sr = srt*1d0
frn = nd - n_m
frel = LINDGEN(frn) + n_m ; -Elements for negative frequencies
fft_freq = LINDGEN(nd)
fft_freq[frel] = (n_m - nd) + DINDGEN(n_m - 2L)
fft_freq = fft_freq*(sr/nd)
;-----------------------------------------------------------------------------------------
; -Determine relevant elements of FFT arrays
;-----------------------------------------------------------------------------------------
lfc1 = lf*1d0
hfc1 = hf*1d0
lowf1 = WHERE(ABS(fft_freq) LE lfc1 AND ABS(fft_freq) GT 0d0,lf1,COMPLEMENT=other_mh)
midf1 = WHERE(ABS(fft_freq) GT lfc1 AND ABS(fft_freq) LE hfc1,mf1,COMPLEMENT=other_lh)
highf1 = WHERE(ABS(fft_freq) GT hfc1,hf1,COMPLEMENT=other_lm)
IF KEYWORD_SET(lowf) THEN lowf1 = 1 ELSE lowf1 = 0
IF KEYWORD_SET(midf) THEN midf1 = 1 ELSE midf1 = 0
IF KEYWORD_SET(highf) THEN highf1 = 1 ELSE highf1 = 0
check = [KEYWORD_SET(lowf1),KEYWORD_SET(midf1),KEYWORD_SET(highf1)]
gcheck = WHERE(check GT 0,gch,COMPLEMENT=bcheck,NCOMPLEMENT=bch)
gelems = {T0:other_mh,T1:other_lh,T2:other_lm}
IF (gch EQ 1L) THEN BEGIN
other = gelems.(gcheck[0])
ENDIF ELSE BEGIN ; -Default setting b/c user entered too many keywords
other = other_lh
ENDELSE
;-----------------------------------------------------------------------------------------
; -Window data
;-----------------------------------------------------------------------------------------
window = hanning_stretch(n_elements(d2x),n_elements(d2x)/8L-2L)
;-----------------------------------------------------------------------------------------
; -Calc FFT
;-----------------------------------------------------------------------------------------
tempx = FFT(window*d2x,/DOUBLE)
tempy = FFT(window*d2y,/DOUBLE)
tempz = FFT(window*d2z,/DOUBLE)
templx = tempx
temply = tempy
templz = tempz
;-----------------------------------------------------------------------------------------
; -Get rid of unwanted frequencies
;-----------------------------------------------------------------------------------------
templx[other] = DCOMPLEX(0d0) ; -Get rid of unwanted frequencies [mid and high]
temply[other] = DCOMPLEX(0d0)
templz[other] = DCOMPLEX(0d0)
;-----------------------------------------------------------------------------------------
; -Calc Inverse FFT
;-----------------------------------------------------------------------------------------
rplx = REAL_PART(FFT(templx,1,/DOUBLE))
rply = REAL_PART(FFT(temply,1,/DOUBLE))
rplz = REAL_PART(FFT(templz,1,/DOUBLE))
;-----------------------------------------------------------------------------------------
; -Keep only useful data [i.e. get rid of the zero-padded elements]
;-----------------------------------------------------------------------------------------
filtered = [[rplx[0:(no-1L)]],[rply[0:(no-1L)]],[rplz[0:(no-1L)]]]
IF (bx GT 0L) THEN filtered[bbx,0] = !VALUES.D_NAN
IF (by GT 0L) THEN filtered[bby,1] = !VALUES.D_NAN
IF (bz GT 0L) THEN filtered[bbz,2] = !VALUES.D_NAN
RETURN,filtered
END