诚心求教Matlab关于用Yule-Walker求解AR模型参数问题

Question

如题所述

Answer 1

function ARMODEL()
Fs = 1000;
t = 0:1/Fs:15;
N = size(t,2) %数据样值点数
randn('state',0);
x = cos(2*pi*t*200) + randn(1,N); % 200Hz cosine plus noise

%计算N个取样数据的取样数据自相关函数
rxx = zeros(1,N); %保存取样数据自相关函数的变量
for m = 0:N-1
sum = 0;
for n= 1:N-m
temp1 = x(n)*x(m+n);
sum = sum + temp1;
end
rxx(m+1) = sum/N;
end

%采用Levison-Durbin算法求解AR模型的Yule-Walker模型
%需要确定AR模型理论公式中的参数：白噪声w(n)的方差、方程系数a1……ap（这里包括了模型的阶次）
PMAX = 100; %设定AR模型最高阶次
atemp1 = zeros(1,PMAX+1); %保存方程系数的中间变量
atemp2 = zeros(1,PMAX+1); %保存方程系数的中间变量
deviationtemp1 = zeros; %保存白噪声w(n)方差的中间变量
deviationtemp2 = zeros; %保存白噪声w(n)方差的中间变量

%AR(1)模型：x(n) + a1*x(n-1) = w(n)
%其Yule-Walker方程： R(0)*1 + R(1)*a1 = deviation1;
% R(1)*1 + R(0)*a1 = 0;
%求解方程确定a1、deviation1
atemp1(1) = 1;
atemp1(2) = -rxx(2)/rxx(1);
atemp2 = atemp1;
deviationtemp1 = ( rxx(1)*rxx(1) - rxx(2)*rxx(2) )/rxx(1);
deviationtemp2 = deviationtemp1;

%利用Levison-Durbin迭代算法计算AR模型参数
%根据FPE准则、AIC准则和BIC准则确定AR模型的阶次
%atemp1、deviation1保存第k次的运算结果
%atemp2、deviation2保存第k+1次的运算结果
FPE(1) = deviationtemp1*(N+2)/N;
AIC(1) = log(deviationtemp1) + 2/N;
BIC(1) = log(deviationtemp1) + log(N)/N;
veriance(1) = deviationtemp1;

criteria = 3

for P = 2:PMAX
sum1 = 0;
sum2 = 0;
for i = 2:(P+1)
sum1 = atemp1(i)*rxx(i) + sum1;
end

for i = 1:(P+1)
sum2 = atemp1(i)*rxx(P+2-i) + sum2;
end

deviationtemp1 = rxx(1) + sum1;
dk = sum2;
ref(P) = dk/deviationtemp1;
deviationtemp2 = ( 1 - ref(P)*ref(P) )*deviationtemp1;

for i = 2:(P+1)
atemp2(i) = atemp1(i) - ref(P)*atemp1(P+2-i);
end
%计算AR(P)模型参数
atemp1 = atemp2;
veriance(P) = deviationtemp2本回答被提问者和网友采纳