Applying an anti-aliasing filter

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.

[Image Analyst]

Attach your data here. No one wants to go to some third party web site.

[L'O.G.]

@Image Analyst Done. Thanks for the pointer.

[Paul]

Hi L'O.G.,

What do you expect the plot to look like for frquency > 3 kHz? Stay relatively flat out to 5 kHz?

What is the purpose of this line?

2*real(1i*2*pi*Freq.*avg_Gw);

The multiplication by frequency kind of looks like an attempt to get the Fourier transform of the derivative of the signal; not sure what taking real part is trying to accomplish.

[L'O.G.]

@Paul I don't know what should happen > 3 kHz, but I imagine the huge dip shouldn't be there. All I expect is the monotonic increase and the plateau that I see.

As for the post-FFT math, that's just a physical relation involved in the formula I am calculating.

[Paul]

Can you post the analytic time-domain function that's being sampled?

[L'O.G.]

Sure, it is of the form where I have calculated all constants τ. I construct it in the time domain with a very high sampling rate of 10 kHz, which is much higher than any τ.

[Star Strider]

When I calculate the Fourier transform and plot it, I get a relatively smooth descending curve and no specific breaks —

LD = load('signal.mat');

tvec = LD.tvec;

avg_Gt = LD.avg_Gt;

L = numel(tvec);

Fs = 1/(tvec(2)-tvec(1));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTsig = fft(avg_Gt-mean(avg_Gt), NFFT)/L;

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

avg_Gt_filt = lowpass(avg_Gt, 1000, Fs, 'ImpulseResponse','iir');

FTsig_filt = fft(avg_Gt_filt-mean(avg_Gt_filt), NFFT)/L;

figure

loglog(tvec, avg_Gt, tvec, avg_Gt_filt)

grid

xlabel('Time')

ylabel('Amplitude')

legend('Original','Lowpass Filtered', 'Location','best')

figure

loglog(Fv, abs(FTsig(Iv))*2)

hold on

plot(Fv, abs(FTsig_filt(Iv))*2)

hold off

grid

xline(3E+3,'-r')

text(3E+3, 10, '3 kHz', 'Horiz','right', 'Vert','bottom', 'Rotation',90, 'Color','r')

xlabel('Frequency')

ylabel('Amplitude')

legend('Original','Lowpass Filtered', 'Location','best')

This is the fft code I routinely use.

Filtering the signal does not appear to make any difference in either the time-domain or frequency-domain plots.

.

[L'O.G.]

@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]

@L'O.G. — I do not understand what you are doing in that assignment, or the reason for it.

Answer 1

[Paul]

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.

[L'O.G.]

@Paul Thank you! That was really excellent, so insightful.