一、最小二乘法

clear
clc
fs=512;
m=0:(fs-1);
xn=sqrt(20)*sin(2*pi*0.2*m)+sqrt(2)*sin(2*pi*0.213*m)+randn(1,fs); %产生含有噪声的序列xn
Rx=xcorr(xn);
q=200;
%以下p1、p2分别是AR的阶数
p1=4;
p2=6;
%p1=4时：
p=p1;
R=zeros(p,p);
r=zeros(1,p);
for i=1:p
    for j=1:p
        R(i,j)=R(i,j)+Rx(1,fs+q-p+i+j-1);
    end
end
for m=1:p
    r(m)=Rx(1,fs+q+m);
end
A=-(inv(R'*R))*R'*r';
A(p+1)=1;
figure(1)
freqz(1,A,fs,1);
title('AR:p=4')
grid on;

%p2=6时：
p=p2;
R=zeros(p,p);
r=zeros(1,p);
for i=1:p
    for j=1:p
        R(i,j)=R(i,j)+Rx(1,fs+q-p+i+j-1);
    end
end
for m=1:p
    r(m)=Rx(1,fs+q+m);
end
A=-(inv(R'*R))*R'*r';
A(p+1)=1;
figure(2)
freqz(1,A,fs,1);
title('AR:p=6')
grid on;

运行结果：

 http://1851.img.pp.sohu.com.cn/images/blog/2008/12/11/20/18/11ece9ecfb6g213.jpg二、SVD_TLS算法

clear;
clc;

fs=512;
m=0:(fs-1);
xn=sqrt(20)*sin(2*pi*0.2*m)+sqrt(2)*sin(2*pi*0.213*m)+randn(1,fs); %产生含有噪声的序列xn
Rx=xcorr(xn);

pe=100;%%%%%%%%%%%%%%%%%    pe>=p;
qe=150;%%%%%%%%%%%%%%%%    qe>=q,(qe-pe)>=(q-p);
M=300;%%%%%%%%%%%%%%%%%    M>=pe;
Re=zeros(M,pe+1);
for i=1:M
    for j=1:(pe+1)
        Re(i,j)=Rx(1,512+qe+i-j+1);
    end
end

[U,S,V]=svd(Re);%%%%%%奇异值分解
%%求p值
p=0;
SS=zeros(1,pe+1)';
for i=1:(pe+1)
    if S(i,i)/S(1,1)>0.05
        p=p+1;
        if p>0
            SS(p)=S(i,i);
        end
    end
end
SS'
p
%%求Sp部分
Sp=zeros(p+1,p+1);
for j=1:p
    for i=1:(pe+1-p)
       Sp=Sp+SS(p)^2*V(i:i+p,j)*V(i:i+p,j)';
    end
end
%%求Sp逆矩阵
inv_Sp=inv(Sp);
%%求x
x=zeros(1,p+1)';
for i=1:(p+1)
    x(i)=inv_Sp(i,1)/inv_Sp(1,1);
end
A=x;
freqz(1,A,512,1);
grid on;
title('SVDTLS')

仿真结果：

http://1802.img.pp.sohu.com.cn/images/blog/2008/12/11/20/21/11ecea5ce7dg214.jpg尽管不是很好，仅是为了共享一下，以供自己以后参考。