arma_hosa.m

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

function [a,Rx,fHOSA,p]=ARMA_HOSA(x,M1,M2,N1,N2,varargin)


[Q,N]=size(x);
%N=1000;
%n=1:N;
%sampling=100;
ftemp=zeros(sampling,2);
% M1=4;
% M2=10;
% N1=-5;
% N2=5;
maxlag=M1+2*M2-1;
overlap=0;
%M=40;
%sampling=100;
%fls=zeros(sampling,2);
%p_average=0;
%for k=1:sampling
%k
% w=randn(1,N+3);
% w3=w(1:N);
% w2=w(2:N+1);
% w1=w(3:N+2);
% w0=w(4:N+3);
% x = sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w0+0.9*w1+0.385*w2+0.771*w3;
Rx=xcorr(x,varargin{:});
%   Rx=xcorr(x,'unbiased');
  temp=zeros(N2-N1+1,2*(M1+2*M2-1)+1);
    Ce=zeros((M2+1)*(N2-N1+1),M2);
    for i=1:(N2-N1+1)
      temp(i,:)=cumest(x,3,maxlag,sampling,overlap,varargin{:},i-1,0);
%        temp(i,:)=cumest(x,3,maxlag,sampling,overlap,'unbiased',i-1,0);
    end
    j=1;
    for i=(M1+M2-1):(M1+2*M2-1)
        Ce(j:j+N2-N1,:)=temp(:,((M1+2*M2-1)+1+i):-1:((M1+2*M2-1)+1+i-M2+1));
        j=j+N2-N1+1;
    end
    [U,S,V]=svd(Ce);
    S=diag(S);
%    S
    threshold=0.90;
    for p=1:min((M2+1)*(N2-N1+1),M2)
        v=sqrt(sum(S(1:p).^2)/sum(S.^2));
        if v>=threshold
            break;
        end
    end
 %   p=5;
    Sp=zeros(p+1,p+1);
    for j=1:p
        for i=1:M2-p
            Sp=Sp+(S(j)^2)*V(i:i+p,j)*(conj(V(i:i+p,j))');
        end
    end
  %  p
    invSp=inv(Sp);
    a=invSp(2:p+1,1)/invSp(1,1);
    a=[1;a];   %add a0=1
 %   a
    z=roots(a');
 
    f=atan(imag(z)./real(z))/2/pi;   
    f=sort(f,'descend');
 %   f

    fHOSA=[f(2),f(1)];   %the frequency obtained by TLS method
    
%     ftemp(k,:)=fHOSA;
%     Cadzow(a,Rx,p,N);
%     hold on;
    %fTLS(k,:)
    %p_average=p_average+p;
    %p
%end
% fu=mean(ftemp)
% fvar=var(ftemp)
% hold off;
%p=round(p_average/sampling); %为了避免干扰，采用统计平均的方式来确定p值