DFT amplitude differs from CFT amplitude

6 views (last 30 days)
AN.
AN. on 31 May 2022
Commented: Paul on 1 Jun 2022
Hi,
I compute the FFT of an odd step signal :
T_ech = 5E-5;
y = [0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0];
t = [T_ech:T_ech:T_ech*length(y)];
L = length(y);
n = 2^nextpow2(L);
Y = fft(y,n);
P3 = abs(Y/n);
P3 = P3(1:n/2+1);
P3(2:end-1) = 2*P3(2:end-1);
f3 = 1/T_ech*(0:n/2)/n;
plot(f3,P3)
I try to compare with the continuous FT:
tau = 4.5E-5; % y is 1 during 9 samples at T_ech.
Sf = tau*sinc(pi*f3*tau);
plot(f3,abs(Sf),'-')
Unfortunately, amplitudes are not the same, and depends of n when calculated by DFT. Do you know why?
Thanks.

Accepted Answer

Paul
Paul on 31 May 2022
Hi An,
A couple of issues with the code, addressed below.
T_ech = 5E-5;
y = [0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0];
t = [T_ech:T_ech:T_ech*length(y)];
It's a good idea to have the first point in y correspond to the sample at t = 0, athough it won't matter here because the question is about only the amplitude of the FT.
t = [0 t];
y = [0 y];
Let's plot it
figure
stem(t,y)
L = length(y);
n = 2^nextpow2(L);
Y = fft(y,n);
So far so good. At this point, thiere is no need to normalize by n.
%P3 = abs(Y/n);
P3 = abs(Y);
Looks lke we are only interested in the DFT for positive frequencies. With n an even number, the positive frequencies (assuming pi maps to a negative frequency consistent with Matlab convention) are 1:(n/2)
%P3 = P3(1:n/2+1);
P3 = P3(1:(n/2));
Because we just want to match the amplitude of the CTFT, there is no need to multiply by 2.
%P3(2:end-1) = 2*P3(2:end-1);
The frequency points, in Hz, that correspond the samples of P3 are then
%f3 = 1/T_ech*(0:n/2)/n;
f3 = 1/T_ech*(0:(n/2-1))/n; % Hz
Now the CTFT of the signal is
syms tt ww real % ww in red/sec
yy(tt) = rectangularPulse(3e-4,7e-4,tt);
figure
fplot(yy(tt),[0 1e-3])
ylim([0 1.1])
YY(ww) = simplify(fourier(yy(tt),tt,ww),100)
YY(ww) = 
Matlabs definiton of sinc(t) is sin(pi*t)/(pi*t). So we see that the CTFT can be expressed as
tau = 4e-4; % window length
% tau = 4.5E-5; % y is 1 during 9 samples at T_ech.
YY1(ww) = tau*sinc(ww*tau/2/pi)*exp(-1j*ww/2000)
YY1(ww) = 
Verify YY1(ww) and YY(ww) are equivalent
YY(ww) - YY1(ww)
ans = 
0
Define a continuous frequency vector for evaluating the CTFT
f3c = linspace(0,f3(end),500);
The function Sf, computed in terms of Hz is
% Sf = tau*sinc(pi*f3*tau);
Sf = tau*sinc(f3c*tau);
Now plot. Note that we have to scale the DFT by the sampling period
% plot(f3,abs(Sf),'-')
figure
plot(f3c,abs(Sf));
hold on
stem(f3,P3*T_ech)
xlabel('Hz')
The reason that the plot still doesn't look like a great match is because two samples of y lie on the leading and trailing edges of the rectangular pulse. Actually, the question didn't explicitly define the underlying continuous signal, so that's an assumption on my part based on what I thought the code in the question was tyring to do. So with this assumption, we can get a much better match by treating the edge samples as 1/2
y(find(y==1,1)) = 0.5;
y(find(y==1,1,'last')) = 0.5;
Y = fft(y,n);
P3 = abs(Y);
P3 = P3(1:(n/2));
figure
plot(f3c,abs(Sf));
hold on
stem(f3,P3*T_ech)
xlabel('Hz')
Now we have a good match at low frequencies, but a small divergence at high frequencies. Feel free to follow up if you'd like more information about that.
  2 Comments
Paul
Paul on 1 Jun 2022
In that example in the doc page for fft, they are scaling by L, the length of the signal, not n, the length of the fft. I believe they are doing that in that example so that the peaks of the fft's for each signal should be the same as the amplitude of the underlying cosine waves, which is unity in this case. But our goal here is to match the CTFT, not an amplitude of a cosine wave.
Recall that the DFT (as implmented by fft()) are frequency domain samples of the DTFT, so what we are really interested in is the relationship between the DTFT and CTFT. This relationship can be derived and in doing so we'd see the scaling by T_ech fall out. A way to see this intuitively is that the CTFT is an integral with respect to time, i.e., an area under a curve. Imagine if we were to approximate that integral using a rectangular rule, with each rectangle of width T_ech. Because each rectangle is the same width, the approximation to the integral becomes T_ech*sum(rectangle heights). That sum is essentially the DTFT. So we have CTFT is apprximated by T_ech*DTFT. This approximation is one of the reasons why the plots above diverge a bit at high frequencies, and this divergence can be eliminated by using the exact relationship between the CTFT and the DTFT. The other problem with the intuitive explanation is that it doesn't really justify why I set the edge samples to 1/2. But I hope the intuitive explanation at least give an idea of why I multlplied the DFT samples by T_ech.

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by