jda.m

盲信号处理牛人cardoso

% Calling  the joint approximate diagonalization function.

M=3;   % dimension
N=3;   % Number of matrices
snr=0.25;
jthresh   = 1.0e-12; % precision on joint diag 
K=1;

for k=1:K
   k;

% drawing a `random' unitary matrix
U= randn(M)+i*randn(M) ; 
[U,to_waste]=eig(U+U'); 
%U=randn(M,M)+i*randn(M,M);
V=U';
%V=inv(U);
% Drawing a random set of commuting matrices
A=zeros(M,M*N);
by=zeros(M,M);
for imat=1:N
   imat1=(imat-1)*M;
   for m=1:M
      by(m,m)=randn+j*randn;
   end
   
   A(1:M,imat1+1:imat1+M) = U*by*V;  
end

% Perturbation of the joint structure ?
 A = A + snr*randn(M,M*N);
AA=A;

%% Do it
tic,[V,D]=joint_diag(AA,jthresh),toc
by1=abs(inv(V)*U);

tic,[V,DD]=schura(A,jthresh),toc

%% should be permutation matrix 
by2=abs(inv(V)*U);
by1
by2
DD=abs(DD);
DD;
end