Fourier Transform with Matlab

TL;DR (summary)

It is very hard to find online examples of fft usage with Matlab that normalise the amplitude/power values correctly (e.g. as can be verified with Parseval's theorem). This is crucial if you want to compare spectra between signals with different lengths. There is also an additional problem with real-valued signals because in that case the spectrum is often computed for positive frequency only, and therefore amplitude or power values need to be scaled to account for frequency-folding. Following the post and answers below, here is a gist which I think scales the coeffcients correctly and consistently for real- and complex-valued inputs.

The take home messages are:

• Do NOT normalise the DFT coefficients directly (eg don't write Y = fft(x)/L);
• If you use an n-points transform (eg fft(x,nfft)), then the normaliser is nfft and not numel(x);
• If you extract a single-sided spectrum, you need to adjust amplitude/power values depending on which correspond to conjugate pairs DFT coefficients;
• If you extract a single-sided spectrum, you should compute the amplitude and power separately (ie don't compute the power from the amplitudes).

---

Amplitude, power and single-sided spectra

As defined and explained on Wikipedia:

• The DFT coefficients are complex and not normalised, while the formula for the inverse DFT carries a 1/N factor in front of the sum. This is natural in some sense, as moving in time-to-frequency direction can be seen as a projection onto a basis of (orthogonal) waves with different frequencies, whereas moving in frequency-to-time direction can be seen as a weighted superposition of waves.
• In that superposition, the overall magnitude should be that of the original time-point (ie, it's an inversion), and the amplitude of each wave in that weighted average is simply the magnitude of the corresponding DFT coefficient divided by the number of waves |Xk| / N. Similarly, the power of each wave is |Xk|^2 / N. Matlab uses that normalisation too (well, FFTW does).
• For real-valued inputs, the DFT coefficients are conjugate pairs for corresponding positive/negative frequencies, apart from the DC component (average term, frequency 0), and for the Nyquist frequency when the number of time-points is even. In practice, most applications avoid this redundancy by extracting the DFT coefficients only for positive frequencies. This leads to complications in the discrete values of the amplitude and power.
• The amplitudes corresponding to paired DFT coefficients (all except the first one and the Nyquist frequency when it exists) can simply be doubled and the negative frequency then discarded. Same for the power.
• Similarly for the power, but note that it is incorrect in that case to compute the discrete power values using the adjusted amplitude values. Indeed N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2 is not equal to 2*|Xk|^2 / N (this is where the square-root of two came from in the OP). Therefore it is necessary to compute independently the amplitude and power values from the DFT coefficients (another good reason not to scale them).

N-points Transform

Many of the examples online use an explicit N-points transform: Y = fft(x,NFFT) where NFFT is typically a power of 2, making the computation more efficient with FFTW.

The effective difference in that case (provided that NFFT >= N) is that x is padded with 0 at its end until it reaches the length of NFFT time-points. This means that the number of frequencies in the decomposition changes, and therefore the normalisation should be done relative to NFFT wave components, and not the original N time-points.

Hence, almost all of the examples found online are wrong in the way they normalise the coefficients. It should not be Y = fft(x,NFFT)/N, but Y = fft(x,NFFT)/NFFT -- yet another good reason to loose that habit of normalising complex coefficients.

Note that this makes no difference then to Parseval's equality, because the added terms in the time-domain are all zero, and therefore their contribution to the now-larger sum is zero too. In the frequency domain though, the added discrete frequencies will have a response to the original signal in general, which gives an intuitive sense of why the obtained coefficients can be quite different indeed between the padded and unpadded transforms.

Code

The code in the OP is therefore incorrect, and instead it appears to be necessary to output both the amplitude and power, as there is no general normalisation coefficient that could accommodate the complex and real cases, with even or odd number of time-points. You can find the Gist here.