JadeMatlab代码

jade的matlab代码，可以用于盲信号分离，具有收敛速度快的特性，比fastica还快

function [ V  D ] =  joint_diag（Ajthresh）
% Joint approximate diagonalization
% 联合近似对角化
%
% Joint approximate of n （complex） matrices of size m*m stored in the
% m*mn matrix A by minimization of a joint diagonality criterion
% 对储存在A中的n个mxm矩阵做联合对角化
%
% Usage:  [ V  D ] =  joint_diag（Ajthresh）
%
% 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 diagonal matrices if A1 ... An are exactly jointly
%   unitarily diagonalizable.
%
% The algorithm finds a unitary matrix V such that the matrices
% V‘*A1*V  ...  V‘*An*V are as diagonal 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 eigenvectors THEN are the column vectors of V and D1 ...Dn are
% diagonal matrices.
% 联合对角化的作用。
%
% 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.2
%
% Copyright 	: Jean-Francois Cardoso.
% Author 	: Jean-Francois Cardoso. cardoso@sig.enst.fr
% Comments bug reports etc are welcome.
%----------------------------------------------------------------

[mnm] = size（A）;
B = [ 1 0 0 ; 0 1 1 ; 0 -1i 1i ] ;
Bt = B‘ ;
V = eye（m）;
encore	= 1;
while encore
    encore=0;
    for p=1:m-1
        Ip = p:m:nm;
        for q=p+1:m
            Iq = q:m:nm;
            
            % Computing the Givens angles 计算旋转角度
            g = [A（pIp）-A（qIq）; A（pIq）; A（qIp）];
            [vcpD] = eig（real（B*（g*g‘）*Bt））;
            [~ K] = sort（diag（D））;
            angles  = vcp（:K（3））;
            if angles（1）<0
                angles= -angles;
            end
            c = sqrt（0.5+angles（1）/2）;
            s = 0.5*（angles（2）-1j*

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4903  2013-06-06 11:03  joint_diag.m
     文件         165  2013-06-06 14:49  joint_diag_Note.txt
     文件        3813  2013-06-06 11:05  joint_diag_r.m
     文件        1021  2013-06-06 11:00  test_joint_diag.m
     文件         970  2013-06-06 11:08  test_joint_diag_r.m