cpca.cpp

该程序包实现了模式识别中的两个特征提取算法

#include "stdlib.h"
#include <iostream>
#include <string.h>
#include <fstream>
#include <math.h>
using namespace std;
#include "global.h"

CPca::CPca()
{
	this->numF = -1;
}

CPca::~CPca()
{
}

bool CPca::extract(char* trainFile, char* newTrainFile, char* testFile, char* newTestFile)
{
	bool ret = true;
	int i, j, k, index;
	if (this->numF <= 0)
		this->numF = 1;
	ret = this->sdata.readFile(trainFile);//读入训练数据
	if (!ret)
		return false;//如果数据格式不正确，退出程序。
	int numClass = this->sdata.numClass;
	int numFeature = this->sdata.numFeature;
	int numSample = this->sdata.numSample;
	CSampleData testData;
	ret = testData.readFile(testFile);
	if (!ret)
		return false;//如果数据格式不正确，退出程序。
	DOUBLE temp1, temp2;
	if (this->numF > this->sdata.numFeature)
		this->numF = this->sdata.numFeature;
	int newNumF = this->numF;
	DOUBLE *sumX = new DOUBLE[numFeature];
	DOUBLE *sumXX = new DOUBLE[numFeature*numFeature];
	int total;
	for (i=0; i<numFeature; i++)
	{
		sumX[i]=0;
	}
	total = numFeature * numFeature;
	for (i=0; i<total; i++)
	{
		sumXX[i]=0;
	}
	for (i=0; i<numSample; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			sumX[j] += this->sdata.xdata[i*numFeature+j];
			for (k=j; k<numFeature; k++)
			{
				sumXX[j*numFeature+k] +=this->sdata.xdata[i*numFeature+j] * this->sdata.xdata[i*numFeature+k];
			}
		}
	}
	for (j=0; j<numFeature; j++)
	{
		sumX[j] = sumX[j] / numSample;
	}
	for (j=0; j<numFeature; j++)
	{
		for (k=j; k<numFeature; k++)
		{
			index = j*numFeature+k;
			sumXX[index] = sumXX[index]/numSample-(sumX[j]*sumX[k]);
		}
	}
	for (j=1; j<numFeature; j++)
	{
		for (k=0; k<j; k++)
		{
			sumXX[j*numFeature+k] = sumXX[k*numFeature+j];
		}
	}
	/*for (j=0; j<numFeature; j++)//限制协方差矩阵的对角元素的下限
	{
		if (sumXX[j*numFeature+j]<0.001)
		{
			sumXX[j*numFeature+j] = 0.001;
		}
	}*/
	//计算特征值
	CMatrix theMat;
	float** Cova;
	// Allocate space for the covariance matrix
	Cova = theMat.allocMat(numFeature);		// n: dimensionality of vector
	for (i=0; i<numFeature; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			Cova[i][j] = (float)sumXX[i*numFeature+j];
		}
	}
	// Compute the covariance matrix......
	// Diagonalization
	// P: matrix with columns as the eigenvectors
	float** P;
	P = theMat.allocMat(numFeature);
	theMat.diagonalize( Cova, numFeature, P );
	//显示计算的特征值
	for (i=0; i<numFeature; i++)
	{
		cout<<endl;
		for (j=0; j<numFeature; j++)
		{
			temp2 = 0;
			for (k=0; k<numFeature; k++)
			{
				temp2 += sumXX[i*numFeature+k]*((DOUBLE)P[k][j]);
			}
			cout<<temp2/P[i][j]<<"  ";
		}
	}
	cout<<endl;

	DOUBLE* featureValue = new DOUBLE[numFeature];
	int* featureIndex = new int[numFeature];
	for (i=0; i<numFeature; i++)
	{
		featureValue[i] = 0;
		featureIndex[i] = i;
	}
	for (j=0; j<numFeature; j++)//计算特征值
	{
		temp1 = 0;
		for (i=0; i<numFeature; i++)
		{
			temp2 = 0;
			for (k=0; k<numFeature; k++)
			{
				temp2 += sumXX[i*numFeature+k]*((DOUBLE)P[k][j]);
			}
			if (P[i][j]>0)
			{
				temp1 += (DOUBLE)P[i][j];
				featureValue[j] += temp2;
			}
			else
			{
				temp1 -= (DOUBLE)P[i][j];
				featureValue[j] -= temp2;
			}
		}
		featureValue[j] = featureValue[j]/temp1;
		cout<<featureValue[j]<<endl;
	}
	for (i=0; i<numFeature-1; i++)
	{
		for (j=numFeature-1; j>i; j--)
		{
			if (featureValue[j-1]<featureValue[j])
			{
				temp1 = featureValue[j-1];
				featureValue[j-1] = featureValue[j];
				featureValue[j] = temp1;
				k = featureIndex[j-1];
				featureIndex[j-1] = featureIndex[j];
				featureIndex[j] = k;
			}
		}
	}
	for (i=0; i<numFeature; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			sumXX[i*numFeature+j] = (DOUBLE)P[i][featureIndex[j]];
		}
	}
	//提取特征
	DOUBLE *temp = new DOUBLE[newNumF];
	for (i=0; i<numSample; i++)//中心化
	{
		for (j=0; j<numFeature; j++)
		{
			this->sdata.xdata[i*numFeature+j] -= sumX[j];
		}
	}
	for (i=0; i<numSample; i++)
	{
		for (j=0; j<newNumF; j++)
		{
			temp[j] = 0;
		}
		for (j=0; j<newNumF; j++)
		{
			for (k=0; k<numFeature; k++)
			{
				temp[j] += (sumXX[k*numFeature+j])*this->sdata.xdata[i*numFeature+k];
			}
		}
		for (j=0; j<newNumF; j++)
		{
			this->sdata.xdata[i*newNumF+j] = temp[j];
		}
	}
	this->sdata.numFeature = newNumF;
	this->sdata.writeFile(newTrainFile);//写新的训练数据文件
	
	for (i=0; i<testData.numSample; i++)//中心化测试数据
	{
		for (j=0; j<numFeature; j++)
		{
			testData.xdata[i*numFeature+j] -= sumX[j];
		}
	}
	for (i=0; i<testData.numSample; i++)
	{
		for (j=0; j<newNumF; j++)
		{
			temp[j] = 0;
		}
		for (j=0; j<newNumF; j++)
		{
			for (k=0; k<numFeature; k++)
			{
				temp[j] += (sumXX[k*numFeature+j])*testData.xdata[i*numFeature+k];
			}
		}
		for (j=0; j<newNumF; j++)
		{
			testData.xdata[i*newNumF+j] = temp[j];
		}
	}
	testData.numFeature = newNumF;
	testData.writeFile(newTestFile);
	delete[] featureValue;
	delete[] featureIndex;
	delete[] Cova;
	delete[] P;
	delete[] temp;
	delete[] sumX;
	delete[] sumXX;
	return ret;
}