jader.m

一种比较常用的盲分离算法JADE(扩展对角化盲分离算法)的源代码.

function [W]=jadeR(X,m)
    verbose=0; %信息显示控制变量
 
    %确定源信号个数
     [n,T]=size(X);
     if nargin==1
           m=n;end;  %源信号数据默认为观测信号个数
     if m>n
         fprintf('jade->源信号个数不能大于传感器数目!\n')
     return
     end
     if verbose
         fprintf('jade->Looking for %d sources\n',m);
     end
     %-------------预处理-----------------------
     if verbose
         fprintf('jade->去均值\n');
     end
     X=X-mean(X')'*ones(1,T);
     %-----------------白化确定抽取信号个数--------------------
     if verbose
         fprinf('jade->Whitening the data\n');
     end
%      Rx=X*X';
%      [U,D]=eig(Rx/T);
%      H=U*(D.^(0.5));
%      B=(Rx-H*H');
%      B=pinv(B);
%        j=inv(H'*B*H)
%       Q=j*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))';  %白化
       iW=U(1:n,k(rangW))*diag(scales);%its pseudo-inverse
       X=W*X;
       %----------------累计量矩阵估计-------------------------
       if verbose
           fprintf('jade->估计累计量矩阵\n');
       end
       dimsymm=(m*(m+1))/2;   %实对称矩阵的维数
       nbcm=dimsymm;    %累计量矩阵的个数
       CM=zeros(m,m*nbcm);  %存储累计量矩阵
      R=eye(m);
      Qij=zeros(m);   %累计量矩阵的临时变量
      Xim=zeros(1,m);   %临时
      Xjm=zeros(1,m);    %临时
      scale=ones(m,1)/T;
      Range=1:m;  %累计量矩阵存储的指数
      for im=1:m
        Xim=X(im,:);
        Qij=((scale*(Xim.*Xim)).*X)*X'-R-2*R(:,im)*R(:,im)';
        CM(:,Range)=Qij;
        Range=Range+m;
         for jm=1:im-1
             Xjm=X(jm,:);
             Qij=((scale*(Xim.*Xim)).*X)*X'-R(:,im)*R(:,jm)'-R(:,jm)*R(:,im)';
             CM(:,Range)=sqrt(2)*Qij;
              Range=Range+m;
         end
      end
      %------------------累计量矩阵的联合对角化-------------------
         %---------------初始化------------------------
       if 1    %对角化信号累计量初始化，节约计算时间
             if verbose
                 fpritf('jade->对角初始化\n');
             end
             [V,D]=eig(CM(:,1:m));
             for u=1:m:m*nbcm
                   CM(:,u:u+m-1)=CM(:,u:u+m-1)*V;
             end
             CM=V'*CM;
       else  V=eye(m);
       end
       
       seuil=1/sqrt(T)/100;   %统计门限
       encore=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('jade->利用联合对角化优化对比函数\n');
    end
    while encore
      encore=0;
         if verbose
             fprintf('jade->Sweep#%d\n',sweep);
         end
         sweep=sweep+1;
         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
           encore=1;
           updates=updates+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
 end
 if  verbose
     fprintf('jade->Total of %d Givens totations\n',updates);
 end
 W=V'*W;   %分离矩阵
 %---------------行置换求第一分量，这个信号被规格化为单位协方差--------------
 if verbose
     fprintf('jade->分量排序\n',updates);
 end
 A=iW*V;
 [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;
%  W=V*W;