/*
???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?	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]?		}
????????if?（?fabs（?fg[1]?）?>=?fabs（?fg[0]?）?）
????????{?
			d?=?fabs（d）;
????????????if?（?fg[1]?		}
????????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版.?清華大學(xué)出版社.?2004

???double?a[m][n]?---?存放mxn維實矩陣A。返回時其對角線給出奇異值（以非遞增次序排列）
??????????????????????其余元素均為0
???int?m?int?n???---?實矩陣A的行數(shù)和列數(shù)
???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?	l?=?m;
	if?（?n-2?	if?（?l?