schura.m

盲信号处理牛人cardoso

function [V,D]=schurb(A,jthresh)
% Joint approximate Schur transformation
% 
% Joint approximate of n (complex) matrices of size m*m stored in the
% m*mn matrix A by minimization of a joint diagonality criterion
%
% Usage:  [V,D]=schur-b(A,jthresh)
%
% Input :
% * the m*nm matrix A is the concatenation of n matrices with size m
%   by m. We denote A = [ A1 A2 .... An ]
% * threshold is an optional small number (typically = 1.0e-8 see the M-file).
%
% Output :
% * V is an m*m unitary matrix.
% * D = V'*A1*V , ... , V'*An*V has the same size as A and is a
%   collection of upper triangular matrices if A1, ..., An are exactly jointly
%   eigenvalue decomposition.
%
% The algorithm finds a unitary matrix V such that the matrices
% V'*A1*V , ... , V'*An*V are as upper triangular as possible, providing a
% kind of `average eigen-structure' shared by the matrices A1 ,...,An.
% If the matrices A1,...,An do have an exact common eigen-structure ie
% a common orthonormal set eigenvectors, then the algorithm finds it.
% The eigenvalues are the diagonal elements of D1, ...,Dn .
% 
% The algorithm implements a properly extended Jacobi algorithm.  The
% algorithm stops when all the Givens rotations in a sweep have sines
% smaller than 'threshold'.
%
% In many applications, the notion of approximate joint
% diagonalization is ad hoc and very small values of threshold do not
% make sense because the diagonality criterion itself is ad hoc.
% Hence, it is often not necessary in applications to push the
% accuracy of the rotation matrix V to the machine precision.
%
% PS: If a numrical analyst knows `the right way' to determine jthresh
%     in terms of 1) machine precision and 2) size of the problem,
%     I will be glad to hear about it.
% 
%
% This version of the code is for complex matrices, but it also works
% with real matrices.  However, simpler implementations are possible
% in the real case.
%
% See more info, references and version history at the bottom of this
% m-file

%
%----------------------------------------------------------------
% Version 1.3
%
% Copyright     : Da-Zheng Feng. 
% Author        : Da-Zheng Feng. [e-mail]
% Comments, bug reports, etc are welcome.
%----------------------------------------------------------------


[M,NM] = size(A);
N=fix(NM/M);

V=eye(M,M);
g=zeros(1,3);
M
N

k=0
while 1
   for m=1:M-1
      for n=m+1:M         
         G=zeros(3,3);
         for p=1:N               
            pp=(p-1)*M;
            mm=pp+m;
            nn=pp+n;
            g(1)=A(m,mm)-A(n,nn);
            bb=A(m,nn);
            cc=A(n,mm);
            g(2)=bb+cc;
            g(3)=j*(cc-bb);
            G=G+g'*g;
         end
         G=real(G);
         [ej,dd]=eig(G);
         [la,KK]=sort(diag(dd));
         v=ej(:,KK(3));   
                  
         if v(1)<0, v=-v; end
         c=sqrt(0.5+v(1)/2);
         s=0.5*(v(2)-j*v(3))/c;
         pair=[m;n];
         G=[c -conj(s);s c];
         %-----------------------------
         V(:,pair)=V(:,pair)*G;
         A(pair,:)=G'*A(pair,:);
         for p=1:N
            pp=(p-1)*M;
            pair1=[pp+m pp+n];
            A(:,pair1)=A(:,pair1)*G;
         end              
         %-----------------------------------------------   
      end
   end
   if abs(s)<jthresh, break, end
end
D = A ;

return