# difference in the spectral densities of autoregressive sequence by FFT and analytical solution

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

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.

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.8*np.cos(2*np.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()


• I wonder if it could be because "a particular instance of a white noise sequence will not have precisely flat response" ref Mar 14, 2022 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. Mar 14, 2022 at 23:11
• Thank you for doing the extra work! Mar 15, 2022 at 2:18

## 1 Answer

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.


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-2*phi*np.cos(2*np.pi*x))
return PSD

def cal_PSD_fft(N, phi):
# covaraince
gamma = phi**np.arange(N)
yf = fft(gamma).real
PSD = 1 / (1 - phi**2) * (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()


• Nice result for the bilateral PSD! Mine was for the unilateral PSD, hence the factor of 2. (+1)
– Ed V
Mar 17, 2022 at 1:56
• Thank you for +1 and your previous comments Mar 17, 2022 at 2:43