AR模型

AR模型的系统函数H(z)可以表示为：

我们的目的就是要求解系统函数的参数a和增益G。

Yule_Walker方程

矩阵形式

根据生成的矩阵，可以解出p个参数，再根据自相关函数，可以求出系统增益G。
Yule_Walker方程也可以写成：

Levinson-Durbin快速递推法

可以看出通过矩阵求逆解Yule-Walker方程的运算量很大，通过观察可以看出Yule_Walker方程的自相关矩阵有一系列良好的性质，可以通过递推的方法求p个参数，从而避免矩阵求逆的运算，减少运算量。
令p=1，可以得到1阶AR模型Yule-Walker方程：

同理，令p=2，3，4…，将不断变化的p用m表示，可以得到m阶AR模型参数Levinson-Durbin递推算法：

MATLAB实现

Yule_Walker方程

clear
clc
% N = 128;
% p = 64;
N = input('请选择采样点数N:');
p = input('请选择AR模型阶数p（建议N/3<p<N/2):');
n = 0:N-1;
dt = 1;t = n*dt;
xn = 2^0.5*sin(2*pi*0.1*t+pi/3) + 10*sin(2*pi*0.13*t+pi/4); %+ randn(1,N);
R = xcorr(xn,'biased');% 求xn的自相关函数
Rx = zeros(1,p+1);% 预分配内存
% 取对称序列的后半部分
for i = 1:p+1Rx(i) = R(i+N-1);
end
Rx
Rmatl = zeros(p,p);
% 生成p个方程对应的矩阵式
for m = 1:pfor n = 1:pRmatl(m,n) = Rx(max(m,n)-min(m,n)+1);end
end
Rmatl
Rmatr = zeros(p,1);
for m = 1:pRmatr(m,1) = -Rx(m+1);
end
% 输出参数ai
ai = (inv(Rmatl)*Rmatr)';% (Rmatl\Rmatr)'
% 根据自相关函数和ai求解系统增益
G = Rx(1);
for i = 1:pG = G+ai(i)*Rx(i);
end
fprintf('系统增益G=%f\n',G^0.5);
[H,w] = freqz(G^0.5,[1,ai],N);% 在2*pi范围内N个等分点求系统函数
% Px = (G^0.5)*(abs(H)).^2;
Hf = abs(H);
f = w/(2*pi);
plot(w/(2*pi),Hf);
[magsor,fsor] = findpeaks(Hf,f,'SortStr','descend');% 降序排列
for i = 1:2fprintf('预测频率f%d=%f\n',i,fsor(i));
end

实验结果

Levinson-Durbin快速递推法

clear
clc
% N = 8;
% p = 4;
N = input('请选择采样点数N:');
p = input('请选择AR模型阶数p（建议N/3<p<N/2):');
n = 0:N-1;
dt = 1;t = n*dt;
xn = 2^0.5*sin(2*pi*0.1*t+pi/3) + 10*sin(2*pi*0.3*t+pi/4);% + randn(1,N);
R = xcorr(xn,'biased');% 求xn的自相关函数
Rx = zeros(1,p+1);% 预分配内存
% 取对称序列的后半部分
for i = 1:p+1Rx(i) = R(i+N-1);
end
% 求出1阶AR模型的系数和预测误差功率
a(1,1) = -Rx(2)/Rx(1);
rou(1) = Rx(1)*(1-(a(1,1))^2);
for m = 2:pkmsum = 0;for i = 1:m-1kmsum = kmsum+a(m-1,i)*Rx(m-i+1);endk(m) = -(Rx(m+1)+kmsum)/rou(m-1);a(m,m) = k(m);for i = 1:m-1a(m,i) = a(m-1,i)+k(m)*a(m-1,m-i);end;rou(m) = rou(m-1)*(1-(k(m))^2);
end
% rou
% plot(rou)
xlabel('p');ylabel('rou');
G = (rou(p))^2;
[H,w] = freqz(G,[1,a(p,1:p)],N);% 在2*pi范围内N个等分点求系统函数
% Px = (G^0.5)*(abs(H)).^2;
Hf = abs(H);
f = w/(2*pi);
plot(w/(2*pi),Hf);
[magsor,fsor] = findpeaks(Hf,f,'SortStr','descend');% 降序排列
for i = 1:2fprintf('预测频率f%d=%f\n',i,fsor(i));
end

