Matlab 中 arburg 函数的理解与实际使用方法

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. 向后预测