Downsample: resample vs antialias fitlering + decimation

Question

I have a discrete signal sampled @Fs. I need to downsample it to Fs/k. Main reason is to reduce signal bandwidth and speed up computation : I'm only interested in a reduced bandwidth < Fs/2k. Some high frequency noise >> Fs/2k can exist : antialiasing is needed. From what I understand, I have several options :

1. brute force : decimation (without taking care of aliasing problem!)
2. antialiasing + decimation : I can use for example numpy decimate function
3. frequency approach : FFT then truncation then IFFT, can be done using resample

Only solution 2- 3- help solving antialiasing issue. 3- filters in frequency domain, supposing signal is periodic. 2- filters with IIR/FIR filter.

I can always use decimate without questioning, but I would like to precisely understand limitations behind. What are drawbacks / advantages of using solution 2- versus 3-? Do you have some practical guidelines about downsampling? Some recommendations?

Answer 1

I will explain why method 2 is often a better choice over method 3.

The frequency domain approach is equivalent to the "Windowing" method of filter design- in that to do that approach correctly you should window your data before taking the FFT. For an anti-alias filter design in the time domain approach, the least squares filter design algorithm outperforms window design approaches. (See this post for a detailed discussion on that: FIR Filter Design: Window vs Parks McClellan and Least Squares). For time-domain filters for decimation and interpolation applications, the least-squares filter design is a better choice over equiripple because of the stop-band roll-off: for equiripple the stop-band is at the same level in each aliased frequency band resulting in more overall noise folding in than you would get with least-squares.

Side Note: If you observe the coefficients for an equiripple design, you will often observe (if the filter isn’t too long) two slightly larger "impulses" toward the beginning and end of the filter impulse response (the coefficients of the filter is the impulse response). Remove those larger coefficients near the tails of the response and the equiripple design will also have the desired feature of stop-band roll-off! For further details on that please see: Convert a Park McClellan FIR Solution to Achieve Stop-band Roll-off

Further the least-squares (and equiripple) design tools in Matlab/Octave/Python feature multi-band filter design which is ideal for decimation (and interpolation) applications since the images are limited to distinct bands. Thus you can optimize the filter rejection to just the frequency locations that would fold in, further optimizing the solution given the same number of taps. Below is an example spectrum I have recently shown for interpolation and the resulting multi-band filter designs for both least squares and equiripple appropriate for eliminating the images (this is the interpolation filter to grow the zeros that are inserted to their interpolated value by eliminating the images, the same would apply to the decimation filter where we want to reject these same image locations prior to throwing samples away). This also converts readily to an efficient polyphase filter structure by mapping the same coefficients row to column in the polyphase filter, as detailed in other posts here.

In this plot the blue is the desired spectrum along with its images, and the red and the black show the multi-band filter response for the two different filter design choices (red is least squares and black is the Parks-McCelllan or equiripple design). This would be equivalent for a decimator except the images would be noise or other signals that could fold in during the process of throwing away samples for decimation. Given the same number of taps observe how the total noise that would fold into band is significantly less with the least squares filter design.

Note: If you have enough samples so as to not truncate the desired response, you could certainly still do the least squares filter design approach in the frequency domain ---- the filtering (convolution) described above that is done in the time domain is equivalent to multiplying in the frequency domain- but to do this properly would necessitate a lot more samples to ensure sufficient tails of the kernel (the frequency transform of the filter's impulse response) are included.

Answer 2

The recommended way of converting an oversampled bandlimited signal into its critical sampling rate (or somehow above that) is to use a time-domain LP filter and decimate approach.

This can be efficiently implemented using a polyphase filterbank architecture as well.

The lowpass filter can be implemented using DFT/FFT frequency domain techniques if the filter length happens to be too long. But that depends on your requirements and it will not improve your compuational accuracy (except that numerical issues may help) but make it more efficient to compute.