实验结果

Python实现

Yule_Walker方程

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
import time# 定义了寻找峰值的函数，返回值是元素在数组中的位置
def findpeaks(data):length = len(data)weizhi = []for i in range(1,length-1):if data[i] > data[i-1] and data[i] > data[i+1]:weizhi.append(i)return weizhi# 定义了自相关函数
def xcorr(data):length = len(data)Rx = []datareverse = data[::-1] # 将数组data从后向前取值，每次步进值为1Rx = scipy.signal.convolve(data,datareverse)/length # 自相关函数return RxN = input('请选择采样点数N:') # input输入的N是str类型
p = input('请选择AR模型阶数p（建议N/3<p<N/2):') # 输入阶数
N = int(N) # 把N转换成int类形
p = int(p)
dt = 1
t = np.arange(0,N,dt)# 数组是从0开始的，t=[0 1 2 ... N-1]# 信号频率
f1 = 0.1
f2 = 0.13
x = np.sqrt(2)*np.sin(2*np.pi*f1*t+(np.pi)/3) + 10*np.sin(2*np.pi*f2*t+(np.pi)/4)# 输入信号
#print(x)
Rx1 = xcorr(x)# 求自相关
Rx = Rx1[N-1:N+p]
#print('Rx=',Rx)#start = time.clock()
#start = time.perf_counter()
# 生成自相关系数矩阵
Rmatl = np.ones([p,p])
for i in range(0,p):for j in range(0,p):Rmatl[i][j] = Rx[max(i,j)-min(i,j)]
#print('Rmatl=',Rmatl)Rmatr = np.ones([p,1])
for i in range(0,p):Rmatr[i][0] = -Rx[i+1]
#print('Rmatr=',Rmatr)# 求系数ai
ai = np.ones([p,1])
Rmatl_inv = np.linalg.inv(Rmatl) # 对左面自相关系数矩阵求逆
ai = np.dot(Rmatl_inv,Rmatr) # 得到各个ai的值
aiT = ai.T # 求ai的转置
#print(aiT)
#end = time.clock()
#end = time.perf_counter()
#print('运行时间=',end-start)# 求系统增益G
G1 = Rx[0]
for i in range(1,p):G1 = G1 + aiT[0][i]*Rx[i]
G = np.sqrt(G1)
#print('系统增益G=',G)a = np.insert(aiT,0,[1])# 系统函数分母的系数
w,h  = scipy.signal.freqz(G,a,worN = N)
Hf = abs(h)
f = w/(2*np.pi)plt.plot(f,Hf)
plt.show()

实验结果

Levinson-Durbin快速递推法

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
import time# 定义了寻找峰值的函数，返回值是元素在数组中的位置
def findpeaks(data):length = len(data)weizhi = []for i in range(1,length-1):if data[i] > data[i-1] and data[i] > data[i+1]:weizhi.append(i)return weizhi# 定义了自相关函数
def xcorr(data):length = len(data)Rx = []datareverse = data[::-1] # 将数组data从后向前取值，每次步进值为1Rx = scipy.signal.convolve(data,datareverse)/length # 自相关函数return RxN = input('请选择采样点数N:') # input输入的N是str类型
p = input('请选择AR模型阶数p（建议N/3<p<N/2):') # 输入阶数
N = int(N) # 把N转换成int类形
p = int(p)
dt = 1
t = np.arange(0,N,dt)# 数组是从0开始的，t=[0 1 2 ... N-1]# 信号频率
f1 = 0.1
f2 = 0.3
# 输入信号
x = np.sqrt(2)*np.sin(2*np.pi*f1*t+(np.pi)/3) + 10*np.sin(2*np.pi*f2*t+(np.pi)/4)
#print(x)
Rx1 = xcorr(x)# 求自相关
Rx = Rx1[N-1:N+p]
#print('Rx=',Rx)#start = time.perf_counter()
#设定初值
a = np.zeros([p,p])
a[0][0] = -Rx[1]/Rx[0]
rou = np.zeros(p)
rou[0] = Rx[0]*(1-np.square(a[0][0]))
k = np.zeros(p)#p阶AR模型参数Levinson_Durbin递推算法
for m in range(1,p):kmsum = 0for i in range(0,m):kmsum = kmsum+a[m-1][i]*Rx[m-i]k[m] = -(Rx[m+1]+kmsum)/rou[m-1]a[m][m] = k[m]for i in range(0,m):a[m][i] = a[m-1][i]+k[m]*a[m-1][m-i-1]rou[m] = rou[m-1]*(1-np.square(k[m]))
#end = time.perf_counter()
#print('运行时间=',end-start)G = np.sqrt(rou[p-1])
a = np.insert(a[p-1],0,[1])# 系统函数分母的系数
w,h  = scipy.signal.freqz(G,a,worN = N)
Hf = abs(h)
f = w/(2*np.pi)plt.plot(f,Hf)
plt.show()

