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是方差。
本文介绍了Levinson递归算法在信号处理中的应用,通过计算自相关函数和使用滤波器,实现了噪声估计和功率谱分析。通过实例展示了如何利用该算法处理随机信号,并估计频率响应和关键频点。

被折叠的 条评论
为什么被折叠?



