jdb2.m

盲信号处理牛人cardoso

% Calling  the joint approximate diagonalization function.

M=3   % dimension
N=3   % Number of matrices
snr=0.025
seuil   = 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);
for imat=1:N
  cols          = 1+(imat-1)*M:imat*M;
  A(:,cols)     = U*diag(randn(M,1)+i*randn(M,1))*V;
  
end

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


%% Do it
[ V , DD ] = joint_diag(A,seuil);

%% should be permutation matrix 
abs(V'*U)
DD=abs(DD);
DD

end