1
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.