python:levinson-durbin算法

import numpy as np
import random
import scipy.signal
import matplotlib.pyplot as plt


class Levinson:
    def __init__(self):
        super(Levinson, self).__init__()
        return

    def forward(self, corr, P):
        p = P
        a = np.zeros([p + 1, p + 1])
        rou = np.zeros(p + 1)
        rou[0] = corr[0]
        a[1][1] = -corr[1] / rou[0]
        rou[1] = rou[0] * (1 - (a[1][1] ** 2))

        for i in range(2, p + 1):
            sum = 0.0
            for m in range(1, i):
                sum = sum + a[i - 1][m] * corr[i - m]
            a[i][i] = -(corr[i] + sum) / rou[i - 1]
            for j in range(1, i):
                a[i][j] = a[i - 1][j] + a[i][i] * a[i - 1][i - j]
            rou[i] = rou[i - 1] * (1 - (a[i][i]) ** 2)

        return a, rou

    def autocorrelation(self, x):
        n = len(x)
        y = np.array(x)
        for m in range(n):
            sum = 0.0
            for i in range(n - m):
                sum = sum + x[i] * x[i + m]
            y[m] = sum / n
        return y


def get_noise(n, m, s):
    x = np.zeros(n)
    for i in range(n):
        x[i] = random.gauss(m, s)
    return x


if __name__ == '__main__':
    N = 256
    p = 90
    mu = 0
    sigma = 1
    random.seed(12)
    x = get_noise(N, mu, sigma)
    Levinson = Levinson()
    corr = Levinson.autocorrelation(x)
    a, rou = Levinson.forward(corr, p)

    G2 = rou[p]
    ap = []
    for k in range(1, p + 1):
        ap.append(a[p][k])
    apk = np.insert(ap, 0, [1])
    G = np.sqrt(G2)

    w, h = scipy.signal.freqz(G, apk, worN=N)
    Hf = abs(h)
    Sx = Hf ** 2
    f = w / (2 * np.pi)

    f1 = scipy.signal.argrelextrema(Sx, np.greater)[0]
    f2 = f1 / (2 * N)
    print('f的估计值:', f2)

    fig = plt.figure(figsize=(8, 8), dpi=100)
    fig.suptitle('N = {} P = {}'.format(N, p), fontsize=30)
    plt.subplots_adjust(hspace=0.35)
    plt.subplot(2, 1, 1)
    plt.plot(x)
    plt.title('random signal', fontsize=20)
    plt.xlabel('n', fontsize=15)
    plt.ylabel('x', fontsize=15)

    plt.subplot(2, 1, 2)
    plt.plot(f, Sx)
    Sy = Sx[f1]
    print('power的估计值：',Sy)
    plt.scatter(f2, Sy, marker='o', edgecolors=['r'], color='none')
    plt.title('power specturm', fontsize=20)
    plt.xlabel('f', fontsize=15)
    plt.ylabel('Sx', fontsize=15)
    plt.show()

N是序列长度。P是滤波器阶数，mu是白噪声均值，sigma是方差。