tls.m

本源码给出了总体最小二乘方法的matlab实现

clear
clc
legth=128;
n=1:legth;
%%%%%%计算x（n）-e（n）
x=sqrt(20)*sin(2*pi*0.2.*n)+sqrt(2)*sin(2*pi*0.215.*n);
%%%%计算e（n）这里用w表示
w=randn([1,legth]);
x=x+w;
%%%%计算自相关矩阵Rx（k）
M=40;p=7;
r=zeros(M+1,2*p+1);
for i=1:M+1;
    for j=1:2*p+1;
        %%% r(i,j)=r(i-j)=r(k)
        k=i-j;
        for l=1:legth-k;
            if l+k<1 || l>legth
               continue;
            end
            r(i,j)=r(i,j)+x(l)*x(l+k);
        end
    end
end
%%%%%%
A=r(:,2:end);
b=-r(:,1);
% %%%%%一般最小二乘解
% a=inv(conj(A)'*A)*conj(A)'*b;
% xishu=[a;1];
% xishu=fliplr(xishu);
% f=roots(xishu);
% f1=atan(imag(f)./real(f))/2/pi
%%%%%%svd总体最小二乘
B=r;
[m,n]=size(B);
[s,d,v]=svd(B);
p=5;
%%%%计算sp矩阵
sp=zeros(p+1,p+1);
for j=1:p;
    for i=1:n-p;
       
        sp=sp+d(i,i)^2*v(i:i+p,j)*conj(v(i:i+p,j))';
    end
end
sp1=inv(sp);

xtls=sp1(2:end,1)./sp1(1,1);
a=[xtls;1];
a=fliplr(a);
c=roots(a);
f2=atan(imag(c)./real(c))/2/pi