2
$\begingroup$

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.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() 

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

1 Answer 1

1
$\begingroup$

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-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()

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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