TL, DR section:
Question 1: Why isn't the scaling of ( 2 / numberOFdataPOINTS ) applied within the FFT algorithm.
Question 2: Is there any reason to use abs(fftSignal).^2 over fftSignal .* conj(fftSignal) to convert from amplitude spectral density to power spectral density?
Question 3: Which portion of the Fourier Space corresponds to the power spectral density, and what is the best way to zero pad the results to maintain resolution in this portion?
-----------------------------------------------------------------------------------------------------------------------------------------------
Intro
I've been trying to develop a script to show what the Autocorrelation Function (ACF) and Power Spectral Density (PSD) of Gaussian White Noise (GWN or WN) looks like. I expected the ACF to only show a strong correlation when there is no time shift, and this is the result I obtain when using the xcorr function. I also expect the PSD of the GWN to look like a (roughly) straight line for all frequencies that has an amplitude of the variance of the GWN divided by two. However, when I try and determine the PSD by taking the Fourier Transform (using the FFT function) of the ACF I do not get the expected result. I was hoping that someone could shed some light on the use of the FFT function because several sources use the same scaling factors without really explaining where they come from (and in different orders).
The following bullet points indicate what I think I know about the FFT algorith, but please correct me if anything is wrong.
1.) The values returned by FFT are the complex fourier coefficients
> ftSignal = fft(signal)
2. Taking the absolute value of the FFT results gives the magnitude of the fourier coefficients
> ftMag = abs(ftSignal)
-----------------------------------------------------------------------------------------------------------------------------------------------
Question 1:
I couldn't find it in the FFT documentation, but this website indicates what equation the FFT algorithm is using (Go to "click here" for a course outline. Then navigate to "week 7", then to "Numerical calculation of Fourier Series Coefficients (FFT)"). The author suggests that the algorithm that MATLAB uses returns the coefficients, but that they are off by a factor of (2 / numPts). Should the output of the Fourier Transform always be scaled by (2 / numPts) to generate the true coefficients?
> realFtSignal = (2 / length(signal)) * fft(signal);
In most of the mathworks documentation that I have viewed sometime between when the FFT is calculated and when it is plotted this scaling is introduced. For instance, in the FFT documentation the values returned from fft(signal) are immediately scaled by the number of data points; however, only after converting to magnitude is it scaled by the factor of 2. In the PSD estimation documentation is scaled by (2 / numPts) after converting to the magnitude. If the scaling always needs to be done, why is it not inculed in the FFT function itself? Also, why wasn't it most appropriate to apply the scaling right after doing the FFT in both articles (as in, why wait until after converting to magnitude)? (I am VERY confused about this 2 / N scaling as it is mentioned in many places, but is absent from the documentation leading me to wonder if this has been corrected in later releases of MATLAB)
-----------------------------------------------------------------------------------------------------------------------------------------------
Quetion 2:
I've seen it stated that without squaring, the magnitude of the fft is simply the amplitude spectral density; if we want the power spectral density then we need to square the magnitude. Is this correct?
Assuming the above is correct, is there any computational efficiency (time or accuracy) gained by using abs(x)^2 vs x .* conj(x)? Or was it simply different programming style associated with whoever wrote the documentation?
-----------------------------------------------------------------------------------------------------------------------------------------------
Question 3:
In the documentation, the second half of the values returned from the FFT are kept, while the first half of the values are discarded. Why is it that the Fourier transform produces a symmetric result that necessitates this? Is this equivalent to doing fftshift(ftSignal) and keeping just the first half?
Now I know that the algorithm performs best if the number of frequency bins is a multiple of 2. Therefore, we can do the following:
> nfft = 2^nextpow2(length(signal)); > ftSignal = fft(signal, nfft);
to ensure that the signal is zero padded before finding the FFT. However, will this affect the scaling that was discussed in question 1? Which of the following 2 commands would be correct:
> ftSignal = (2 / length(signal)) * fft(signal, nfft);
or
> ftSignal = (2 / nfft) * fft(signal, nfft);
I'm thinking the first one, because padding with zeros shouldn't affect the number of nonzero terms in the summation, but I just wanted to doublecheck to be sure.
These two articles on stack exchange discuss the benefits of padding the signal with zeros to be double its length (the next power of 2 of double its length, actually). They claim that doing so will leave the half of the PSD that is kept with the same resolution as the ACF that it was determined from. Is this an acceptable way to improve resolution? Or will this have other costs associated with it (like degraded numerical accuracy)?
I realize this has been a really long question, and I appreciate any time that you take to help me out!
