arma_svdtls.m

基于累计量的奇异值－总体最小二乘法求AR参数 用奇异值－总体最小二乘法求AR参数 一般最小二乘法求AR参数 根据AR参数和自相关函数以及AR阶数用Cadzow谱估计子求出频谱密度

function [a,Rx,fTLS,p]=ARMA_SVDTLS(M,x,pe,varargin)

%pe=60;
[Q,N]=size(x);
%M=40;
%sampling=100;
%fls=zeros(sampling,2);
%p_average=0;
%for k=1:sampling
%     w=randn(size(n));
%     x = sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w;
    Rx=xcorr(x,varargin{:});
  % Rx=xcorr(x,'unbiased');
    R=zeros(M,pe);
    for i=1:M
        R(i,:)=Rx(pe+i-1+N:-1:i+N);
    end
    b=-Rx(pe+1+N:pe+M+N)';
    B=[-b,R];
    [U,S,V]=svd(B);
    S=diag(S);
    threshold=0.9999;
    for p=1:min(M,pe+1)
        v=sqrt(sum(S(1:p).^2)/sum(S.^2));
        if v>=threshold
            break;
        end
    end
    Sp=zeros(p+1,p+1);
    for j=1:p
        for i=1:pe+1-p
            Sp=Sp+(S(j)^2)*V(i:i+p,j)*(conj(V(i:i+p,j))');
        end
    end

    invSp=inv(Sp);
    a=invSp(2:p+1,1)/invSp(1,1);
    a=[1;a];   %add a0=1
    z=roots(a');
 
    f=atan(imag(z)./real(z))/2/pi;   
    f=sort(f,'descend');

    fTLS=[f(2);f(1)];   %the frequency obtained by TLS method
    %fTLS(k,:)
    %p_average=p_average+p;
    %p
%end
%p=round(p_average/sampling); %为了避免干扰，采用统计平均的方式来确定p值