主成分分析PCA源代码（C++，matlab）.zip_matlab）_neighbor1rt_wheatgva_主成分分析pc

/* 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 <stdio.h> # include <stdlib.h> # include <math.h> void ppp( double a[], double e[], double s[], double v[], int m, int n ) { int i,j,p,q; 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 = 0; ll = k; if ( l > k ) ll = l; if ( ll >= 1 ) { for ( kk = 1; kk <= ll; kk++ ) { if ( kk <= k ) { d = 0.0; for ( i = kk; i <= m; i++ ) { ix = (i-1)*n+kk-1; d = d+a[ix]*a[ix]; } s[kk-1] = sqrt(d); if ( s[kk-1] != 0.0 ) { ix = (kk-1)*n+kk-1; if ( a[ix] != 0.0 ) { s[kk-1] = fabs( s[kk-1] ); if ( a[ix] < 0.0 ) s[kk-1] = -s[kk-1]; } for ( i = kk; i <= m; i++ ) { iy = (i-1)*n+kk-1; a[iy] = a[iy]/s[kk-1]; } a[ix] = 1.0+a[ix]; } s[kk-1] = -s[kk-1]; } if ( n >= kk+1 ) { for ( j = kk+1; j <= n; j++ ) { if (( kk <= k ) && ( s[kk-1] != 0.0 ) ) { d = 0.0; for ( i = kk; i <= m; i++ ) { ix = (i-1)*n+kk-1; iy = (i-1)*n+j-1; d = d+a[ix]*a[iy]; } d = -d/a[(kk-1)*n+kk-1]; for ( i = kk; i <= m; i++ ) { ix = (i-1)*n+j-1; iy = (i-1)*n+kk-1; a[ix] = a[ix]+d*a[iy]; } } e[j-1] = a[(kk-1)*n+j-1]; } } if ( kk <= k ) { for ( i = kk; i <= m; i++ ) { ix = (i-1)*m+kk-1; iy = (i-1)*n+kk-1; u[ix] = a[iy]; } } if ( kk <= l ) { d = 0.0; for ( i = kk+1; i <= n; i++ ) d = d + e[i-1]*e[i-1]; e[kk-1] = sqrt(d); if ( e[kk-1] != 0.0 ) { if ( e[kk] != 0.0 ) { e[kk-1] = fabs(e[kk-1]); if ( e[kk] < 0.0 ) e[kk-1] = -e[kk-1]; } for ( i = kk+1; i <= n; i++ ) e[i-1] = e[i-1]/e[kk-1]; e[kk] = 1.0+e[kk]; } e[kk-1] = -e[kk-1]; if (( kk+1 <= m ) && ( e[kk-1] != 0.0 ) ) { for ( i = kk+1; i <= m; i++ ) w[i-1] = 0.0; for ( j = kk+1; j <= n; j++) for ( i = kk+1; i <= m; i++ ) w[i-1] = w[i-1]+e[j-1]*a[(i-1)*n+j-1]; for ( j = kk+1; j <= n; j++ ) for ( i = kk+1; i <= m; i++ ) { ix = (i-1)*n+j-1; a[ix] = a[ix]-w[i-1]*e[j-1]/e[kk]; } } for ( i = kk+1; i <= n; i++ ) v[(i-1)*n+kk-1] = e[i-1]; } } } mm = n; if ( m+1 < n ) mm = m+1; if ( k < n ) s[k] = a[k*n+k]; if ( m < mm ) s[mm-1] = 0.0; if ( l+1 < mm) e[l] = a[l*n+mm-1]; e[mm-1] = 0.0; nn = m; if ( m > n ) nn = n; if ( nn >= k+1 ) { for ( j = k+1; j <= nn; j++ ) { for ( i = 1; i <= m; i++ ) u[(i-1)*m+j-1] = 0.0; u[(j-1)*m+j-1] = 1.0; } } if ( k >= 1) { for (ll = 1; ll <= k; ll++ ) { kk = k-ll+1; iz = (kk-1)*m+kk-1; if ( s[kk-1] != 0.0 ) { if ( nn >= kk+1 ) for ( j = kk+1; j <= nn; j++ ) { d = 0.0; for ( i = kk; i <= m; i++ ) { ix = (i-1)*m+kk-1; iy = (i-1)*m+j-1; d = d+u[ix]*u[iy]/u[iz]; } d = -d; for ( i = kk; i <= m; i++ ) { ix = (i-1)*m+j-1; iy = (i-1)*m+kk-1; u[ix] = u[ix]+d*u[iy]; } } for ( i = kk; i <= m; i++ ) { ix = (i-1)*m+kk-1; u[ix] = -u[ix]; } u[iz] = 1.0+u[iz]; if ( kk-1 >= 1 ) for ( i = 1; i <= kk-1; i++ ) u[(i-1)*m+kk-1] = 0.0; } else { for ( i = 1; i <= m; i++ ) u[(i-1)*m+kk-1] = 0.0; u[(kk-1)*m+kk-1] = 1.0; } } } for ( ll = 1; ll <= n; ll++ ) { kk = n-ll+1; iz = kk*n+kk-1; if ( (kk<=l) && (e[kk-1] != 0.0) ) { for ( j = kk+1; j <= n; j++ ) { d = 0.0; for ( i = kk