DFT amplitude differs from CFT amplitude

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.

Answer 1

[Paul]

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.

[AN.]

Wahou, thanks a lot for this precise answer ! I see it was wrong to normalize by n the fft(), and not to multiply by T_ech.

However, I do not understand why mathematically we have to do it (n and T_ech). Why in doc fft Cosine Wave example, this is not taken into account?

[Paul]

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.