%自編函數實現AR模型的最小二乘估計
%-----------------------------------------
close all;
clear all;
clc;
N=128;
n=1:128;
a1=sqrt(20);
a2=sqrt(2);
w=randn(128,1);
x=a1*sin(2*pi*0.2*n)+a2*sin(2*pi*0.213*n)+w';%數據源序列
%xcorr
xxcorr=xcorr(x,'unbiased');

%-----------------------------------------
%ls
% Matrix---R4x4階自相關矩陣
p=4;%AR階數4
q=20;%q>=p
%matrix
for i=1:p
    for j=1:p
        recorr(i,j)=xxcorr(q+i+j-1-p);%生成自相關序列
    end
end
disp('自相關序列的維數=');%4維
recorrlegnth=size(recorr);
%r

for i=1:p %i=1,2,3,4
    r(i)=xxcorr(q+i);
end
disp('r的維數=');
rlength=size(r)
r=-r';%因為Ra=-r
yyy=(inv(recorr'*recorr)*recorr')*r;
disp('求解序列的維數');
yylength=length(yyy)
for i=1:yylength
    %--------------a=[ap,...,a1]
    yyy1(i)=yyy(yylength+1-i);%序列反轉
    %--------------a=[a1,...,ap]
end
figure(1);
freqz(1,yyy1,1024,1);
title('自編函數實現AR模型的最小二乘估計');
legend('AR階數=4');
grid on;