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

时间:2024-03-07 17:30:08

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