www.pudn.com > JnS-1.2.rar > JnS.c
/* ================================================================== */
/*
Implements the Jade and the Shibbs algorithms
Copyright: JF Cardoso. cardoso@tsi.enst.fr
This is essentially my first C program. Educated comments are more
than welcome.
version 1.2 Jun. 05, 2002.
Version 1.1 Jan. 20, 1999.
Changes wrt 1.1
o Minor fix for new versions of mex (see History)
Changes wrt 1.0
o Switched to Matlab-wise vectorization of matrices
o Merged a few subroutines into their callers
o Implemented more C tricks to make the code more unscrutable
o Changed the moment estimating routines to prepare for a
read-from-file operation (the sensor loops are nested inside
the sample loops)
o Limited facility to control verbosity levels. Messages directed
to sterr.
To do:
o Address the convergence problem of Shibbs on (e.g.) Gaussian data.
o Control of convergence may/should be based on the variation of the objective
rather than one the size of the rotations (see above item).
o Smarter use of floating types: short for the data, long when
during moment estimation (issue of error accumulation).
o An `out of memory' should return an error code rather than exiting.
*/
/* ================================================================== */
#include
#include
#include
#include "Matutil.h"
#define VERBOSITY 0
#define RELATIVE_JD_THRESHOLD 1.0e-4
/* A `null' Jacobi rotation for joint diagonalization is smaller than
RELATIVE_JD_THRESHOLD/sqrt(T) where T is the number of samples */
#define RELATIVE_W_THRESHOLD 1.0e-12
/* A null Jacobi rotation for the whitening is smaller than
RELATIVE_W_THRESHOLD/sqrt(T) where T is the number of samples */
void OutOfMemory()
{
printf("Out of memory, sorry...\n");
exit(EXIT_FAILURE) ;
}
#define SPACE_PER_LEVEL 3
void Message0(int level, char *mess) {
int count ;
if (level < VERBOSITY) {
for (count=0; count A(mxn) x R where R=[ c s ; -s c ] rotates the (p,q) columns of R */
void RightRotSimple(double *A, int m, int n, int p, int q, double c, double s )
{
double nx, ny ;
int ix = p*m ;
int iy = q*m ;
int i ;
for (i=0; i Ak(mxn) x R where R rotates the (p,q) columns R =[ c s ; -s c ]
and Ak is the k-th M*N matrix in the stack */
void RightRotStack(double *A, int M, int N, int K, int p, int q, double c, double s )
{
int k, ix, iy, cpt, kMN ;
int pM = p*M ;
int qM = q*M ;
double nx, ny ;
for (k=0, kMN=0; k R * A(mxn) where R=[ c -s ; s c ] rotates the (p,q) rows of R
*/
void LeftRotSimple(double *A, int m, int n, int p, int q, double c, double s )
{
int ix = p ;
int iy = q ;
double nx, ny ;
int j ;
for (j=0; j R * Ak(mxn) where R rotates the (p,q) rows R =[ c -s ; s c ]
and Ak is the k-th matrix in the stack
*/
void LeftRotStack(double *A, int M, int N, int K, int p, int q, double c, double s )
{
int k, ix, iy, cpt ;
int MN = M*N ;
int kMN ;
double nx, ny ;
for (k=0, kMN=0; kqq)
return 0.5 * atan2(-pq-qp, pp-qq) ;
else
return 0.5 * atan2(pq+qp, qq-pp) ;
}
/* Givens angle for the pair (p,q) of a stack of K M*M matrices */
double GivensStack(double *A, int M, int K, int p, int q)
{
int k ;
double diff_on, sum_off, ton, toff ;
double *cm ; /* A cumulant matrix */
double G11 = 0.0 ;
double G12 = 0.0 ;
double G22 = 0.0 ;
int M2 = M*M ;
int pp = p+p*M ;
int pq = p+q*M ;
int qp = q+p*M ;
int qq = q+q*M ;
for (k=0, cm=A; k0) {
encore = 0 ;
for (p=0; p threshold ) {
c = cos(theta);
s = sin(theta);
LeftRotSimple (A, m, m, p, q, c, s) ;
RightRotSimple(A, m, m, p, q, c, s) ;
LeftRotSimple (R, m, m, p, q, c, s) ;
encore = 1 ;
rots++ ;
}
}
}
return rots ;
}
/* Joint diagonalization of a stack K of symmetric M*M matrices by Jacobi */
/* returns the number of plane rotations executed */
int JointDiago (double *A, double *R, int M, int K, double threshold)
{
int rots = 0 ;
int more = 1 ;
int p, q ;
double theta, c, s ;
Identity(R, M) ;
while ( more > 0 ) /* One sweep through a stack of K symmetric M*M matrices. */
{
more = 0 ;
for (p=0; pthreshold) {
c = cos(theta);
s = sin(theta);
LeftRotStack (A, M, M, K, p, q, c, s) ;
RightRotStack(A, M, M, K, p, q, c, s) ;
LeftRotSimple(R, M, M, p, q, c, s) ;
rots++ ;
more = 1 ; /* any pair rotation triggers a new sweep */
}
}
}
return rots ;
}
/* W = sqrt(inv(cov(X))) */
void ComputeWhitener (double *W, double *X, int n, int T)
{
double threshold_W = RELATIVE_W_THRESHOLD / sqrt((double) T) ;
double *Cov = (double *) calloc(n*n, sizeof(double)) ;
double rescale ;
int i,j ;
if (Cov == NULL) OutOfMemory() ;
EstCovMat (Cov, X, n, T) ;
Diago (Cov, W, n, threshold_W) ;
for (i=0; i0)
{
Message0(2, "Computing cumulant matrices...\n") ;
EstCumMats (CumMats, X, nbc, nbs) ;
Message0(2, "Joint diagonalization...\n") ;
rots = JointDiago (CumMats, Transf, nbc, nbc, threshold_JD) ;
MessageI(3, "Total number of plane rotations: %6i.\n", rots) ;
MessageF(3, "Size of the total rotation: %10.7e\n", NonIdentity(Transf,nbc) );
Message0(2, "Updating...\n") ;
Transform (X, Transf, nbc, nbs ) ;
Transform (B, Transf, nbc, nbc ) ;
}
free(Transf) ;
free(CumMats) ;
}
/* _________________________________________________________________ */
void Jade (
double *B, /* Output. Separating matrix. nbc*nbc */
double *X, /* Input. Data set nbc x nbs */
int nbc, /* Input. Number of sensors */
int nbs /* Input. Number of samples */
)
{
double threshold_JD = RELATIVE_JD_THRESHOLD / sqrt((double)nbs) ;
int rots = 1 ;
double *Transf = (double *) calloc(nbc*nbc, sizeof(double)) ;
double *CumTens = (double *) calloc(nbc*nbc*nbc*nbc, sizeof(double)) ;
if ( Transf == NULL || CumTens == NULL ) OutOfMemory() ;
/* Init */
Message0(2, "Init...\n") ;
Identity(B, nbc);
MeanRemoval(X, nbc, nbs) ;
Message0(2, "Whitening...\n") ;
ComputeWhitener(Transf, X, nbc, nbs) ;
Transform (X, Transf, nbc, nbs) ;
Transform (B, Transf, nbc, nbc) ;
Message0(2, "Estimating the cumulant tensor...\n") ;
EstCumTens (CumTens, X, nbc, nbs) ;
Message0(2, "Joint diagonalization...\n") ;
rots = JointDiago (CumTens, Transf, nbc, nbc*nbc, threshold_JD) ;
MessageI(3, "Total number of plane rotations: %6i.\n", rots) ;
MessageF(3, "Size of the total rotation: %10.7e\n", NonIdentity(Transf,nbc) ) ;
Message0(2, "Updating...\n") ;
Transform (X, Transf, nbc, nbs ) ;
Transform (B, Transf, nbc, nbc ) ;
free(Transf) ;
free(CumTens) ;
}