自回归（AR）模型的功率谱估计（实现）

上一部分介绍了AR模型的理论知识，这一部分将介绍AR模型的各种估计方法。点击这里，快速查看理论知识。在这些实验中，均假设N=1024，P=512。

Yule-Walker法

流程图

python 代码

import scipy
import numpy as np
from matplotlib import pyplot as plt


# 自相关函数
def x_corr(x):
    x_inv = x[::-1]
    # r_x = np.correlate(x, x)
    r_x = np.convolve(x, x_inv)
    r_x /= len(x)
    return r_x


if __name__ == "__main__":
    N = int(input("采样点数N："))
    P = int(input("模型阶数P："))
    # a_N [0, 1, ..., N-1]
    a_N = np.arange(0, N, 1)

    # 信号频率
    f1 = 0.1
    f2 = 0.13
    a_X = 2 * np.sin(2 * np.pi * f1 * a_N + np.pi / 3) + 10 * np.sin(2 * np.pi * f2 * a_N + np.pi / 4)

    # 自相关函数Rxx
    # Rxx(0), Rxx(1), ..., Rxx(N-1)
    Rxx = x_corr(a_X)[N-1:]
    # 取前P+1个
    # Rxx(0), Rxx(1), Rxx(P-1), Rxx(P)
    Rxx = Rxx[: P+1]

    # 自相关矩阵a_Rxx
    # [Rxx(0)   Rxx(1)   ... Rxx(p-1)]
    # [Rxx(-1)  Rxx(0)   ... Rxx(p-2)]
    # ...
    # [Rxx(1-p) Rxx(2-p) ... Rxx(0)  ]
    a_Rxx = np.zeros([P, P])
    for i in range(P):
        for j in range(P):
            a_Rxx[i, j] = Rxx[abs(i-j)]

    # 结果矩阵a_rx
    # [-Rxx(1)]
    # [-Rxx(2)]
    # ...
    # [-Rxx(p)]
    a_rx = np.zeros([P, 1])
    for i in range(P):
        a_rx[i, 0] = -Rxx[i + 1]

    # 系数矩阵a_apk
    # [a_p1]
    # [a_p2]
    # ...
    # [a_pp]
    a_Rxx_inv = np.linalg.inv(a_Rxx)
    a_apk = np.matmul(a_Rxx_inv, a_rx)

    # 标准差sigma
    sigma = Rxx[0]
    for i in range(P):
        sigma += a_apk[i, 0] * Rxx[i + 1]
    sigma = np.sqrt(sigma)

    # [ap1, ..., app]->[1, ap1, ..., app]
    a_ap = np.insert(a_apk, 0, [1])
    w, h = scipy.signal.freqz(sigma, a_ap, worN=N)

    Hf = abs(h)
    psd = Hf**2
    f = w / (2 * np.pi)

    plt.figure(1)
    plt.subplot(211)
    plt.plot(f, Hf)
    plt.xlabel("f/hz")
    plt.ylabel("Hf")

    plt.subplot(212)
    plt.plot(f, psd)
    plt.xlabel("f/hz")
    plt.ylabel("psd")
    plt.show()

实验结果

Levinson-Durbin法

流程图

python代码

import numpy as np
import scipy
from matplotlib import pyplot as plt


def x_corr(a_x):
    a_x_inv = a_x[::-1]
    a_rx = np.convolve(a_x, a_x_inv)
    a_rx /= len(a_x)
    return a_rx


