Implementing an anti-alliasing filter in Mathematica

Question

I'm experimenting with FFT simulations. I'm generating transient data at the some sampling rate, down-sampling this (mimicking a device's behavior) to 2048 points. Then take the FFT using Fourier[Data,FourierParamters->{-1, -1}]. However I realise I need to apply an anti-aliasing filter to the transient signal to remove higher frequency components that violate the Shanning-Nyquist theorem. Does anyone know of an approach to implement this in Mathematica?

Here is my sample code

TransientSignal[A_, ω_, t_] := A Cos[ω t]

A0 = 0.4;
ν0 = 21.;
ω0 = 2. π ν0; 
C0 = 0;


Sr = 262100;
dt = 1. / Sr;

FFTSpan = 50;
FFTLines = 800.;
FFTSamples = 2048;
FFTCentre = 24000.;

Δf = FFTSpan / FFTLines;
τAvg = 1. / Δf;

TransPts = τAvg / dt;

DownSample = Round[128/FFTSpan FFTLines];



TransientSignalData = TransientSignal[A0, ω0, Subdivide[0., τAvg, TransPts - 1]];
    DownSampledData = Downsample[TransientSignalData + TransientNoise * 0, DownSample];

FFT = Fourier[DownSampledData,FourierParameters->{-1, -1}];
Freq = Table[i Δf, {i, 0, FFTSamples - 1}];
    AbsFFTwithFreqs = Transpose[{Freq, 2. Abs[FFT]}];

ListLinePlot[AbsFFTwithFreqs,PlotRange->All,Frame->True,FrameLabel->{{"AMPLITUDE ARB UNITS",""},{"FREQ. HZ",""}},ImageSize->500]

The reason I think I need this filter is that, if one keeps increasing the frequency, ν0, of the input transient signal you notice that the two peaks of the resultant FFT converge and wrap around cyclically as ν0 is increased. This doesn't make sense as once ν0 moves out of the frequency span I shouldn't see any peak features because of Shanning-Nyquist. So I need an AA filter to remove these higher terms for things to make sense.

Answer 1

The aliasing is more easily described by animating your code. The only changes that I've made are:

• Set the amplitude to 1.0 rather than 0.4 (for a nicer plot).
• Change PlotRange -> All to PlotRange -> {0,1} (for animation on fixed axes).
• Replace the line \[Nu]0 = 21; by \[Nu]0 = freq; (to define plot[freq]).

The animation:

plot[freq_] := (   (*your code here*)  ) 

Animate[(plot[freq]), {freq, 0, 512, 0.001}, DefaultDuration -> 40]

It does make sense (64 is the Nyquist frequency for this example), but you are correct to suppose that your results will not make sense unless you filter. Beyond the Nyquist frequency, sampling aliasing kicks in, which is exactly what you see in the animation: the peaks cross (or reflect?) at freq = 64, and wrap (or bounce?) at freq = 128, ... , and this happens all the way to freq = Infinity.

If TransientSignal contains frequencies above freq = 64, then the spectrum you get from (discrete) FFT is going to be contaminated with aliased Fourier components as soon as you take the FFT. So these need to be removed before taking FFT by use of a down-sampling digital filter. In fact, you are already using one, but Downsample is rather a poor choice, because, as you have already observed, aliased Fourier components do not get suppressed by Downsample. But at least they don't get amplified :)

A very much better (but far from optimal) down-sampling filter is obtained by taking the average of the samples that you threw away. Eg, for this data,

o . . . . . . o . . . . . . o . . . . . . o . . . . . . o

you could down-sample by averaging each "o . . . . . . o" rather than just keeping "o" and throwing away all the useful information in " . . . . . . ".

BTW, the above is clearly N times more efficient than a moving average, where N + 1 is the number of samples in "o . . . . . . o".

To try out this averaging idea (which should be regarded as only a first step towards something better suited to your particular problem), let's define that filter,

SimpleDownsampleFilter[s0_List, n_Integer] := Module[{},
len = Length[s0];
size = len - Mod[len, n];
s = Take[s0, size];
Map[Total, Partition[s, n]]/n]

To us this filter instead of Downsample, replace the line of code DownSampledData = Downsample[ ... ] (which is in your code, inside plot), by this,

DownSampledData = SimpleDownsampleFilter[TransientSignalData+TransientNoise*0,DownSample];

(* Small fudge here: got 2047 samples but want 2048, so assume periodic and add another one *)
AppendTo[DownSampledData, DownSampledData[[1]]];

Now run the animation again, but for a longer time and to higher frequencies, so that we can see how the aliased Fourier components now slowly die away,

Animate[(plot[freq]), {freq, 0, 1024, 0.0001}, DefaultDuration -> 120]

How to improve on this averaging idea?

It would be nice if the aliased components died away (as freq increased) a lot faster than they do for SimpleDownsampleFilter. There is a lot of freedom here that can be exploited to achieve that:

1. Instead of using an ordinary average, one can use a weighted average. In this particular example, that gives us 1024 weights that we can choose arbitrarily.

2. Even if one was constrained to ordinary averaging (e.g., for speed), one could still apply digital filtering to the down-sampled data.

If you have a particular system that you need to deal with, then the weights and coefficients for averaging and filtering should be optimized to deal with the noise sources in your system. Otherwise, IMHO the best way to go would be to use my Convolution Splines (or at least formulate the problem in a similar way).

Convolution Splines are just very cool digital filters, so don't let the name put you off. Here is a .zip of the 5 papers that I wrote (many years ago) that explain what they are and how to use them,

https://ahnorton.com/uploads/Convolution-Splines-Digital-Filters-and-Aliasing.zip

The context for their development was the problem of suppressing unwanted aliased frequencies that are generated by non-linear terms in PDEs while using spectral collocation methods (as opposed to filtering out noise from a signal). But you can just skip those bits. The notebook I used for this post is here,

https://ahnorton.com/uploads/FFT-anti-aliasing.nb