Applying an anti-aliasing filter

95 views (last 30 days)
L'O.G.
L'O.G. on 10 Jul 2022
Edited: L'O.G. on 10 Jul 2022
I've been trying to understand how to design and implement an anti-aliasing filter. I'm really stuck. I have an analytic function in the time doman that I calculate according to a model, and then I want to process my results numerically, so I sample at a high rate that I choose to be 10 kHz. Then, I take the FFT like:
dt = tvec(2) - tvec(1);
nfft = 2^nextpow2(length(avg_Gt));
Freq = (1/dt/2)*linspace(0,1,nfft/2+1);
avg_Gw = fft(avg_Gt,nfft)/length(avg_Gt);
avg_Gw = avg_Gw(1:numel(Freq));
I then do a bit of math on the FFT:
avg_result = 2*real(1i*2*pi*Freq.*avg_Gw);
On a log-log plot in the frequency domain, I get a huge dip around 3 kHz that doesn't make physical sense to me:
and that I assume is due to aliasing. The rest of the function makes sense based on what I physically expect.
The signal isn't bandlimited, but since the time domain signal is constructed from an analytic function I feel that it should be treated as exact and that I should be able to filter out (and detect?) the aliased components. Is that indeed the case? What would be an effective way of robustly implementing an anti-aliasing filter in this case?
One thing (among many) that I've tried is a Butterworth filter, e.g.,
[b,a] = butter(4,1000/(10000/2));
dataOut = filter(b,a,avg_Gt);
Then, I take the FFT, etc. as before and get:
I chose the order of 4 arbitrarily, 1000 Hz because that seems to be roughly the cutoff that I want (at least based on plot #1 where the function plateaus), and 10000 because that's the sampling rate of 10 kHz when I take the analytic function and construct the time domain signal. But the plot doesn't seem right because now I don't have the plateau that I was expecting (and seeing in plot #1). It also still has the huge dip that I am suspicious about as being unphysical. Can somebody please help me out?
Attached is a .mat file with my signal (and the time vector) in case somebody can take a closer look. Plotting the time domain data on a log-log plot, like with the FFT data, is what I recommend if you want to visualize it.
  9 Comments
L'O.G.
L'O.G. on 10 Jul 2022
@Star Strider Thanks. The non-monotonic plots that I show are after doing a bit of math in the frequency domain. From my original post:
avg_result = 2*real(1i*2*pi*Freq.*avg_Gw);
where avg_Gw is the FFT of avg_Gt.
Star Strider
Star Strider on 10 Jul 2022
@L'O.G. — I do not understand what you are doing in that assignment, or the reason for it.

Sign in to comment.

Accepted Answer

Paul
Paul on 10 Jul 2022
I think all that we are seeing is the effect of sampling the contiuous signal, nothing to do with aliasing.
Let's look at just one element of ths signal, an exponential decay. For illustration, assume the time constant is 0.1 seconds
syms t real
tau = 0.1;
f(t) = exp(-t/tau)*heaviside(t);
The CTFT of the continuous signal is (all frequencies in rad/sec)
syms Omega real
Fc(Omega) = fourier(f(t),t,Omega)
Fc(Omega) = 
For sampling, select a sampling period of
Ts = 1/1000;
After sampling, the sampled signal is of the form a^n*h[n] where (assuming an infinite number of samples, which can't be done with the DFT via fft() )
a = exp(-Ts/tau);
The DTFT of the sampled signal is (from a DTFT table)
syms omega real
Fd(omega) = 1/(1-a*exp(-1i*omega));
The DTFT approximatly scales to the CTFT by
Fd1(Omega) = Ts*Fd(Omega*Ts);
This approximation is only appropriate for 0 < Omega < pi/Ts.
Compare the two
figure;
fplot(abs([Fc(Omega) Fd1(Omega)]),[0 pi/Ts])
So the amplitudes are a reasonable match. Zoom in at high frequency
figure;
fplot(abs([Fc(Omega) Fd1(Omega)]),[2500 pi/Ts])
The curves don't have quite the same shape.
Now multiply by 1i*Omega as in the Question and plot
Gc(Omega) = 1i*Omega*Fc(Omega)
Gc(Omega) = 
Gd1(Omega) = 1i*Omega*Fd1(Omega)
Gd1(Omega) = 
figure;
fplot(abs([Gc(Omega) Gd1(Omega)]),[0.01 pi/Ts])
set(gca,'YScale','log');
set(gca,'XScale','log');
We see that multiplication amplifies the difference as Omega approaches the Nyquist frequency
Check the limit of the real part as it approaches the Nyquist frequency.
limit(simplify(real([Gc(Omega) Gd1(Omega)])),Omega,pi/Ts)
ans = 
The limit of the CTFT is basically 1, but sampling makes the limit zero.
Plot the real part as in the question
figure;
fplot(real([Gc(Omega) Gd1(Omega)]),[0.01 pi/Ts])
ylim([1e-6 10])
set(gca,'YScale','log');
set(gca,'XScale','log');
We see the CTFT-based plot show a plateau, but the DTFT-based plot decays to zero.
Because the Fourier transform is linear, I would would expect the sum of exponentials to display similar behavior.
  1 Comment
L'O.G.
L'O.G. on 10 Jul 2022
Edited: L'O.G. on 10 Jul 2022
@Paul Thank you! That was really excellent, so insightful.

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!