Complex filters and data smoothing

Because the data in the Fourier domain is in frequency order, it is often simpler to create a filter in the Fourier domain than to produce the corresponding convolution function in normal space and then transform it. In fact, determining the optimum filter for data is often best done by examining the power spectrum of the data to be filtered and then designing the filter around it. Brault and White (1971) spend some time on this topic.

At present, Figaro provides only one function that generates complex filters, and this is `cmplxfilt'. This produces filters of the type described by Hunstead (1980) for use in cross-correlation measurements. These are mid-pass filters formed by taking a Gaussian that peaks at zero frequency and drops to a half height at a specified high cut value, and subtracting a narrower Gaussian that rises from zero reaching its half height at a specified low cut value. That is, the functional form of the data is given by

$$f(x) = \exp{\left({{-x^2}\over{2 v^2}}\right)} - \exp{\left({{-x^2}\over{2 u^2}}\right)}$$

where u and v determine the low frequency and high frequency cutoff values. From this definition it is clear that u and v are just the sigma values for the two Gaussians, and that what has been rather loosely termed a `half height' is in fact the point at which the Gaussian has a value of exp(-1/2) ~ 0.6. They are specified for `cmplxfilt' in terms of the Nyquist frequency, so a filter might have, say, u = 0.1 and v = 0.3. The best way to get a feel for these filters is to generate them and plot them superimposed on the power spectrum of the data to be filtered.

If the low cutoff value is specified as zero, `cmplxfilt' generates a low pass filter whose functional form is just

$$f(x) = \exp{\left({{-x^2}\over{2 v^2}}\right)}$$

i.e. a Gaussian that drops from unity at zero frequency, having a half width of v. Note that the cyclic nature of the Fourier domain data means that the filter generated is actually symmetrical about the mid point of the data-this is something to be remembered if you want to try to produce other, more specific, filters by generating them yourself. The imaginary parts of the filters generated by `cmplxfilt' are always zero. This means that `cmplxmult' will have the effect of just multiplying both the real and imaginary parts of the data to be filtered by the real part of the filter (obviously, since [a + i b] * [c + i d] = a c + i b c, if d = 0).

`cmplxfilt' requires a template complex data structure, which it will use as the basis for the filter it produces. Normally, this will be the data to be filtered, although any data of the same dimensions will do. If the template data is n-dimensional, so will be the resulting filter. The low and high cutoff frequencies will be the same in all dimensions. Since they are specified in terms of the Nyquist frequency, this is usually fairly satisfactory.

So, to take an easy example, an elaborate way of smoothing a spectrum, `myspect' say, by Gaussian convolution, would be

   ICL> istat myspect stat_mean=(x) reset accept
   ICL> icsub myspect (x) myspect
   ICL> cosbell myspect 10 bell
   ICL> imult myspect bell myspect
   ICL> r2cmplx myspect cmplxspect
   ICL> fft cmplxspect cmplxspect
   ICL> cmplxfilt cmplxspect 0 0.3 cfilt
   ICL> cmplxmult cmplxspect cfilt cmplxspect
   ICL> bfft cmplxspect cmplxspect
   ICL> cmplx2r cmplxspect myspect
   ICL> idiv myspect bell myspect
   ICL> icadd myspect (x) myspect

where the 0,0.3 parameters for `cmplxfilt' indicate that a low pass filter is to be created. The 0.3 indicates a cut at around a third of the Nyquist frequency for the data, and the actual value is one best determined by comparison with the power spectrum of `myspect', obtained by performing a `splot' of `cmplxspect' just after the FFT is obtained. Being able to do this is really the main justification for indulging in such a complex procedure when `ixsmooth' could do the same job far faster.

Note that the sequence given above is a pretty complete one. The use of `istat' at the start is to determine the mean level of the data in the spectrum and the result of the `icsub' is to reduce the spectrum to a zero mean level. This reduces the constant component of the power spectrum and should produce a better result. (Note that `istat' sets its output parameter `stat_mean' to the mean value it determines.) The data is apodised by the application of a 10% cosine bell prior to the transform. Both the effects of the cosine bell and the subtraction of the mean level are reversed at the end of the sequence.