Fast Fourier Transform Issue

Hello.

In the Matlab description of the FFT implementation (https://www.mathworks.com/help/matlab/ref/fft.html). There is this part where the vector Y is divided by the length L in order to compute the spectrum. I cannot understand why it is required to divide by L, if someone can maybe enlighten this point.

Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.

P2 = abs(Y/L);

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

Thank you.

[dpb]

In short, it's simply to normalize the PSD peak to the amplitude of the input time series -- if don't then the magnitude is proportional to the length of the input.

Just had long discussion (amongst many over the years) a short time ago. The final comment I added at https://www.mathworks.com/matlabcentral/answers/847665-why-do-i-receive-results-mismatch-in-the-fft-signal-while-using-matlab-built-in-fft-function-how-ca?s_tid=srchtitle#comment_1593465 illustrates what happens if one doesn't normalize;

Answer 1

[David Goodmanson]

Hi Anas,

suppose you have N = 1000 points over an interval of one second. The time array is t = (0:999)/1000. [see note below]. Let f0 = 1 Hz and consider the complex waves

g_n(t) = exp(2*pi*i*n*f0*t) for n = 1,2,3, ....

Each g_n(t) oscillates n times in the time window.

The fft takes the signal, and for each m = 1,2,3 ... [1] multiplies the signal by g_m(t)* = exp(-2*pi*i*m*f0*t), [2] does a sum over all the array points and [3] reports out the answer.

Suppose your signal is a single oscillatory wave of amplitude 1 for a particular n0. After step [1] the fft will do a sum over all the points of the expression

    B = exp(2*pi*i*(n0-m)*f0*t).

The result is 0 except when m = n0. In that case B = contant = 1 and the sum over the N points gives the result N. So if you want to recover the original amplitude 1 you have to divide the fft result by N.

[note] The array has 1000 points and 1000 intervals, including the interval from .999 sec to 1 sec, but does not include the point at 1 sec.

[Anas Abid]

I get it better now aside from the fft algorithm part, I can't seem to find a detailed explanation of how it works, if you could maybe help me with that it would be very much appreciated.

Thank you very much for your kind help !

[dpb]

That's the deal -- the algorithm essentially sums N values of unity.

NB: however -- to normalize back to the input scale if augment the FFT by nextpow2 or other way to add extra points to the measured input time series, still divide by ONLY L, the length of the original time series, NOT nFFT since the added points are all just padded zeros and thus contain no power and so don't add to the overall scaling.

[David Goodmanson]

Hi Anas,

do you mean how the fft works in principle, or how it is implemented in software?

[Anas Abid]

Both please

[David Goodmanson]

Every signal can be represented as the sum of complex waves with amplitude coefficients cn:

sig = Sum{n} cn*exp(2*pi*i*n*f0*t) 

The fft determines the cn coefficients, which is the frequency spectrum. To do that multiply the signal by g_m(t)* = exp(-2*pi*i*m*f0*t), take the sum over points and divide by N, which gives, on the right hand side,

Sum{n} cn*exp(2*pi*i*(n-m)f0*t) /N  

AsI mentioned before, the sum is N when n = m and zero otherwise. You get N*cm/N = cm. Do that for each m. That's the principle. As for the underlying software, that's proprietory with Matlab, but if you look up the Cooley-Tukey algorithm you will get the basic idea.