上一部分介绍了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才能正常运行。有会的兄弟姐妹私信我一下，谢谢啦！文章来源：https://www.toymoban.com/news/detail-531950.html
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()
实验结果
文章来源地址https://www.toymoban.com/news/detail-531950.html
到了这里，关于自回归（AR）模型的功率谱估计（实现）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
