calc_k.cpp

来自「本程序将独立分量分析技术和数字水印技术有机地结合在一起」· C++ 代码 · 共 57 行

#include <malloc.h>
#include <math.h>
void mat_transpose (float *, int , int , float *);
void mat_mult (float *, int , int , float *, int , int , float *);
int svd(float *, int , int , float *, float *);
/*calc_K-求白化矩阵子程序
 *白化的过程为:1.求去均值后的数据矩阵的协方差阵 c=x*x'；2.对c做SVD分解 3. k=inv(sqrt(D))*u'
 */

void calc_K(float *x, int n, int p, float *K)
{
    int i,j;
    float *xxt, *xt, *u, *d, *v, *temp1, *temp2;
    
    xt = (float *)malloc(n * p* sizeof(float));
    xxt = (float *)malloc(n * n* sizeof(float));
   
	/* 求x的转置矩阵 xt*/
    mat_transpose (x, n, p, xt); 
    /* 计算样本的协方差阵xxt */
    mat_mult (x, n, p, xt, p, n, xxt); 
	
	for (i = 0; i < n; i++) {
	    for (j = 0; j < n; j++) {
		    xxt[n * i + j] = xxt[n * i + j] / p;
	    }
    }	
    

		    
    /* 对xxt做svd分解，得到矩阵u，d，v*/ 
    u = (float *)malloc(n * n* sizeof(float));
    d = (float *)malloc(n* sizeof(float));
    v = (float *)malloc(n * n* sizeof(float));
    svd (xxt, n, n, d, v); 


    /* 计算白化矩阵K*/
    temp1 = (float *)malloc(n * n* sizeof(float));
    temp2 = (float *)malloc(n * n* sizeof(float));
    
	
    for (i = 0; i < n; i++) {
		for(j = 0; j < n; j++){
			if(i==j)
				temp1[n * i + j] = (float)(1 / sqrt (d[i]));
			else
				temp1[n * i + j] = 0.0;
		}
	    
	
    }

    mat_transpose (xxt, n, n, temp2);
    mat_mult (temp1, n, n, temp2, n, n, K);
}