clda.cpp

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

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

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

CLda::~CLda()
{
}

bool CLda::extract(char* trainFile, char* newTrainFile, char* testFile, char* newTestFile)
{
	bool ret = true;
	int i, j, k, index, cnum;
	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[numClass*numFeature];
	DOUBLE *sumXX = new DOUBLE[numClass*numFeature*numFeature];
	int* classSample = new int[numClass];
	int total;
	for (i=0; i<numClass; i++)
	{
		classSample[i]=0;
	}
	total = numClass * numFeature;
	for (i=0; i<total; i++)
	{
		sumX[i]=0;
	}
	total = numClass * numFeature * numFeature;
	for (i=0; i<total; i++)
	{
		sumXX[i]=0;
	}
	for (i=0; i<numSample; i++)
	{
		cnum = this->sdata.ydata[i];
		classSample[cnum] ++;
		for (j=0; j<numFeature; j++)
		{
			sumX[cnum*numFeature+j] += this->sdata.xdata[i*numFeature+j];
			for (k=j; k<numFeature; k++)
			{
				sumXX[(cnum*numFeature+j)*numFeature+k] +=this->sdata.xdata[i*numFeature+j] * this->sdata.xdata[i*numFeature+k];
			}
		}
	}
	for (i=0; i<numClass; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			sumX[i*numFeature+j] = sumX[i*numFeature+j] / ((DOUBLE) classSample[i]);
		}
	}
	for (i=0; i<numClass; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			for (k=j; k<numFeature; k++)
			{
				index = (i*numFeature+j)*numFeature+k;
				sumXX[index] = sumXX[index]/((DOUBLE)classSample[i])-(sumX[i*numFeature+j]*sumX[i*numFeature+k]);
			}
		}
	}
	DOUBLE *U0 = new DOUBLE[numFeature];
	DOUBLE *SW = new DOUBLE[numFeature*numFeature];
	DOUBLE *SB = new DOUBLE[numFeature*numFeature];
	for (j=0; j<numFeature; j++)
	{
		U0[j] = 0;
		for (k=0; k<numClass; k++)
		{
			U0[j] += sumX[k*numFeature+j]*((DOUBLE)classSample[k]);
		}
		U0[j] = U0[j]/((DOUBLE)numSample);
	}
	for (i=0; i<numFeature; i++)
	{
		for (j=i; j<numFeature; j++)
		{
			index = i*numFeature+j;
			SW[index] = 0;
			for (k=0; k<numClass; k++)
			{
				SW[index] += sumXX[((k*numFeature)+i)*numFeature+j]*classSample[k];
			}
			SW[index] = SW[index]/((DOUBLE)numSample);
		}
	}
	for (i=0; i<numFeature; i++)
	{
		for (j=i; j<numFeature; j++)
		{
			index = i*numFeature+j;
			SB[index] = 0;
			for (k=0; k<numClass; k++)
			{
				SB[index] += (sumX[k*numFeature+i]-U0[i])*(sumX[k*numFeature+j]-U0[j])*((DOUBLE)classSample[k]);
			}
			SB[index] = SB[index]/((DOUBLE)numSample);
		}
	}
	for (j=1; j<numFeature; j++)
	{
		for (k=0; k<j; k++)
		{
			SW[j*numFeature+k] = SW[k*numFeature+j];
			SB[j*numFeature+k] = SB[k*numFeature+j];
		}
	}
	/*for (j=0; j<numFeature; j++)//限制协方差矩阵的对角元素的下限
	{
		if (SW[j*numFeature+j]<0.001)
		{
			SW[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)SW[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 += SW[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 += SW[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;
	}
	DOUBLE* W1 = new DOUBLE[numFeature*numFeature];//计算W1
	for (i=0; i<numFeature; i++)
	{
		if (featureValue[i]<0.001)//限制特征值最小值为0.001
		{
			temp1 = sqrt(0.001);
		}
		else
		{
			temp1 = sqrt(featureValue[i]);
		}
		for (j=0; j<numFeature; j++)
		{
			W1[j*numFeature+i] = ((DOUBLE)P[j][i])/temp1;
		}
	}
	for (i=0; i<numFeature; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			sumXX[i*numFeature+j] = 0;
			for (k=0; k<numFeature; k++)
			{
				sumXX[i*numFeature+j] += W1[k*numFeature+i]*SB[k*numFeature+j];
			}
		}
	}
	for (i=0; i<numFeature; i++)
	{
		for (j=i; j<numFeature; j++)
		{
			SB[i*numFeature+j] = 0;
			for (k=0; k<numFeature; k++)
			{
				SB[i*numFeature+j] += sumXX[i*numFeature+k]*W1[k*numFeature+j];
			}
		}
	}
	for (i=1; i<numFeature; i++)
	{
		for (j=0; j<i; j++)
		{
			SB[i*numFeature+j] = SB[j*numFeature+i];
		}
	}
	cout<<endl;
	for (i=0; i<numFeature; i++)
	{
		cout<<endl;
		for (j=0; j<numFeature; j++)
		{
			Cova[i][j] = (float)SB[i*numFeature+j];
			cout<<Cova[i][j]<<"  ";
		}
	}
	cout<<endl;
	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 += SB[i*numFeature+k]*((DOUBLE)P[k][j]);
			}
			cout<<temp2/P[i][j]<<"  ";
		}
	}
	cout<<endl;

	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 += SB[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* W = new DOUBLE[numFeature*numFeature];//计算W
	for (i=0; i<numFeature; i++)
	{
		for (j=0; j<numFeature; j++)
		{
			W[i*numFeature+j] = 0;
			for (k=0; k<numFeature; k++)
			{
				W[i*numFeature+j] += W1[i*numFeature+k]*sumXX[k*numFeature+j];
			}
		}
	}
	//提取特征
	DOUBLE *temp = new DOUBLE[newNumF];
	for (i=0; i<numSample; i++)//中心化
	{
		for (j=0; j<numFeature; j++)
		{
			this->sdata.xdata[i*numFeature+j] -= U0[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] += (W[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] -= U0[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] += (W[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[] Cova;
	delete[] P;
	delete[] temp;
	delete[] featureValue;
	delete[] featureIndex;
	delete[] W1;
	delete[] W;
	delete[] sumX;
	delete[] sumXX;
	delete[] U0;
	delete[] SW;
	delete[] SB;
	return ret;
}