shibbsr.m

将基于最大峰值和ＪＡＤＥ盲分离算法相结合的一种盲分离算法．

function [W]=shibbsR(X,m) 
    [n,T]=size(X);
    verbose=0; %信息显示控制变量  verbose=0不显示提示信息，verbose=！0显示提示信息
      seuil=0.01/sqrt(T);%联合对角上的统计意义上的门限
%--------------确定要恢复的源信号个数-------------------------------------
     if nargin==1
           m=n;
     end  %源信号数据默认为观测信号个数
     if m>n
         fprintf('shibbs->源信号个数不能大于传感器数目!\n')
     return
     end
     if verbose
         fprintf('shibbs->分离的源信号数目m\n',m);
     end
%----------------------去均值----------------------------------------
 if verbose
         fprintf('shibbs->去均值\n');
 end
   meanX=mean(X,2)
   X=X-meanX(:,ones(1,T));
   clear meanX;    %清内存临时变量
%--------------------白化，把信号投影到主分量子空间（提取主分量）-------
 if verbose
         fprinf('shibbs->Whitening the data\n');
 end
%      Rx=X*X';
%      [U,D]=eig(Rx);
%      H=((D)^(-0.5))*U';
%      B=diag(Rx-H*H');
%      B=pinv(B);
%      m=pinv(H'*B*H)
%      Q=m*H'*B;
%      X=Q*X;
     [U,D]=eig((X*X')/T);
     [puiss,k]=sort(diag(D));
     rangW=n-m+1:n;  %indices to the m most significant directions
     scales=sqrt(puiss(rangW));  %scales
     W=diag(1./scales)*U(1:n,k(rangW))';  %白化
       X=W*X;
%        clear W;
%-------------------估计累计量矩阵------------------------------------
     nbcm=m;  %累计量矩阵个数
     CM=zeros(m,m*nbcm);%分配累计量矩阵的存储空间
%--------------------临时变量-----------------------------
     Rk=zeros(m);
     R=eye(m);
     Xk=zeros(m,T);
     xk=zeros(1,T);
     Uns=ones(m,1);
     onemorestep=1;
     while onemorestep
           if verbose
               fprintf('shibbs->估计累计量矩阵\n');
           end
           range=1:m;   % will index the columns of CM where to store the cum.mats
           for k=1:m
               xk=X(k,:);
               Xk=X.*xk(Uns,:);
               Rk=(Xk*Xk')/T-R;
               Rk(k,k)= Rk(k,k)-2;
               CM(:,range)=Rk;
               range=range+m;
           end
%-------------------累计量矩阵的对角化---------------------------------
   %-----------------初始化---------------------------
   V=eye(m);     %旋转阵初化
   nbrs=1;       %扫面时旋转次数变量
  sweep=0;       %扫描次数
  updates=0;
  g=zeros(2,nbcm);
  gg=zeros(2,2);
  G=zeros(2,2);
  c=0;
  s=0;
  ton=0;
  toff=0;
  theta=0;
%--------------------联合对角化处理--------------------------
   if verbose
       fprintf('shibbs->联合对角化优化对比函数\n');
   end
   while nbrs
           nbrs=0;
       sweep=sweep+1;
        if verbose
             fprintf('shibbs->Sweep#%d\n',sweep);
        end
        for p=1:m-1
             for q=p+1:m
                 Ip=p:m:m*nbcm;
                 Iq=q:m:m*nbcm;
 %---------------计算Givens变换角度--------------------
         g=[CM(p,Ip)-CM(q,Iq);CM(p,Iq)-CM(q,Ip)];
         gg=g*g';
         ton=gg(1,1)-gg(2,2);
         toff=gg(1,2)+gg(2,1);
         theta=0.5*atan2(toff,ton+sqrt(ton*ton+toff*toff));
         
        %---------------Givens变换---------------------
       if abs(theta)>seuil
           nbrs=nbrs+1;
           c=cos(theta);  
           s=sin(theta);
           G=[c -s;s c];
           pair=[p;q];
           V(:,pair)= V(:,pair)*G;
           CM(pair,:)=G'*CM(pair,:);
           CM(:,[Ip Iq])=[c*CM(:,Ip)+s*CM(:,Iq)-s*CM(:,Ip)+c*CM(:,Iq)];
      end
     end
  end
  if  verbose
     fprintf('%d步旋转完成\n',nbrs);
  end
  updates=updates+nbrs;
 end
rotsize=norm(V-eye(m),'fro');
if verbose
    fprintf('shibbs->旋转量:%14.6f\n',rotsize);
end
X=V'*X;
W=V'*W;

if rotsize<(m*seuil)
    onemorestep=0;
  end
end

   %---------------行置换求第一分量，这个信号被规格化为单位协方差--------------
 if verbose
     fprintf('shibbs->正存储数据\n',updates);
 end
     A=pinv(W);
    [vars,keys]=sort(sum(A.*A));
    W=W(keys,:);
       W=W(m:-1:1,:);
 %-------------------固定符号，使W的第一行非负-------------------------------
 if verbose
     fprintf('shibbs->固定符号\n',updates);
 end
 b=W(:,1);
 signs=sign(sign(b)+0.1);
 W=diag(signs)*W;