信息处理 - FFT和解析解的自回归序列的光谱密度差异 - 吾爱随笔录

FFT和解析解的自回归序列的光谱密度差异

信息处理 fft 功率谱密度时频时间序列自回归模型

2022-02-24 11:13:21

我想获得自回归序列 AR(1) 的功率谱密度 (PSD)。根据此参考资料（第 12 页）的分析解决方案是

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

f (w) = \frac{σ_{w}^{2}}{1 - 2 ϕ_{1} \cos (2 π w) + ϕ_{1}^{2}}

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

设 ,，则 [0, 2pi] 中使用 fft 和解析解对应的 PSD 曲线为 $\sigma_w = 1$ $\phi_1 = 0.4$

从下图可以看出，fft得到的PSD远不是解析解得到的U型曲线。

谁能告诉我为什么结果不同？

我用来生成情节的代码是：

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

1个回答

我得到了正确的结果。我之前的比较有两个错误。

概念错误：PSD是通过对协方差进行傅里叶变换得到的，而不是原始时间序列数据。

\begin{aligned} f (v) & = \sum_{h = - \infty}^{\infty} \frac{σ_{w}^{2}}{1 - ϕ_{1}^{2}} \cdot ρ (h) e^{- 2 π i v h} \\ = \frac{σ_{w}^{2}}{1 - ϕ_{1}^{2}} [2 \sum_{h \in (0, \infty)} ϕ_{1}^{h} e^{- 2 π i v h} - ϕ_{1}^{0} e^{0}] \end{aligned}

$\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}$

编码：“scipy.fft”将向量作为输入。但是，前面代码中的“y”是一个矩阵。

下面是正确的结果。


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

其它你可能感兴趣的问题

上一篇什么时候是 Schwartz 空间中的连续小波？下一篇在频域中使用互相关计算延迟/滞后