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 × 10⁶ 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.

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),

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

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),

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.

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,

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.

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,

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

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,

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.

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.

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.

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).