Questions about FFT (and applying it to determine power spectral density)

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!

[Yongcui Mi]

Regarding Quetion 2, both the real and complex parts of the amplitude should be considered, which is done by x .* conj(x)

Answer 1

[Daniel kiracofe]

"Question 1: Why isn't the scaling of ( 2 / numberOFdataPOINTS ) applied within the FFT algorithm."

Matlab uses what I like to call the "mathematician's FFT". The mathematical definition of the discrete fourier transform (i.e. if you look up the theory on wikipedia) does not include that scaling. However, what we really want is the "engineer's FFT". With the right scaling and whatnot to make it useful. You'll need a wrapper script to make matlab's fft() useful. To me, it is a serious shortcoming of matlab that no such wrapper script is provided as part of the standard package. Here is my standard wrapper script which you main find useful: http://www.mechanicalvibration.com/Making_matlab_s_fft_functio.html

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?

Not that I know of. But if you want to do power spectral density, you are much better off using something other than just abs(fftSignal).^2. The key thing to keep in mind is that practically speaking you can never ever calculate power spectral density, you can only estimate it. And using a plain fft(), while it does give an estimate, it is not a good estimate. There are much better estimators available. Some tutorials I wrote on power spectral density and better estimation method: http://www.mechanicalvibration.com/Introduction_power_spectral.html

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?

Mathematically, you can take the fourier transform of a complex (i.e. real + imaginary) signal. In that case, the positive and negative frequency components are not redundant. But if you have a real valued signal, then the negative frequency components are always redundant and give exactly the same information as the positive frequency components. i.e. both sides of the fourier space contain the exact same set of information. Therefore one side can be discarded.

zero padding is a separate question in itself. If you are trying to estimate power spectral density, I don't think you should be using zero padding to improve the resolution. Instead you should use a better method (like welch's method, described in my website above), rather than zero padding.

Hope some of this helps.

[Nicholas Trunfio]

Thank you very much, the wrapper is a big help and the explanation (plus some background in the links) really cleared up the issues I was having with the FFT implementation.

Now that I'm more confident that I'm applying the FFT correctly, I can start to figure out why the PSD still appears incorrect (even when using a built in function like periodogram).

[maryam]

Hi Daniel,

Many thanks for your explanations, they were pretty helpful to me.. I just have one question, in your wrapper, when converting to one sided spectrum, you take half of data and multiply all the points except first one by two % convert from 2 sided spectrum to 1 sided %(assuming that the input is a real signal) amp(1) = abs(z(1)) ./ (n); amp(2:halfn) = abs(z(2:halfn)) ./ (n / 2); In Mathworks help though, they took xdft(1:N/2+1) and then 2*psdx(2:end-1) as is in (<http://www.mathworks.com/help/signal/ug/psd-estimate-using-fft.html>) I was wondering what is the difference between these two approaches? Thanks

Answer 2

[Rick Rosson]

These questions have been asked and answered dozens of times already. Please search MATLAB Answers for fourier and/or fft.