Power spectrum incorrectly yielding negative values

Question

I have a real signal in time given by:

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:

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:

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.

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?

Answer 1

This is due to an implicit time shift, which corresponds to a phase shift, perceived here as a sign inversion.

The signal you are Fourier transforming is symmetric around t=0 and this is why you should expect a positive power spectrum. However, the discrete Fourier transform (DFT) of a time-series $x_n$, which is what the FFT alorithm implements, is defined as $$ X_k = \sum_{n=0}^{N-1}x_n e^{-2\pi i k n}, $$ where $N$ is the number of data points.

Ignoring some normalization issues, this can be interpreted as an approximation of the Fourier transform of a sampled signal, where $$\omega_k = 2\pi k$$ and $$t_n = \frac{n}{T},$$ $T$ being the width of the time window. This is the source of the implicit time shift. The "correct" time sample values are $$\hat t_n = t_n - T/2,$$ so that the phase corrected power spectrum, which has real and non-negative values, is given by $$ \hat X_k = \sum_{n=0}^{N-1}x_n e^{-2\pi i k (n-1/2)} = X_k \cdot e^{ i \pi k},$$ or $$X_k=(-1)^k \hat X_k$$, which explains the fast sign change that you observe.