2

I want to obtain the power spectral density (PSD) of an autoregressive sequence, AR(1). The analytical solution according to this reference (page 12) is

For $X_t = \phi_1X_{t-1}+W_t, W_t \sim N(0,\sigma^2_w)$，

$$f(w) = \frac{\sigma_w^2}{1-2\phi_1 \cos(2 \pi w)+\phi_1^2}$$

Let $\sigma_w = 1$, $\phi_1 = 0.4$, then the corresponding PSD curve in [0, 2pi] using fft and analytical solution is

enter image description here

As can be seen from the figure below, the PSD obtained from fft is far from being a U-shaped curve obtained by the analytical solution.

enter image description here

Can anyone enlighten me on why the results are different?

The code I used to generate the plot is:

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq

def analy(x):
    PSD = 1/(1.16-0.8np.cos(2np.pi*x))
    return PSD


def my_fft(N):

    T = 1  # sample spacing
    t = np.linspace(0.0, N*T, N, endpoint=True)
    y = np.zeros([N,1])


rho = 0.4
for i in range(len(t)-1):
    y[i+1] = rho*y[i] + (1-0)*np.random.normal(0,1) 

yf = fft(y)   
xf = fftfreq(N, T)[:N//2]  # sample frequency points
PSD = 2.0/N * np.abs(yf[0:N//2])

return xf*np.pi*2, PSD




def compare_PSD():
    N = 100     # Number of sample points
    UB = np.pi  # sample upper bound
    x_analy = np.linspace(0, UB, num=N)
    PSD_analy = analy(x_analy)
    x_fft, PSD_fft = my_fft(N)
    plt.plot(x_analy, PSD_analy, label='Analytical')
    plt.plot(x_fft, PSD_fft, label='scipy.fft')
    plt.xlabel('Frequency')
    plt.ylabel('PSD')
    plt.legend()
    plt.show()


compare_PSD()
Jayyu
  • 131
  • 5
  • I wonder if it could be because "a particular instance of a white noise sequence will not have precisely flat response" ref – Jayyu Mar 14 '22 at 18:55
  • Try zero padding your FFT out to ten times longer or more. This will interpolate the frequency spectrum which may provide a result closer to what you are expecting. – Dan Boschen Mar 14 '22 at 23:11
  • Thank you for doing the extra work! – Jayyu Mar 15 '22 at 02:18

1 Answers1

1

I got the correct results. There are two mistakes in my previous comparison.

  1. Conceptual error: PSD is obtained by taking the Fourier transformation on the covariance, not the original time series data.

$$ \begin{aligned} f(v) &=\sum_{h=-\infty}^{\infty} \frac{\sigma_{w}^{2}}{1-\phi_{1}^{2}} \cdot \rho(h) e^{-2 \pi i v h} \\ &=\frac{\sigma_{w}^{2}}{1-\phi_{1}^{2}}\left[2 \sum_{h \in(0, \infty)} \phi_{1}^{h} e^{-2 \pi i v h}- \phi_1^{0} e^{0}\right] \end{aligned} $$

  1. Coding: "scipy.fft" takes a vector as its input. However, "y" in the previous code is a matrix.

Below is the correct result. enter image description here


import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq

def cal_PSD_anal(x, phi):
    PSD = 1/(1 + phi*2-2phinp.cos(2np.pi*x))
    return PSD


def cal_PSD_fft(N, phi):

    # covaraince 
    gamma = phinp.arange(N)
    yf = fft(gamma).real
    PSD = 1 / (1 - phi2) * (2*np.abs(yf[0:N//2]) - 1)

    return PSD


def compare_PSD():
    N = 100     # Number of sample points
    UB = 0.5 #np.pi  # sample upper bound
    T = 1  # sample spacing
    phi = 0.4


x = fftfreq(N, T)[:N//2]  # sample frequency points
PSD_fft = cal_PSD_fft(N, phi)    

PSD_anal = cal_PSD_anal(x, phi)

plt.figure(figsize=(4, 3), dpi=300) 
plt.plot(x, PSD_anal, 'bo', label='Analytical', markersize=4)
plt.plot(x, PSD_fft, 'r', label='FFT') 

plt.title('PSD of AR(1, {})'.format(phi))
plt.xlabel('Frequency')
plt.ylabel('PSD')
plt.legend()
plt.grid(axis='y')
plt.tight_layout()
plt.savefig('./PSD_AR_1_{}.png'.format(phi), dpi=300)
plt.show()




compare_PSD()
Jayyu
  • 131
  • 5