function lowpass_3d, data, period

; Purpose
; -------
; Uses the FFT function to implement a low-pass filter on a three-dimensional
; dataset
;
; Inputs
; ------
; data		raw data
; period	minimum period, below which to filter out variability
;
; Outputs
; -------
; lowpass	filtered data
;
; History
; -------
; 2005 Aug 22	Steve Phipps	Original version
; 2005 Aug 23	Steve Phipps	Modified to use the averages of the first and
;				last five values when calculating the linear
;				trend within the data. This avoids the ringing
;				that arises when either the first and/or last
;				values are outliers.
; 2006 May 11	Steve Phipps	Modified to calculate the linear trend by
;				calculating the line of best fit

; Derive the dimensions of the data array
nx = n_elements(data(*, 0, 0))
ny = n_elements(data(0, *, 0))
nt = n_elements(data(0, 0, *))

; Declare local variables
trend = fltarr(nx, ny, nt)

; Remove any linear trend from the data
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    ab = linfit(findgen(nt), data(i, j, *))
    trend(i, j, *) = ab(0) + ab(1) * findgen(nt)
  endfor
endfor
data = data - trend

; Calculate the Fourier transform of the raw data
data_fft = fft(data, dimension=3, /double)

; Apply a low-pass filter, by setting the Fourier transform equal to zero for
; periods shorter than the specified minimum, and then calculating the inverse
; Fourier transform.
;
; data_fft(*, *, 0) contains the data for zero frequency, while the other
; elements contain the data for the following frequencies:
;
; data_fft(*, *, 1)            ->   f = 1/nt
; data_fft(*, *, 2)            ->   f = 2/nt
; ...                               ...
; data_fft(*, *, nt/period)    ->   f = 1/period     MINIMUM PERIOD
; ...                               ...
; data_fft(*, *, nt/2)         ->   f = 1/2          NYQUIST FREQUENCY
; ...                               ...
; data_fft(*, *, nt-nt/period) ->   f = -1/period    MINIMUM PERIOD
; ...                               ...
; data_fft(*, *, nt-2)         ->   f = -2/nt
; data_fft(*, *, nt-1)         ->   f = -1/nt
max_freq = nt / period
data_fft(*, *, max_freq+1:nt-max_freq-1) = dcomplex(0.0)
lowpass = fft(data_fft, dimension=3, /double, /inverse)
lowpass = float(lowpass)

; Add any linear trend back onto the data
lowpass = lowpass + trend

return, lowpass

end