if __name__ == "__main__":
    N = int(input("采样点数N："))
    P = int(input("模型阶数P："))

    a_N = np.arange(N)
    # 信号频率
    f1 = 0.1
    f2 = 0.3
    # 信号
    a_X = 2 * np.sin(2 * np.pi * f1 * a_N + np.pi / 3) + 10 * np.sin(2 * np.pi * f2 * a_N + np.pi / 4)

    # 自相关函数Rx
    Rxx = x_corr(a_X)[N - 1:]
    Rxx = Rxx[: P + 1]

    # 初始化参数apk, sigma
    # p=1
    a_apk = np.zeros([P + 1, P + 1])
    a_apk[1, 1] = -Rxx[1] / Rxx[0]
    a_sigma = np.zeros(P + 1)
    a_sigma[1] = (1 - np.square(a_apk[1, 1])) * Rxx[0]

    # p>1
    for p in range(2, P + 1):
        kp_sum = 0
        for i in range(1, p):
            kp_sum += a_apk[p - 1, i] * Rxx[p - i]
        a_apk[p, p] = - (Rxx[p] + kp_sum) / a_sigma[p - 1]
        for k in range(1, p):
            a_apk[p, k] = a_apk[p - 1, k] + a_apk[p, p] * a_apk[p - 1, p - k]
        a_sigma[p] = (1 - np.square(a_apk[p, p])) * a_sigma[p - 1]
    
    a_ap = a_apk[- 1][1:]

    # a_apk = np.zeros([P, P])
    # a_apk[0, 0] = -Rxx[1] / Rxx[0]
    # a_sigma = np.zeros(P)
    # a_sigma[0] = (1 - np.square(a_apk[0, 0])) * Rxx[0]
    #
    # for p in range(1, P):
    #     kp_sum = 0
    #     for i in range(0, p):
    #         kp_sum += a_apk[p - 1, i] * Rxx[p - i]
    #     a_apk[p, p] = - (Rxx[p + 1] + kp_sum) / a_sigma[p - 1]
    #     for k in range(p):
    #         a_apk[p, k] = a_apk[p - 1, k] + a_apk[p, p] * a_apk[p - 1, p - k - 1]
    #     a_sigma[p] = a_sigma[p - 1] * (1 - np.square(a_apk[p, p]))

    # a_ap = a_apk[- 1]
    
    sigma = a_sigma[- 1]
    sigma = np.sqrt(sigma)
    a_ap = np.insert(a_ap, 0, [1])
    w, h = scipy.signal.freqz(sigma, a_ap, worN=N)

    Hf = abs(h)
    psd = Hf**2
    f = w / (2 * np.pi)

    plt.figure(1)
    plt.subplot(211)
    plt.plot(f, Hf)
    plt.xlabel("f/hz")
    plt.ylabel("Hf")

    plt.subplot(212)
    plt.plot(f, psd)
    plt.xlabel("f/hz")
    plt.ylabel("psd")
    plt.show()

实验结果

Burg法

流程图

python代码

import numpy as np
import scipy
from matplotlib import pyplot as plt


def x_corr(a_x):
    a_x_inv = a_x[::-1]
    a_rx = np.convolve(a_x, a_x_inv)
    a_rx /= len(a_x)
    return a_rx


if __name__ == "__main__":
    N = int(input("采样点数N："))  # 1024
    P = int(input("模型阶数P："))  # 512
    a_N = np.arange(N)
    # 信号频率
    f1 = 0.1
    f2 = 0.13
    a_X = 2 * np.sin(2 * np.pi * f1 * a_N + np.pi / 3) + 10 * np.sin(2 * np.pi * f2 * a_N + np.pi / 4)

    a_epn = np.zeros([P+1, N])
    a_bpn = np.zeros([P+1, N])

    # 初始化e0n, b0n
    a_epn[0] = a_X
    a_bpn[0] = a_X
    # for i in range(N):
    #     a_epn[0, i] = a_X[i]
    #     a_bpn[0, i] = a_X[i]

    a_Kp = np.zeros(P+1)
    a_apk = np.zeros([P+1, P+1])
    a_sigma = np.zeros(P+1)

    # 初始化sigma1
    # a_sigma[0] = x_corr(a_X)[N-1:][0]
    sum1 = 0
    for n in range(N):
        sum1 += np.square(a_X[n])
    a_sigma[0] = sum1 / N

    for p in range(1, P+1):
        Kp_sum1 = 0
        Kp_sum2 = 0
        for n in range(p, N):
            Kp_sum1 += a_epn[p-1, n] * a_bpn[p-1, n-1]
            Kp_sum2 += np.square(a_epn[p-1, n]) + np.square(a_bpn[p-1, n-1])
        a_Kp[p] = -(2 * Kp_sum1) / Kp_sum2

        for n in range(p, N):
            a_epn[p, n] = a_epn[p-1, n] + a_Kp[p] * a_bpn[p-1, n-1]
            a_bpn[p, n] = a_bpn[p-1, n-1] + a_Kp[p] * a_epn[p-1, n]

    a_apk[1, 1] = a_Kp[1]
    a_sigma[1] = (1 - np.square(a_Kp[1])) * a_sigma[0]
    for p in range(2, P+1):
        a_apk[p, p] = a_Kp[p]
        for i in range(1, p):
            a_apk[p, i] = a_apk[p-1, i] + a_apk[p, p] * a_apk[p-1, p-i]
        a_sigma[p] = (1 - np.square(a_Kp[p])) * a_sigma[p - 1]

    a_ap = a_apk[-1][1:]
    sigma = np.sqrt(a_sigma[-1])
    a_ap = np.insert(a_ap, 0, [1])
    w, h = scipy.signal.freqz(sigma, a_ap, worN=N)
    Hf = abs(h)
    psd = Hf ** 2
    f = w / (2 * np.pi)

    plt.figure(1)
    plt.subplot(211)
    plt.plot(f, Hf)
    plt.xlabel("f/hz")
    plt.ylabel("Hf")

    plt.subplot(212)
    plt.plot(f, psd)
    plt.xlabel("f/hz")
    plt.ylabel("psd")
    plt.show()

