利用c++实现主成分分析算法

用c++实现的主成分分析,适合遥感技术应用.

/*
   PCA program （using SVD）
   Written by Y. Bin Mao
   Video and Image Processing and Analysis Group （VIPAG）
   School of Automation NJUST
   Jan. 8 2008

   All rights reserved. （c）

   Here is a matlab description of the algorithm
   % 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
   %
   function [signals PC V] = pca2（ data ）

   [M N] = size（ data ）;

   % subtract off the mean for each dimension
   mn = mean（ data 2 ）;
   data = data - repmat（ mn 1 N ）;

   % construct the matrix Y
   Y = data‘ / sqrt（N-1）;

   % SVD does it all
   [u S PC] = svd（ Y ）;

   % calculate the variances
   S = diag（ S ）;
   V = S .* S;

   % project the original data
   signals = PC‘ * data;
*/

# include 
# include 
# include 

void ppp（ double a[] double e[] double s[] 
		  double v[] int m int n ）
{ 
	int ijpq;
	double d;
	
	if （ m >= n ） i = n;
	else i = m;
	for （ j = 1; j <= i-1; j++ ）
	{ 
		a[（j-1）*n+j-1] = s[j-1];
		a[（j-1）*n+j] = e[j-1];
	}
	a[（i-1）*n+i-1] = s[i-1];
	if （ m < n ） a[（i-1）*n+i] = e[i-1];
	for （ i = 1; i <= n-1; i++ ）
		for （ j = i+1; j <= n; j++ ）
		{ 
			p = （i-1）*n+j-1; q = （j-1）*n+i-1;
			d = v[p]; v[p] = v[q]; v[q] = d;
		}
		return;
}

void sss（ double fg[2] double cs[2] ）
{
	double r d;
    if （ （ fabs（fg[0]） + fabs（fg[1] ） ） == 0.0 ）
	{ 
		cs[0] = 1.0; cs[1] = 0.0; d = 0.0;
	}
    else
	{ 
		d = sqrt（ fg[0]*fg[0]+fg[1]*fg[1] ）;
		if （ fabs（ fg[0] ） > fabs（ fg[1] ） ）
		{
			d = fabs（d）;
			if （ fg[0] < 0.0 ） d = -d;
		}
        if （ fabs（ fg[1] ） >= fabs（ fg[0] ） ）
        { 
			d = fabs（d）;
            if （ fg[1] < 0.0 ） d = -d;
		}
        cs[0] = fg[0]/d; cs[1] = fg[1]/d;
    }
    r = 1.0;
    if （ fabs（ fg[0] ） > fabs（ fg[1] ） ） 
		r = cs[1];
    else
		if （ cs[0] != 0.0 ） r = 1.0/cs[0];
		fg[0] = d; fg[1] = r;
		
		return;
}

/*
   一般实矩阵奇异值分解
   徐士良. 常用算法程序集（C语言描述），第3版. 清华大学出版社. 2004

   double a[m][n] --- 存放mxn维实矩阵A。返回时其对角线给出奇异值（以非递增次序排列）
                      其余元素均为0
   int m int n   --- 实矩阵A的行数和列数
   double u[m][m] --- 返回左奇异向量U
   double v[n][n] --- 返回右奇异向量V‘
   double eps     --- 给定的精度要求
   int ka         --- 其值为max（m n） + 1
   返回值：
           若返回标志值小于0，则表示程序工作失败；
		   若返回标志值大于0，则表示正常返回   
*/
int svd（ double a[] int m int n double u[] double v[]
	     double eps int ka ）
{ 
	int i j k l it ll kk ix iy mm nn iz m1 ks;
    double d dd t sm sm1 em1 sk ek b c shh fg[2] cs[2];
    double * s * e * w;

    s = （ double * ）malloc（ ka * sizeof（double） ）;
    e = （ double * ）malloc（ ka * sizeof（double） ）;
    w = （ double * ）malloc（ ka * sizeof（double） ）;
    it = 60; k = n;
	if （ m-1 < n ） k = m-1;
	l = m;
	if （ n-2 < m ） l = n-2;
	if （ l < 0 ） l 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      16249  2011-08-24 18:41  pca2.cpp

----------- ---------  ---------- -----  ----

                16249                    1