Hi Francesco, For an N-point fft, the usual way to set up the time and frequency arrays of step size dt and df is
dt = 1/Fs; % inverse of sample rate
df = Fs/N; % total span of frequency array is Fs
dt and df are usually thought of in terms of the comments above, and are consistent with
df*dt = 1/N. [1]
This is a fundamental relationship for the fft, and I think it deserves more attention.
For signal y, fft(y) / N gives the correct amplitudes. For simplicity I will use y = exp(2*pi*i*f0*t). Its fft has a single peak of amplitude 1 at f0, compared to sin or cos which have two frequency peaks of amplitude 1/2 at +-f0.
In the first example below (fig 1), N = 1000 and everything works (as with the first version of your code). There are exactly 60 cycles in the time window. df is 1 Hz. The fft has a single peak at 60 Hz, zeroes everywhere else.
Fig 2 shows zerofilling of the time domain signal, with Nfft = 3000. The spacing dt stays the same, but N has changed. From [1], df becomes 1/3 Hz.
For normalization, the fft is divided by 3000.
At least 60 Hz is a multiple of 1/3 Hz, and the expanded view of the fft shows a peak right at 60 Hz, but the amplitude is only 1/3. However, the amplitude is perfectly correct. That's because the wave only covers a third of the new time array, i.e. the duty cycle is 1/3. That’s not so bad, but there are other problems. Since the wave is cut off, you don't have a continuous periodic wave anymore. So there are contributions at a lot of frequencies other than 60 Hz as shown.
Sometimes to get out of this, people divide by the original N rather than the new N. That gets the amplitude back up to 1, but does nothing about the nonzero amplitudes at other frequencies. It's a bit of a con job.
The last example in Fig. 3 is the highly regarded nextpow2 method, with Nfft = 1024. Now there are all the same problems as before, but from [1] the new df is 1/(Nfft*dt) = 1000/1024 = .9765625 Hz. There is no way that 60 Hz is a multiple of that, and that frequency array can't even represent the peak correctly. You could reverse engineer things to find the factor needed to make the largest subpeak equal to 1, but the result is meaningless.
So why zerofill to 1024? Supposedly it's all worthwhile because of the big speed advantage, but in my opinion any speed advantage is highly oversold. I think that in many many situations the speed gained is minuscule and unnecessary, and the disadvantages of nextpow2 are evident.
I ran each version of your code (except plotting) 50k times in a for loop with tic and toc. The first version took 7 seconds, about 140 microsec each. If it’s that fast anyway, how much useful time will you save by using nextpow2? It’s still worthwhile, though, to see how much you improve on the time. Running the nextpow2 version 50k times takes 8.7 seconds. The fast version is actually slower. And that’s not atypical of what can happen in the real world. Evidently the added complication costs more than it gains.
Sometimes one might need to zerofill to make sample lengths match up, or to avoid having N be a large prime number since the fft is slow in that case. And if the frequency grid is dense enough and the signal is somewhat of a mess anyway, then sometimes a bit of zerofilling doesn’t make much difference. But I do not understand why people overuse it as much as they do.
N = 1000;
dt = .001;
t = (0:N-1)*dt; % time span 1 sec
df = 1/(N*dt); % df = 1 Hz
f = (0:N-1)*df; % freq span 1 kHz
f0 = 60;
y = exp(2*pi*i*f0*t);
figure(1)
subplot(2,1,1)
plot(t,real(y))
z = fft(y)/N;
subplot(2,1,2)
plot(f,abs(z),'o')
% zerofill cases
Nfftval = [nan 3000 1024];
for k = 2:3
Nfft = Nfftval(k);
t = (0:Nfft-1)*dt; % longer time span
df = 1/(Nfft*dt)
f = (0:Nfft-1)*df;
y = exp(2*pi*i*f0*t);
y(1001:end) = 0; % zerofill remaining signal
figure(k)
subplot(2,1,1)
plot(t,real(y))
z = fft(y)/Nfft;
subplot(2,1,2)
plot(f,abs(z),'o')
xlim([40 80])
end
