Frequency Filtering in IDL Using the Butterworth Function

When IDL released version 6.3 they included a function called BUTTERWORTH. This function generates the Butterworth filter function commonly used in frequency filtering of time series data. I have previously used a Butterworth filter by generating my own array. Since it is almost always simpler to use provided functions, I decided to use this new one. The help guide is not very helpful, however, so I made some notes about it while figuring it out. The following discussion is a brief overview of filtering data using fast Fourier transforms (FFTs) and a Butterworth filter.

My test data set is shown below. This is a time series of 4096 points with a sampling rate of 1.5625 MS/s (1.5625 × 106 samples per second). These parameters are chosen just to more closely match those in my work, but the filtering process is the same regardless. The time series is just over 2.5 ms in length and it features an offset (DC) value in addition to a clear oscillation of approximately 5 kHz. On top of this 5 kHz oscillation there is a much faster oscillation.

Sample data for filtering.

The goal is to remove the higher frequency fluctuations so that a cleaner signal comprised of solely the 5 kHz fluctuation can be studied. This removal of the undesired fluctuations will be achieved by transforming the raw data into its frequency representation by way of the FFT.

Before the FFT can be performed, the offset of the data must be removed. In the plot below, the offset of the data is removed by subtracting a smoothed version of the data from itself. From the command line in IDL this might be achieved in the following manner (“data” is the variable name of the data, “time” is the time axis),

;;plot the original signal (black trace in figure)
IDL> plot, time, data
;;plot the smoothed version of the signal (red trace)
IDL> plot, time, smooth( data, 1201, /edge_truncate )
;;plot the fluctuating component of the signal
IDL> plot, time, data − smooth( data, 1201, /edge_truncate )

where the amount of smoothing (the 1201 points) is dependent on the signal and your value may be entirely different.

The original data, its smoothed version, and the zero mean part (which is achieved by subtracting the smoothed signal from the original).

At this point the fluctuating component of the signal, f(t), still has some offset. The remaining offset can be removed by subtracting the mean value of f(t). Our new working signal is s(t),

IDL> s = f – mean( f )

and the resulting signal, s(t), will not be very different from the fluctuating component, f(t).

The real part of the FFT of s(t) and the corresponding frequency array are shown in the figure below. The IDL Reference Guide has a good example of generating the frequency array for use with their FFT function (see the entry for FFT). Notice that the frequency array goes from zero to the largest positive value, and then from the largest negative value back to near zero. Almost all of the FFT amplitude is located at the edges of the frequency array that contain the lowest values because the test signal does not have much in the way of high frequency noise.

Real part of the FFT and the corresponding frequency array.

The basic idea behind frequency filtering with the FFT is that we can reduce the values of the FFT shown above (the full complex value, however, not just the real part) and then perform the inverse transform to return to a regular time signal. The next figure shows the result from the standard usage of the Butterworth function. This is achieved using,

IDL> filter = butterworth( 4096 )

where 4096 is the number of points in the signal (or just the number of points passed to the FFT if you are working with a subset of a larger signal). The figure below shows that the Butterworth filter is near zero for all frequencies except those at the edges. The edge frequencies are also the lowest absolute values, therefore, this filter function would result in a low pass filter and remove most of the higher frequencies in the signal.

Standard Butterworth filter function returned from IDL compared to the standard frequency array of an FFT.

The Butterworth function has a keyword that serves to make this a high pass filter. The following IDL command results in the filter function shown in the next figure,

IDL> filter = butterworth( 4096, /origin )

and the plot shows that the filter function is near zero for all but the largest absolute value frequencies.

Butterworth filter with the ORIGIN keyword specified to return a high pass result.

If the FFT of the true zero mean signal is given as FFT(s(t)), then filtering data using these standard Butterworth filters may be done with,

IDL> lowpass = fft( FFT(s(t)) * low, 1 )
IDL> highpass = fft( FFT(s(t)) * high, 1 )

where “lowpass” and “highpass” are the filtered time signals. The variables “high” and “low“ are the Butterworth functions (recall that the high-pass filter is achieved with the ORIGIN keyword). An example of this filtering is found under the Butterworth entry in the IDL Reference Guide. The next figure plots the time signals that result from these standard filters. Notice that the low-pass signal is very smooth because most of the high frequency noise has been removed. Some low frequency behavior is still seen in the high-pass signal, but it is of very low amplitude (the high-pass signal has been multiplied by 10 in order for it to show up compared to the low-pass amplitude). At this point it is desirable to learn how to set the Butterworth function to return a filter that cuts out frequencies as set by the user.

Results from applying the standard low-pass and high-pass Butterworth filters.

Setting a Frequency Cutoff

The CUTOFF keyword in the Butterworth function can be used to set the lowest (or highest) frequency that will remain in the filtered signal. This keyword must be set to an index number, however, and not a frequency so it is necessary to determine which index (seen in some of the FFT frequency plots above) corresponds to the target frequency.

In the example below I am generating a low-pass filter function.

;; set the maximum frequency to survive (Hz)
IDL> maxFreq = 5e3
;; determine frequency resolution of FFT
;; acq = data acquisition rate (samples / second)
;; FFTsize = number of data points to be passed into FFT function
IDL> deltaFreq = acq / FFTsize
;; determine index value of cutoff frequency
IDL> cutFreq = fix( maxFreq / deltaFreq )
;; generate Butterworth function with the desired cutoff
IDL> filterFun = butterworth( FFTsize, cutoff=cutFreq )
;; filter data (as shown in IDL Reference Guide)
IDL> low = fft( fft( original , -1 ) ∗ filterFun, 1 )

Finally, the results with a 5 kHz and a 100 kHz cutoff setting are shown in the final figure below. The 100 kHz cutoff result still displays a lot of noise, while the 5 kHz cutoff is considerably cleaner. In this case it was easy to set the filter because I could see the strong ≈ 5 kHz oscillation before trying to do any filtering.

Comparison between original signal and data filtered with a set low-pass cutoff.

The point here is not to discuss filtering in general, but only to see how to apply the relatively new Butterworth function provided in IDL. Now, if anyone can explain to me what the XDIM and other keywords mean in the Butterowrth function I would be happy to hear it (you can leave a comment if you wish).

Leave a Reply