ogwe.m

基于累计量的独立分量分析的改进算法

function w=ogwe(x,m,wkur,wxi)
[n,T]=size(x);
verbose=0; %信息显示控制变量  
%-----------判断输入参数个数---------------------
switch(nargin)
case 1
     wxi=3/7;
     wkur=0;
     m=n;
case 2
     if m-floor(m)>0 | m<=0
         fprintf('m should be a positive integer!\n')
         return
     end
     if m>n
         fprintf('no more sources than sensors here!\n')
         return
     end
     if verbose
         fprintf('looking for %d sources\n',m)
     end
     wxi=3/7;
     wkur=0;
case 3
     wxi=m;
     wkur=wxi;
     m=n;
case 4
      if m-floor(m)>0 | m<=0
         fprintf('m should be a positive integer!\n')
         return
     end
     if m>n
         fprintf('no more sources than sensors here!\n')
         return
     end
     if verbose
         fprintf('looking for %d sources\n',m)
     end
end
%-----------------去均值-----------------------------------
 if verbose
         fprintf('zero mean\n');
 end
 x=x-mean(x')'*ones(1,T); %去均值
 %---------------白化-------------------------------------------
 if verbose
         fprintf('白化\n');
 end
 [F,D]=eig((x*x')/T);
 [vals,idxs]=sort(diag(D));
 mmvar=n-m+1:n;
 F=F(:,idxs(mmvar));
 D=D(idxs(mmvar),idxs(mmvar));
 v=sqrt(diag(1./diag(D)))*F;
 x=v*x;
 %--------------最小化对比函数到单位矩阵R-----------------------
  if verbose
         fprintf('优化对比函数：计算单位矩阵\n');
  end
 if m==2
     k=1;
 else
     k=1+round(sqrt(m));
 end
 seuil=pi/360;  %角度门限
 burden=1;
 sweep=0;
 givensrotations=0;
 R=eye(m);  %变换矩阵的初值
 %------------选择不初始化SICA或初始化SICA---------------------
 a=123.75;
 b=-237.5;
 y=a*n+b;
 if n<15 & T>y
     initalization=1; %'ISICA'; %disp('ISICA');
 else
     initalization=0; %'ISICA'; %disp('ISICA');
 end
 if  initalization
     nm=m*(m+1)/2;
     x4=zeros(nm,nm);
     ix1=1;
     ix2=m;
%--------------旋转向量和累计量矩阵初始化-----------------------------------
 if verbose
         fprintf('初始化旋转向量和累计量矩阵\n');
 end
 Rpp=zeros(1,nm);
 Rqq=zeros(1,nm); 
 Rpq=zeros(1,nm);
 Xpp=zeros(1,T);
 p1=[];
 q1=[];
 ix=1;
 for p_1=1:m;
     sccero=0.5;
     scuno=1;
     for q_1=p_1:m
         sc1(ix)=scuno;
         sc0(ix)=sccero;
         p1=[p1 p_1];
         q1=[q1 q_1];
         scuno=2;
         sccero=1;
         ix=ix+1;
     end
 end
 %-----------------块数-----------------------------------------------
 Nb=m;
 dim=sum(m:-1:1);
 ixq=[];vq=1;
 ixk1=[];
 ixr=[];
 for l=Nb:-1:1
     ixq=[ixq vq*ones(1,l)];
     ixk1=[ixk1 [vq:Nb]];
     vq=vq+1;
     ixr=[ixr sum(Nb:-1:l+1)+1];
 end
 ixr=[ixr ixr(m)];
 Xqk1=zeros(ixr(m),T);
 Xqk1=x(ixq,:).*x(ixk1,:);
 for p=1:Nb
     k2=p:m;
     idc=ixr(p)+k2-p;
     idr=ixr(k2(1)):dim;
     X4(idr,idc)=Xqk1(idr,:)*(Xqk1(idc,:))'/T;
     X4(idc,idr)=X4(idr,idc)';
     idc=ixr(p);
     for k2=p:m
         for k3=idc:ixr(k2)-1
             ixcoj=[sort([ixk1(k3) k2])];
             col=ixr(p)-p+ixq(k3);
             row=ixr(ixcoj(1))-ixcoj(1)+ixcoj(2);
             X4(k3,idc)=X4(row,col);
             X4(idc,k3)=X4(k3,idc);
         end
         idc=idc+1;
     end
 end
 %--------------------------对比函数最小化-----------------------------
while burden & sweep<k
    burden=0;
    if verbose
         fprintf('sweep #%d\n',sweep);
    end
 sweep=sweep+1;
 for p=1:m-1
     for q=p+1:m
%------------------------计算旋转向量需要的累计量矩阵--------------------
       Rpp=sc1.*R(p,p1).*R(p,q1);
       Rqq=sc1.*R(q,p1).*R(q,q1);
       Rpq=sc0.*(R(p,p1).*R(q,q1)+R(p,q1).*R(q,q1));
%-------------------计算累计量-------------------------------------------
        RppX4= Rpp*X4;
        X4Rqq=X4*Rqq';
        Mppqq= RppX4*Rqq';
        Mpppp= RppX4*Rpp';
        Mqqqq= Rqq*X4Rqq;
        Mpqqq= Rpq*X4Rqq;
        Mpppq= RppX4*Rpq';
        r4=Mpppp+Mqqqq+2*Mppqq;
        sxpxqr2=2*(Mpqqq+Mpppq);
        r4c4phi=r4-8*Mppqq;
        r4s4phi=8*Mpppq-2*sxpxqr2;
        r4c2phi=Mpppp-Mqqqq;
        r4s2phi=sxpxqr2;
%-------------------givens角度计算---------------------------------------

        r4c2phi=r4c2phi^2;
        r4s2phi=r4s2phi^2;
        r4_c=r4-8;
        xi4=r4c4phi+j*r4s4phi;
        xi2_2=(r4c2phi+j*r4s2phi)^2;
        theta=angle(wkur*xi4+(1-abs(wkur))*(wxi*r4_c*xi4+(1-wxi)*xi2_2))/4;
 %--------------------------givens变换---------------------------------------
        givensrotations=givensrotations+1;
        if abs(theta)>seuil
            burden=1;
            c=cos(theta);
            s=sin(theta);
            G=[c s;-s c]
            pair=[p;q]
            R(pair,:)=G*R(pair,:);
        end
     end
 end
end
else
while burden & sweep<k
    burden=0;
    if verbose
         fprintf('sweep #%d\n',sweep);
    end
 sweep=sweep+1;
 for p=1:m-1
     for q=p+1:m
         Mpp=x(p,:).*x(p,:);
         Mqq=x(q,:).*x(q,:);
         Mpq=x(p,:).*x(q,:);
         Mpppp=Mpp*Mpp'/T;
         Mqqqq=Mqq*Mqq'/T;
         Mpqqq=Mpq*Mqq'/T;
         Mpppq=Mpp*Mpq'/T;
         Mppqq=Mpp*Mqq'/T;
        r4=Mpppp+Mqqqq+2*Mppqq;
        sxpxqr2=2*(Mpqqq+Mpppq);
        r4c4phi=r4-8*Mppqq;
        r4s4phi=8*Mpppq-2*sxpxqr2;
        r4c2phi=Mpppp-Mqqqq;
        r4s2phi=sxpxqr2;
%-------------------givens角度计算---------------------------------------

        r4c2phi=r4c2phi^2;
        r4s2phi=r4s2phi^2;
        r4_c=r4-8;
        xi4=r4c4phi+j*r4s4phi;
        xi2_2=(r4c2phi+j*r4s2phi)^2;
        theta=angle(wkur*xi4+(1-abs(wkur))*(wxi*r4_c*xi4+(1-wxi)*xi2_2))/4;%GWE
 %--------------------------givens变换---------------------------------------
        givensrotations=givensrotations+1;
        if abs(theta)>seuil
            burden=1;
            c=cos(theta);
            s=sin(theta);
            G=[c s;-s c]
            pair=[p;q]
            R(pair,:)=G*R(pair,:);
            x([p q],:)=G*x([p q],:);
        end
     end
  end
 end
end
if verbose
         fprintf('%d givens angle computed\n',givensrotations);
end
 %------------------最终的结果-----------------------------------------
 w=R*v;
 if verbose
         fprintf('sorting components\n');
 end
 P=inv(D)*R';
 [vals,idxs]=sort(sum(P.*P));
 w=w(idxs,:);