实验结果

协方差法

python代码

个人感觉代码有点问题，但是又不知道哪里有问题，麻烦各位老哥知道的赐教一下。

import numpy as np
import scipy
from matplotlib import pyplot as plt


if __name__ == "__main__":
    N = int(input("采样点数N："))
    P = int(input("模型阶数P："))
    a_N = np.arange(N)
    # 信号频率
    f1 = 0.1
    f2 = 0.13
    a_X = 2 * np.sin(2 * np.pi * f1 * a_N + np.pi / 3) + 10 * np.sin(2 * np.pi * f2 * a_N + np.pi / 4)

    Cxx = np.zeros([P+1, P+1])
    for i in range(P+1):
        for k in range(P+1):
            sum_c = 0
            for n in range(P, N):
                sum_c += a_X[n-k] * a_X[n-i]
            Cxx[i, k] = sum_c / (N - P)

    a_Cxx = Cxx[1:, 1:]
    a_Cxx_inv = np.linalg.inv(a_Cxx)

    a_cx = Cxx[:, 0][1:]
    a_apk = np.matmul(a_Cxx_inv, -a_cx[:, None])

    sigma = Cxx[0, 0]
    for k in range(P):
        sigma += a_apk[k, 0] * Cxx[0, k + 1]
    sigma = np.sqrt(sigma)

    # [ap1, ..., app]->[1, ap1, ..., app]
    a_ap = np.insert(a_apk, 0, [1])
    w, h = scipy.signal.freqz(1, a_ap, worN=N)

    Hf = abs(h)
    psd = Hf**2
    f = w / (2 * np.pi)

    plt.figure(1)
    plt.subplot(211)
    plt.plot(f, Hf)
    plt.xlabel("f/hz")
    plt.ylabel("Hf")

    plt.subplot(212)
    plt.plot(f, psd)
    plt.xlabel("f/hz")
    plt.ylabel("psd")
    plt.show()

实验结果

修正协方差法

流程图

python代码

这个代码也感觉有问题，当设置p=512,N=1024时，程序一直报错。没办法只有改成p=512才能正常运行。有会的兄弟姐妹私信我一下，谢谢啦！

import numpy as np
import scipy
from matplotlib import pyplot as plt


if __name__ == "__main__":
    N = int(input("采样点数N："))
    P = int(input("模型阶数P："))
    a_N = np.arange(N)
    # 信号频率
    f1 = 0.1
    f2 = 0.13
    a_X = 2 * np.sin(2 * np.pi * f1 * a_N + np.pi / 3) + 10 * np.sin(2 * np.pi * f2 * a_N + np.pi / 4)

    Cxx = np.zeros([P+1, P+1])
    for i in range(P+1):
        for k in range(P+1):
            sum1 = 0
            sum2 = 0
            for n in range(P, N):
                sum1 += a_X[n-k] * a_X[n-i]
            for m in range(0, N-P):
                sum2 += a_X[m+k] * a_X[m+i]
            Cxx[i, k] = (sum1 + sum2) / (2 * (N - P))

    a_Cxx = Cxx[1:, 1:]
    a_Cxx_inv = np.linalg.inv(a_Cxx)

    a_cx = Cxx[:, 0][1:]
    a_apk = np.matmul(a_Cxx_inv, -a_cx[:, None])

    sigma = Cxx[0, 0]
    for k in range(P):
        sigma += a_apk[k, 0] * Cxx[0, k + 1]
    sigma = np.sqrt(sigma)

    # [ap1, ..., app]->[1, ap1, ..., app]
    a_ap = np.insert(a_apk, 0, [1])
    w, h = scipy.signal.freqz(sigma, a_ap, worN=N)

    Hf = abs(h)
    psd = Hf**2
    f = w / (2 * np.pi)

    plt.figure(1)
    plt.subplot(211)
    plt.plot(f, Hf)
    plt.xlabel("f/hz")
    plt.ylabel("Hf")

    plt.subplot(212)
    plt.plot(f, psd)
    plt.xlabel("f/hz")
    plt.ylabel("psd")
    plt.show()

实验结果