Matlab 中 arburg 函数的理解与实际使用方法
阅读原文时间:2023年07月09日阅读:1

1. 理解

  1.1 Matlab 帮助:

  a = arburg(x,p)返回与输入数组x的p阶模型相对应的归一化自回归(AR)参数。 如果x是一个向量,则输出数组a是一个行向量。 如果x是矩阵,则参数沿模型的第n行位于x的第n列。 a有p + 1列。 p必须小于x的元素(或行)数。

  [a,e] = arburg(x,p)返回白噪声输入的估计方差e。

  [a,e,rc] = arburg(x,p)返回rc中的反射系数。

  1.2 自我的一些理解

AR(P) 模型作用: 使用前几个数据来预测后面的数据,p是阶数,对应p个历史数据

  In an AR model of order p, the current output is a linear combination of the past p outputs plus a white noise input.

  The weights on the p past outputs minimize the mean-square prediction error of the autoregression.

  If y(n) is the current value of the output and x(n) is a zero mean white noise input, the AR(p) model is:

  在阶数为p的AR模型中,当前输出是过去p输出的线性组合加上白噪声

  p过去的输出上的权重使自回归的均方预测误差最小。

  如果y(n)是输出的当前值并且x(n)是零平均白噪声输入,则AR(p)模型为:

  y(n) =  -(负号)  p个累加a(k)*y(n-k) + 0均值的高斯白噪声x(n)  [通常其他表达里使用 w(n)]

  -----------------------------------

  [ ar_coeffs,NoiseVariance] = arburg(data,order)

  ar_coeffs 第一位不需要用,是用于归一化,后面是对应的系数

  NoiseVariance 是方差

  举例:

  p=3: [a E]=arburg(x,3)

  a =1.0000 -0.6982 -0.2626 0.0739

  E =0.4567

  p=12: [a E]=arburg(x,12)

  a =1.0000 -0.6495 -0.3066 -0.0934 0.0987 0.4076 -0.1786 -0.0126 -0.0805 -0.0899 0.0382 0.1628 -0.2501

  E =0.3237

  1.3 网络

https://blog.csdn.net/weixin_43165881/article/details/106878784

  1.4 python

  Python 中 librosa 库中有 lpc 函数是使用的 burg 方法。

2使用:

  向前向后预测

  寻找一个案列,进行AR拟合,与原函数一致

  这个AR(1)模型在当初陀螺仪卡尔曼滤波中用过,w(n) 是高斯白噪声,先要通过一组数据算出方差。

  burg的具体的应用方法还需要再组织一下想法

  2021-03-25 日更新

  恰巧这日用python 实现了下, python种的librosa 库种lpc 函数是使用 burg法的线性预测 ,代码如下

# %% 0_0.import libs
import numpy as np
import librosa as lb
import matplotlib.pyplot as plt

#%% 0_1.def funcs
def next_predict(sig,M_nexts,p):

a= lb.lpc(sig,p) # 计算模型系数

len\_sig = len(sig)  
len\_sig\_pred = M\_nexts+ sig  
sig\_predict = np.zeros(len\_sig\_pred)  
sig\_predict\[:len\_sig\]=sig

for i in range(M\_nexts):  
    sig\_predict\[ len\_sig + i \] = -np.dot(a\[1:\],sig\_predict\[len\_sig + i :len\_sig+ i-p:-1\]) # + error  

sig\_predict\_out = sig\_predict\[len\_sig:\]

return sig\_predict\_out

def forward_predict(sig,M_forward,p):

a= lb.lpc(sig,p) # 计算模型系数

len\_sig = len(sig)  
len\_sig\_pred = M\_forward + len\_sig  
sig\_predict = np.zeros(len\_sig\_pred)  
sig\_predict\[M\_forward:\]=sig

for i in range(M\_forward):  
    sig\_predict\[ M\_forward -1 - i \] = - np.dot(a\[1:\],sig\_predict\[ M\_forward-i : M\_forward-i+ p\]) # + error  

sig\_predict\_out = sig\_predict\[:M\_forward\]  

return sig\_predict\_out

# %% 1. 数据的ar-burg 预测系数

t = np.arange(0,10,0.01)
sig1 = np.sin(2*np.pi*t) # sin 函数
sig2 = 0.2 *np.sin(10*np.pi*t)
noise = np.random.randn(len(sig1))
sig = sig1 + sig2 + noise

%% 2. 求arburg 预测的数据

ar,err = arburg(sig,n_Order)

librosa.lpc 使用了burg法计算 lpc系数

p=300
a= lb.lpc(sig,p)

%% 根据得到的系数进行向前向后预测

M_nexts = 200 # 向后预测的个数
M_befores = 200 # 向后预测的个数

设第 n个数的预测为,预测M个数

sig_predict(n)= - a[1]*sig(n-1) -a[2]*sig(n-2) - a[3]*sig(n-3)… -a[p]*sig(n-p)

len_sig = len(sig)
len_sig_pred = M_nexts+ len_sig
sig_predict = np.zeros(len_sig_pred)
sig_predict[:len(sig)]=sig

预测

for i in range(M_nexts):
sig_predict[ len_sig + i ] = -np.dot(a[1:],sig_predict[len_sig + i :len_sig+ i-p:-1]) # + error

如果我们知道err,计算出err的方差,也可以进行预测

plt.figure()
plt.plot(sig_predict)
plt.plot(sig )

sig_f_pred = forward_predict(sig,M_befores,p)

sig_all = np.hstack([sig_f_pred,sig])
plt.figure()
plt.plot(sig_all)
plt.plot(sig_f_pred)

  

图1.向前预测

图2. 向后预测

手机扫一扫

移动阅读更方便

阿里云服务器
腾讯云服务器
七牛云服务器

你可能感兴趣的文章