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

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