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.
9 个评论
Image Analyst
2022-7-10
Attach your data here. No one wants to go to some third party web site.
L'O.G.
2022-7-10
Paul
2022-7-10
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.
2022-7-10
Paul
2022-7-10
Can you post the analytic time-domain function that's being sampled?
L'O.G.
2022-7-10
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 τ.
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
2022-7-10
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.
2022-7-10
Star Strider
2022-7-10
@L'O.G. — I do not understand what you are doing in that assignment, or the reason for it.
采纳的回答
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)
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)
Gd1(Omega) = 1i*Omega*Fd1(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)
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Simulation, Tuning, and Visualization 的更多信息
Translated by 