实验结果

随机信号处理AR模型Yule_Walker方程直接解法和Levinson_Durbin递推法的MATLAB与Python实现相关推荐

建模股票价格数据并进行预测（统计信号模型）：随机信号AR模型+Yule-Walker方程_Python...
1.背景: 针对股票市场中AR 模型的识别.建立和估计问题,利用AR 模型算法对股票价格进行预测. 2.模型选取: 股票的价格可视为随机信号,将此随机信号建模为:一个白噪声通过LTI系统的输出,通过原 ...
python财务报表预测股票价格_建模股票价格数据并进行预测（统计信号模型）：随机信号AR模型+Yule-Walker方程_Python...
1.背景: 针对股票市场中AR 模型的识别.建立和估计问题,利用AR 模型算法对股票价格进行预测. 2.模型选取: 股票的价格可视为随机信号,将此随机信号建模为:一个白噪声通过LTI系统的输出,通过原 ...
现代信号处理——AR模型谱估计
AR谱估计方法可归结为求解AR模型系数或线性预测器系数的问题. AR模型参数估计方法:信号预测误差最小原则(或预测误差功率最小) 自相关法(Levison递推法) Burg法协方差法修正协方差法( ...
随机数字信号处理实验报告三——Levinson和Burg递推法MATLAB实现
完整的实验报告下载连接https://download.csdn.net/download/LIsaWinLee/14884452 一.实验原理随机信号的功率谱密度用来描述信号的能量特征随频率的变化 ...
matlab pburg,现代数字信号处理——AR模型
1. AR模型概念观 AR模型是一种线性预测,即已知N个数据,可由模型推出第N点前面或后面的数据(设推出P点),所以其本质类似于插值,其目的都是为了增加有效数据,只是AR模型是由N点递推,而插值是由两 ...
现代数字信号处理——AR模型
本文目标:分析AR模型并求解AR模型的输出x(n)的功率谱. 1. AR模型概念观数字信号处理功率谱估计方法分经典功率谱估计和现代功率谱估计,现代功率谱估计以参数模型功率谱估计为代表,参数功率谱模型 ...
1308 方程的解（组合计数--隔板法）
1. 问题描述: 佳佳碰到了一个难题,请你来帮忙解决.对于不定方程 a1+a2+⋯+ak−1+ak=g(x),其中 k ≥ 1 且 k∈N∗,x 是正整数,g(x)=x ^ x mod 1000(即 ...
随机信号的参数建模法--AR模型及Matlab实现
目录一.随机信号参数模型二.AR模型三.AR模型参数的估计 1. AR 模型参数和自相关函数的关系 2. Y-W 方程的解法--L-D 算法 2-1 前向预测器 2-2 建立更高阶的AR模型 2 ...
随机信号的参数建模法( AR 模型)
引言为随机信号建立参数模型是研究随机信号的一种基本方法,其含义是认为随机信号x(n)x(n)x(n)是由白噪w(n)w(n)w(n)激励某一确定系统的响应(如图7.5).只要白噪的参数确定了,研究随 ...

随机信号处理AR模型Yule_Walker方程直接解法和Levinson_Durbin递推法的MATLAB与Python实现

AR模型

Yule_Walker方程

Levinson-Durbin快速递推法

MATLAB实现

Yule_Walker方程

实验结果

Levinson-Durbin快速递推法

实验结果

Python实现

Yule_Walker方程

实验结果

Levinson-Durbin快速递推法

实验结果

随机信号处理AR模型Yule_Walker方程直接解法和Levinson_Durbin递推法的MATLAB与Python实现相关推荐

最新文章

热门文章