Calculating a Power Spectrum in Python via Autocorrelation

Question

I want to calculate the Power Spectrum of a 8192 data vector through the usage of cumulants. I calculated autocorrrelation with 128 max shiftings, reduced it by the signal's mean and performed an fft. The result was complex instead of real and positive. Where did I go wrong?

This is my code.

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
import pandas as pd

#Creating the random variables
φ=[]
for i in range(0,6):
    if(i==2):
        φ.append((φ[0]+φ[1]))
    elif(i==5):
        φ.append((φ[3]+φ[4]))
    else:
        φ.append(np.random.uniform(0,2*np.pi))
        
#Creating the λ variables
λ=[0.12,0.3,0.42,0.19,0.17,0.36]

#Building the x process
x=[]
samples=8192
for k in range(0,samples):
    x.append(0)
    for i in range(0,6):
        x[k]+=np.cos(k*2*λ[i]*np.pi+φ[i])

#Preparing The Plot and adding x to it
fig, [ax1,ax2,ax3,ax4]=plt.subplots(nrows=4,ncols=1)
ax1.plot(x)
ax1.set_title("Time Signal")

#Building The Autocorellation function
lags=128
corr=[]
temp=[]
for i in range(0,lags):
    for k in range(0,samples):
        if(i+k>=samples):
            temp.append(0)
        else:
            temp.append(x[k+i]*x[k])
    corr.append(np.mean(temp))
    corr[i]-=(np.mean(x))**2
    temp.clear()
  
#Calculating The Power Spectrum
C2=np.fft.fft(corr)
print(C2)

Maybe this can help: https://stackoverflow.com/questions/10304532/why-does-fft-produce-complex-numbers-instead-of-real-numbers — Alessandro Togni, Jun 08 '22 at 19:03
Because in theory, the autocorrelation function is symmetric but you have calculated autocorrelation function for only positive lags. — Ömer Şayli, Oct 08 '22 at 07:38

score 0 · Answer 1 · answered Jun 08 '22 at 19:08

0

When performing the Fourier transform with np.fft.fft, the result is always a complex signal which contains the amplitude and phase of each frequency.

What you want is the power spectrum, which can be found by squaring the Fourier transformed signal:

C2= np.abs(np.fft.fft(corr))**2

answered Jun 08 '22 at 19:08

Endre Dåvøy

1
1

But second-order statistics theory states that the Power Spectrum of a discrete x[k] signal IS the Fourier Transform of (autocorrelation-mean), not the squared absolute of it – Wreckord Jun 08 '22 at 20:29
Ah, yes you are correct. The autocorrelation function has a even symmetry so the Fourier transform only contains cosines, i.e the real part. However, when using the FFT, directly squaring the Fourier transform is a faster and more commonly used technique, apparently. https://www.sciencedirect.com/topics/engineering/power-spectrum – Endre Dåvøy Jun 09 '22 at 06:59
The FFT function in numpy will return a complex number, but if the autocorrelation was correct the imaginary components would be negligible. I suspect that since you only calculate the forward delay of the correlation it is no longer even symmetric and thus the Fouier Transform will be complex. You could try to use the hermettian FFT, "np.fft.hfft". It returns a real FFT, assuming that the time function is symmetric – Endre Dåvøy Jun 09 '22 at 07:04

Calculating a Power Spectrum in Python via Autocorrelation

1 Answers1