pca2.cpp

实现pca功能的完整matlab程序

			ppp( a, e, s, v, m, n );
            free( s ); 
			free( e ); 
			free( w ); 
			return( -1 );
        }
        kk = mm-1;
		while ( (kk != 0) && (fabs(e[kk-1]) != 0.0) )
		{ 
			d = fabs(s[kk-1])+fabs(s[kk]);
			dd = fabs(e[kk-1]);
			if ( dd > eps*d ) 
				kk = kk-1;
			else 
				e[kk-1] = 0.0;
		}
		if ( kk == mm-1 )
		{ 
			kk = kk+1;
			if ( s[kk-1] < 0.0 )
			{ 
				s[kk-1] = -s[kk-1];
				for ( i = 1; i <= n; i++ )
				{ 
					ix = (i-1)*n+kk-1; 
					v[ix] = -v[ix];
				}
			}
			while ( (kk!=m1) && (s[kk-1] < s[kk]) )
			{ 
				d = s[kk-1]; 
				s[kk-1] = s[kk]; 
				s[kk] = d;
				if ( kk < n )
					for (i=1; i<=n; i++)
					{ 
						ix = (i-1)*n+kk-1; iy = (i-1)*n+kk;
						d = v[ix]; v[ix] = v[iy]; v[iy] = d;
					}
					if ( kk < m )
						for ( i = 1; i <= m; i++ )
						{ 
							ix = (i-1)*m+kk-1; iy = (i-1)*m+kk;
							d = u[ix]; u[ix] = u[iy]; u[iy] = d;
						}
					kk = kk+1;
			}
			it = 60;
			mm = mm-1;
		}
		else
		{ 
			ks = mm;
			while ( (ks>kk) && (fabs(s[ks-1]) != 0.0 ) )
			{ 
				d = 0.0;
				if ( ks != mm ) d = d+fabs(e[ks-1]);
				if ( ks != kk+1 ) d = d+fabs(e[ks-2]);
				dd = fabs(s[ks-1]);
				if ( dd > eps*d ) ks = ks-1;
				else s[ks-1] = 0.0;
			}
			if ( ks == kk )
			{ 
				kk = kk+1;
				d = fabs(s[mm-1]);
				t = fabs(s[mm-2]);
				if ( t > d ) d = t;
				t = fabs( e[mm-2] );
				if ( t > d ) d = t;
				t = fabs( s[kk-1] );
				if ( t > d ) d = t;
				t = fabs( e[kk-1] );
				if ( t > d ) d = t;
				sm = s[mm-1]/d; sm1 = s[mm-2]/d;
				em1 = e[mm-2]/d;
				sk = s[kk-1]/d; ek = e[kk-1]/d;
				b = ((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
				c = sm*em1; c = c*c; shh = 0.0;
				if ( (b!=0.0) || (c!=0.0) )
				{ 
					shh = sqrt( b*b+c );
					if ( b < 0.0 ) shh = -shh;
					shh = c/(b+shh);
				}
				fg[0] = (sk+sm)*(sk-sm)-shh;
				fg[1] = sk*ek;
				for ( i = kk; i <= mm-1; i++ )
				{ 
					sss( fg, cs );
					if ( i != kk ) e[i-2] = fg[0];
					fg[0] = cs[0]*s[i-1]+cs[1]*e[i-1];
					e[i-1] = cs[0]*e[i-1]-cs[1]*s[i-1];
					fg[1] = cs[1]*s[i];
					s[i] = cs[0]*s[i];
					if ( (cs[0] != 1.0) || (cs[1] != 0.0) )
						for ( j = 1; j <= n; j++ )
						{ 
							ix = (j-1)*n+i-1;
							iy = (j-1)*n+i;
							d = cs[0]*v[ix]+cs[1]*v[iy];
							v[iy] = -cs[1]*v[ix]+cs[0]*v[iy];
							v[ix] = d;
						}
						sss( fg, cs );
						s[i-1] = fg[0];
						fg[0] = cs[0]*e[i-1]+cs[1]*s[i];
						s[i] = -cs[1]*e[i-1]+cs[0]*s[i];
						fg[1] = cs[1]*e[i];
						e[i] = cs[0]*e[i];
						if ( i < m )
							if ( ( cs[0] != 1.0 ) || ( cs[1] != 0.0 ) )
								for ( j = 1; j <= m; j++ )
								{ 
									ix = (j-1)*m+i-1;
									iy = (j-1)*m+i;
									d = cs[0]*u[ix]+cs[1]*u[iy];
									u[iy] = -cs[1]*u[ix]+cs[0]*u[iy];
									u[ix] = d;
								}
				}
				e[mm-2] = fg[0];
				it = it-1;
			}
			else
			{ 
				if ( ks == mm )
				{ 
					kk = kk+1;
					fg[1] = e[mm-2]; e[mm-2] = 0.0;
	                for ( ll = kk; ll <= mm-1; ll++ )
		            { 
						i = mm+kk-ll-1;
				        fg[0] = s[i-1];
					    sss( fg, cs );
						s[i-1] = fg[0];
					    if ( i != kk )
						{ 
							fg[1] = -cs[1]*e[i-2];
							e[i-2] = cs[0]*e[i-2];
						}
						if ( (cs[0]!=1.0)||(cs[1]!=0.0) )
							for ( j = 1; j <= n; j++ )
							{ 
								ix = (j-1)*n+i-1;
								iy = (j-1)*n+mm-1;
								d = cs[0]*v[ix]+cs[1]*v[iy];
								v[iy] = -cs[1]*v[ix]+cs[0]*v[iy];
								v[ix] = d;
							}
					}
				}
				else
				{ 
					kk = ks+1;
					fg[1] = e[kk-2];
					e[kk-2] = 0.0;
					for ( i = kk; i <= mm; i++ )
					{ 
						fg[0] = s[i-1];
						sss( fg, cs );
						s[i-1] = fg[0];
						fg[1] = -cs[1]*e[i-1];
						e[i-1] = cs[0]*e[i-1];
						if ( (cs[0]!=1.0)||(cs[1]!=0.0) )
							for (j=1; j<=m; j++)
							{ 
								ix = (j-1)*m+i-1;
								iy = (j-1)*m+kk-2;
								d = cs[0]*u[ix]+cs[1]*u[iy];
								u[iy] = -cs[1]*u[ix]+cs[0]*u[iy];
								u[ix] = d;
							}
					}
				}
			}
		}
    }

    return( 1 );
}

