fastica.cpp

本程序将独立分量分析技术和数字水印技术有机地结合在一起

#include <malloc.h>
#include <stdio.h>
#include <math.h>

void rowcentre (float *, int, int, float * );
void mat_transpose (float *, int , int , float *);
void mat_mult (float *, int , int , float *, int , int , float *);
int svd(float *, int , int , float *, float *);
void calc_K(float *, int , int , float *);
void mat_orth (float *, int , float *);
void gramsch (float *, int , int , int );
void row_std (float *, int , int , int );
void Def_logcosh (float *, int , float *, int , int , float *);
void Def_logcosh1 (float *, int , float *, int , int , float , float *);

/*FastICA算法子程序
 *该子程序应用的是A.Hyvarinen和E.Oja在Independent component analysis: algorithms and applications
 *提出的快速定点ICA算法,论文可以在 www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf 获得
 */

void FastICA(float *X, float *W_orig, int nn, int pp, int ee,int maxit,float lim, float *X_Pre,
			 float *xmean, float *B)
/*参数解释
 *X          : 要分离的数据矩阵X
 *W_orig     : 初始分解矩阵W_orig
 *nn         : 数据矩阵的行
 *pp         : 数据矩阵的列
 *ee         : 独立分量的个数
 *maxit      : 最大的循环次数
 *lim        : 收敛值
 *X_Pre      : 中心化后的数据矩阵
 *xmean      : 中心化时的平均值
 *B          : 分离矩阵B=W_final*K

 */
{
	
	int i, j, k, n, p, e,usedNlinearity;
	float tol,tol1,step,stroke,Long; 
	float *W_init,*temp_w1, *temp_w2, *wOld, *wOld2;
	float *W_final,*W1;
	/*Z白化后的数据*/
	float *Z, *K1, *temp1;
	
	
	
	n = nn;
	p = pp;
	e = ee;
	

	/* 预处理矩阵Z的初始化，直接复制原始数据*/
	//X_Pre =(float *) malloc(n * p* sizeof(float));
	for (i = 0; i < n; i++) {
		for (j = 0; j < p; j++) { 
			X_Pre[i * p + j] = X[i * p + j];
		}
	}
	
	/* 行中心化*/
	rowcentre (X_Pre, n, p, xmean);
	printf ("Centering\n");

	
	/*计算白化矩阵K */

    printf ("Whitening\n");
	K1 = (float *)malloc(n * n * sizeof(float));
    Z  = (float *)malloc(n * p * sizeof(float));
	calc_K(X_Pre, n, p, K1); 
	mat_mult(K1, e, n, X_Pre, n, p, Z);
	
	/* 初始分离矩阵W_orig正交化为W_init*/
	temp1 = (float *)malloc(e * e* sizeof(float));	
	W_init = (float *)malloc(e * e* sizeof(float));
	for (i = 0; i < e; i++) {
	    for (j = 0; j < e; j++) {
			temp1[i * e + j] = W_orig[i * e + j];
		}
	}
	mat_orth(temp1, e, W_init);     
	
	/*数据预处理完成*/
	printf("Completed pre-processing\n");
	
	
	/* 主要的ICA代码 
	 * 目标函数：负熵最大化估计，函数选择：tanh
	 * 使用单个估计
	 */
	
    temp_w1 = (float *)malloc(e* sizeof(float));
	temp_w2 = (float *)malloc(e* sizeof(float));
	    
	printf ("Deflation FastICA using logcosh approx. to neg-entropy function\n");
    
	//搜索基向量的过程重复numOfIC次
      i = 0;
      while (i < e)
	  { 
	    //步长定义 
	    step = 1;
        //非线性函数选择
        usedNlinearity = 0;
	    //步长修正
        stroke = 0;
        Long = 0;
        // 显示进程...
        printf("IC %d ", i+1);
        wOld = (float *)malloc(e* sizeof(float));
        wOld2 = (float *)malloc(e* sizeof(float));
        for (j = 0; j < e; j++){
		    wOld[j]=0;
		    wOld2[j]=0;
		}
        //每个基向量的搜索循环...
        k = 1;
        while (k <= maxit + 1)
		{ 
	      //克莱姆-施密特正交化
	      gramsch (W_init, e, e, i + 1); 
	      //行归一化
	      row_std  (W_init, e, e, i + 1);
          // 显示进程...
          printf("."); 
          //测试搜索终止条件。当w和wOld的方向一致，算法收敛。
	      tol = 0.0;
	      for (j = 0; j < e; j++) {
              tol += ((wOld[j]) * (W_init[i * e + j]));
		  }
	      tol = (float)(fabs (fabs (tol) - 1));
          if (tol < lim )
		  { 
          
              //显示进程...
              printf("computed ( %d steps ) \n", k);
              break;//跳出循环
		  }   
          else  
		  {
		      tol1 = 0.0;
	          for (j = 0; j < e; j++) {
                  tol1 += ((wOld2[j]) * (W_init[i * e + j]));
			  }
	          tol1 = (float)(fabs (fabs (tol1) - 1));
		      if ((!stroke) & (tol1<lim))
			  {
			      stroke = step;
	              printf("Stroke!"); 
	              step = 0.5*step;
	              if (usedNlinearity%2 == 0)
                      usedNlinearity = usedNlinearity + 1;
			  }
	          if (stroke)
			  {
			      step = stroke;
	              stroke = 0;
	              if ((step == 1) & (usedNlinearity%2) != 0)
	                  usedNlinearity = usedNlinearity - 1;
			  }
	          if ((!Long) & (i > maxit / 2))
			  {   
			      printf("Taking long (reducing step size) "); 
	              Long = 1;
	              step = 0.5*step;
			      if( (usedNlinearity%2) == 0 )
	                  usedNlinearity = usedNlinearity + 1;
			  }

		  }
	      for (j = 0; j < e; j++) {
			    wOld2[j] = wOld[j];
				    }
		  for (j = 0; j < e; j++) {
			    wOld[j] = W_init[i * e + j];
				    }
      
          switch (usedNlinearity)
          {
	      // 非线性函数选择tanh
          case 0:
	      for (j = 0; j < e; j++) {
			    temp_w1[j] = W_init[i * e + j];
		  }
		  Def_logcosh (temp_w1, e, Z, e, p, temp_w2);
		  for (j = 0; j < e; j++) { 
		        W_init[i * e + j] = temp_w2[j];
		  }
		  break;
          case 1:
		  for (j = 0; j < e; j++) {
			    temp_w1[j] = W_init[i * e + j];
		  }
		  Def_logcosh1 (temp_w1, e, Z, e, p, step, temp_w2);
		  for (j = 0; j < e; j++) { 
		        W_init[i * e + j] = temp_w2[j];
		  }
		  break;
	      
          default: printf("Code for desired nonlinearity not found!");
          }

          // 归一化新的w.
          row_std  (W_init, e, e, i + 1);
          k = k + 1;
		}
        i = i + 1;
	  }
      printf("Done.\n"); 


	  W_final = (float *)malloc(16*16*sizeof(float));	
	  W1 = (float *)malloc(16*16*sizeof(float));
	  
	  for (i = 0; i < e; i++) {
	    for (j = 0; j < e; j++) {
			    W_final[i * e + j] = W_init[i * e + j];
		}
	} 
	    
    //计算估计的分离矩阵B
	
	mat_mult (W_final, e, e, K1, e, n, W1);
    for (i = 0; i < e; i++) {
	    for (j = 0; j < e; j++) {
			    B[i * e + j] = W1[i * e + j];
		}
	} 

 }