Power spectrum incorrectly yielding negative values

1
votes

I have a real signal in time given by:

enter image description here

And I am simply trying to compute its power spectrum, which is the Fourier transform of the autocorrelation of the signal, and is also a purely real and positive quantity in this case. To do this, I simply write:

import numpy as np
from scipy.fftpack import fft, arange, rfftfreq, rfft 
from pylab import *

lags1, c1, line1, b1 = acorr(((Y_DATA)), usevlines=False, normed=True, maxlags=3998, lw=2)
Power_spectrum = (fft(np.real(c1)))
freqs = np.fft.fftfreq(len(c1), dx)
plt.plot(freqs,Power_spectrum)
plt.xlabel('f (Hz)')
plt.xlim([-20000,20000])
plt.show()

But the output gives:

enter image description here

which has negative-valued output. Although if I simply take the absolute value of the data on the y-axis and plot it (i.e. np.abs(Power_spectrum)), then the output is:

enter image description here

which is exactly what I expect. Although why is this only fixed by taking the absolute value of my power spectrum? I checked my autocorrelation and plotted it—it seems to be working as expected and matches what others have computed.

enter image description here

Although what appears odd is the next step when I take the FFT. The FFT function outputs negative values which is contrary to the theory discussed in the link above and I don't quite understand why. Any thoughts on what is going wrong?

2
numpyscipysignal-processingfftspectrum
maybe its an fft.fftshift issue? – Paul Panzer
@PaulPanzer Do you mind expanding? For example, are you saying the FFT is not considering the time-axis of the autocorrelation correctly? In the above code, is it not simply assuming constant interval spacing between the items in np.real(c1) with the same indexing as c1 (i.e. element 0 is the first one, element 1 is the second, etc.)? – Mathews24
Well, since you seem to be saying that the sizes of the spectrum are ok but not the complex arguments; that is precisely the effect of a shift in the time domain. If your autocorrelogram is centred then that is wrong because if I'm not mistaken zero in "fft coordinates" is not in the centre but on the very left. Since this discrepancy is a common headache there is this convenience function fftshift. – Paul Panzer
@PaulPanzer Thank you for the clarification. So it does appear to work and give me the same answer as simply taking the absolute value when I change the x-axis to np.fft.fftshift(freqs) and the y-axis to fftshift(fft(ifftshift(c1))). The documentation does provide some information, although I'm still a little uncertain on why exactly applying ifftshift and then fftshift fixes the problem. Any thoughts? – Mathews24

2 Answers

0
votes

The power spectrum is the FFT of the autocorrelation, but that's not an efficient way to calculate it.

The autocorrelation is probably calculated with an FFT and iFFT, anyway.

The power spectrum is also just the squared magnitude of the FFT coefficients.

Do that instead so that the total work will be one FFT instead of 3.

0
votes

An fft produces a complex result (real and imaginary components to represent both magnitude and phase of the spectrum). You have to take the (squared) magnitude of the complex vector to get the power spectrum.