# define EPS  0.000001

/*
  PCA2: Perform PCA using SVD.
  data     --- MxN matrix of input data ( M dimensions, N trials )
  signals  --- MxN matrix of projected data 
  PC       --- each column is a PC
  V        --- Mx1 matrix of variances
  row = M dimensions, col = N trials 
*/
int pca2( double * data, int row, int col,
	      double * signals, double * PC, double * V )
{
	int x, y, i, ka, rvalue;
	double * mean;
	double * u;
    double * d;
    double * v;
	double * sd;
	double * data1;
	double sqrt_col_1;
	int col1, row1;

	if ( row <= 1 || col <= 1 ) return( -1 );

	/* subtract off the mean for each dimension */
	mean = ( double * )malloc( sizeof( double )*row );
	data1 = ( double * ) malloc( sizeof( double )*row*col );
	for ( y = 0; y < row; y++ ) /* calculate mean */
	{
		mean[y] = 0;
		for ( x = 0; x < col; x++ )
			mean[y] += data[y*col+x];
	}
	for ( y = 0; y < row; y++ ) mean[y] = mean[y]/col;
	for ( y = 0; y < row; y++ ) /* subtract mean */
		for ( x = 0; x < col; x++ )
		{
			data1[y*col+x] = data[y*col+x]-mean[y];
		}
	free( mean );

	/* construct the matrix Y: Y = data' / sqrt(col-1); */

	/* construct the matrix Y: Y = data' / sqrt(col-1); */
	sqrt_col_1 = sqrt(col-1);
	row1 = col;
	col1 = row;
	sd = (double *)malloc( sizeof(double)*col1*row1 );
	u = (double *)malloc( sizeof(double)*row1*row1 );
	v = (double *)malloc( sizeof(double)*col1*col1 );
	for ( y = 0; y < row1; y++ )
		for ( x = 0; x < col1; x++ )
		{
			sd[y*col1+x] = data1[x*row1+y]/sqrt_col_1;
		}

	/* SVD does it all: [u, S, PC] = svd( Y ); */
	if ( row1 >= col1 ) 
		ka = row1+1;
	else 
		ka = col1+1;
	rvalue = svd( sd, row1, col1, u, v, EPS, ka ); /* svd decomposition */
	d = (double *)malloc( sizeof(double) * col1 );
	for ( i = 0; i < col1; i++ ) d[i] = 0;
	if ( row1 <= col1 )
	{
		for ( i = 0; i < row1; i++ )
			d[i] = sd[i*col1+i];
	}
	else
	{
		for ( i = 0; i < col1; i++ )
			d[i] = sd[i*col1+i];
	}

	/* calculate the variances: S = diag( S ); V = S .* S; */
	for ( x = 0; x < col1; x++ )
		V[x] = d[x] * d[x];
	/* get PC */
	for ( y = 0; y < col1; y++ )
		for ( x = 0; x < col1; x++ )
		{
			PC[y*col1+x] = v[x*col1+y];
		}

	/* project the original data: signals = PC' * data; */
	for ( y = 0; y < row; y++ )
		for ( x = 0; x < col; x++ )
		{
			signals[y*col+x] = 0;
			for ( i = 0; i < row; i++ )
				signals[y*col+x] += PC[i*row+y] * data1[i*col+x];
		}

	free( u );
	free( v );
	free( sd );
	free( d );
	free( data1 );

	return( 1 );
}

/*
   main function
*/
int main()
{
	double data[20] = { 2, 3, 8, 2, 3,
		               7, 9, 29, 3, 5,
					   3, 8, 22, 12, 1,
					   3, 12, 12, 33, 2};
	int row = 4, col = 5;
	double signals[20], PC[16], V[4];
	int x, y;

	pca2( data, row, col, signals, PC, V );

	printf( "Project to Principal Component: \n" );
	for ( y = 0; y < row; y++ )
	{
		for ( x = 0; x < col; x++ )
		{
			printf( "%7.4f ", signals[y*col+x] );
		}
		printf( "\n" );
	}
	printf( "Eigen Values (in deceasing order): \n" );
	for ( y = 0; y < row; y++ )
		printf( "%7.4f ", V[y] );
	printf( "\n" );
	printf( "Eigen Vectors: \n" );
	for ( y = 0; y < row; y++ )
	{
		for ( x = 0; x < row; x++ )
		{
			printf( "%7.4f ", PC[y*row+x] );
		}
		printf( "\n" );
	}

	return( 1 );
}