⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 qrmethod.cpp

📁 mimo多天线blast系统qr算法在matlab上的实现
💻 CPP
字号:
🔍 
#include "math.h"
#include "iostream.h"
#include "stdio.h"
#include <iomanip.h>

#define BOOL int
#define FALSE 0
#define TRUE 1

double ee = 0.00000000001;  //精度水平ee=1e-12;

//显示矩阵
void ShowMatrix(double * P_Matrix,int N)
{
	cout<<endl;
	for(int i = 0; i < N; i++)
	{
		cout<<endl<<i+1<<" row numbers are:"<<endl;
		for(int j = 0; j < N; j++)
		{
	 		cout<<setiosflags(ios::scientific)<<setprecision(11)<<P_Matrix[i*N+j]<<"  ";
		}
	}
	
}

//显示列向量
void ShowVector(double * P_Vector,int N)
{
	cout<<endl;
	
	for(int i = 0; i < N; i++)
	{
		cout<<setiosflags(ios::scientific)<<setprecision(11)<<P_Vector[i]<<endl;
	}	
}

//显示所有特征值Lmuda
void ShowVectorLumda(double * pf_LmudaRe,double * pf_LmudaIm,int N)
{
	cout<<endl;
	
	for(int i = 0; i < N; i++)
	{
		if(pf_LmudaIm[i] == 0)
		{
			cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[i]<<endl;
		}
		if(pf_LmudaIm[i] < 0)
		{
			cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[i]<<" - "<<-pf_LmudaIm[i]<<" i"<<endl;
		}
		if(pf_LmudaIm[i] > 0)
		{
			cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[i]<<" + "<<pf_LmudaIm[i]<<" i"<<endl;
		}
	}
	
	
}

//矩阵乘法——相同维数矩阵相乘
double * Multiply(double * pf_MatraixA, double * pf_MatraixB, int N)
{
	int i,j,k;
	double * pf_Result = new double[N*N];
	
	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			pf_Result[i*N+j] = 0;
			for(k = 0; k < N; k++)
			{
				pf_Result[i*N+j] += pf_MatraixA[i*N+k]*pf_MatraixB[k*N+j]; 
			}
		}
	}
	
	return pf_Result;

}


//矩阵乘法——矩阵乘以列向量
double * MultiplyMV(double * pf_MatraixA, double * pf_Vector, int N)
{
	int i,k;
	double * pf_Result = new double[N];
	
	for(i = 0; i < N; i++)
	{
		pf_Result[i] = 0;
		for(k = 0; k < N; k++)
		{
			pf_Result[i] += pf_MatraixA[i*N+k]*pf_Vector[k]; 
		}
	}
	
	return pf_Result;

}

//矩阵乘法——数乘以矩阵
double * MultiplyNumM(double * pf_MatraixA, double num, int N)
{
	int i,j;
	double * pf_Result = new double[N*N];
	
	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			pf_Result[i*N+j] = pf_MatraixA[i*N+j] * num; 
		}
	}
	
	return pf_Result;
	
}

//行向量乘以列向量
double MultiplyRC(double * pf_VectorA, double * pf_VectorB, int N)
{
	int i;
	double temp = 0;

	for(i = 0; i < N; i++)
	{
			temp += pf_VectorA[i]*pf_VectorB[i]; 
	}
	
	return temp;

}

//列向量乘以行向量
double * MultiplyCR(double * pf_VectorA, double * pf_VectorB, int N)
{
	int i,j;
	double temp = 0;
	double * pf_Result = new double[N*N];
	
	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			pf_Result[i*N+j] = pf_VectorA[i]*pf_VectorB[j];
		}
	}
	
	return pf_Result;
	
}


//矩阵转置
double * MatrixT(double * pf_Matraix, int N)
{
	int i,j;
	double * pf_MatraixT = new double[N*N];

	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			pf_MatraixT[i*N+j] = pf_Matraix[j*N+i];
		}
	}

	return pf_MatraixT;
	
}

//向量乘以数
double * VectorMultiNum(double * pf_Vector, double num, int N)
{
	int i;

	double * pf_Result = new double[N];

	for(i = 0; i < N; i++)
	{
		pf_Result[i] = pf_Vector[i]*num;
	}

	return pf_Result;
}

//向量相减
double * VectorSub(double * pf_VectorA, double * pf_VectorB, int N)
{
	int i;
	double temp;
	double * pf_Result = new double[N];

	for(i = 0; i < N; i++)
	{
		pf_Result[i] = pf_VectorA[i] - pf_VectorB[i];
	}

	return pf_Result;
}


//矩阵相减
BOOL MatrixSub(double * pf_MatrixA, double * pf_MatrixB, double * pf_Result, int N)
{
	int i,j;

	for(i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			pf_Result[i*N+j] = pf_MatrixA[i*N+j] - pf_MatrixB[i*N+j];
		}
	}
	
	return TRUE;
}



//对矩阵第r列进行上三角变换
BOOL HessOneStep(double * pf_Matraix,int r, int N)
{
	int i;
	double dr,cr,hr,temp;

	temp = 0;
	for(i = r+1; i < N; i++)
	{
		temp += pf_Matraix[i*N+r]*pf_Matraix[i*N+r];
	}

	dr = sqrt(temp);
	
	if (pf_Matraix[(r+1)*N+r] <= 0)
	{
		cr = dr;
	}
	else
	{
		cr = -dr;
	}

	hr = cr*cr - cr*pf_Matraix[(r+1)*N+r];

	double * ur = new double [N];
	double * pr;
	double * qr;
	double tr = 0;
	double * wr;

	for(i = 0; i < N; i++)
	{
		if (i <= r) 
		{
			ur[i] = 0;
		}
		
		if (i == r+1)
		{
			ur[i] = pf_Matraix[i*N+r] - cr;
		}
		if( i > r+1)
		{
			ur[i] = pf_Matraix[i*N+r];
		}
	}
	
	double * pf_MatraixT;

	pf_MatraixT = MatrixT(pf_Matraix,N);

	pr = VectorMultiNum(MultiplyMV(pf_MatraixT,ur,N),1/hr,N);
	qr = VectorMultiNum(MultiplyMV(pf_Matraix,ur,N),1/hr,N);
	tr = MultiplyRC(pr,ur,N)/hr;
	wr = VectorSub(qr,VectorMultiNum(ur,tr,N),N);

	MatrixSub(pf_Matraix,MultiplyCR(wr,ur,N),pf_Matraix,N);
	MatrixSub(pf_Matraix,MultiplyCR(ur,pr,N),pf_Matraix,N);

	delete ur;
	delete pr;
	delete qr;
	delete wr;
	delete pf_MatraixT;
	return TRUE;
}

//矩阵拟上三角化
double * Hessenberg(double * pf_Matraix,int N)
{
	double * pf_HessMatraix = new double[N*N];
	int i,j;
	
	//矩阵：A1 = A
	for(i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			pf_HessMatraix[i*N+j] = pf_Matraix[i*N+j];
		}
	}	

	int r;
	BOOL flag;

	for(r = 0; r < N-2; r++)
	{
		flag = FALSE;
		for(i = r+2; i < N; i++)
		{
			if (fabs(pf_HessMatraix[i*N+r]) > ee)
			{
				flag = TRUE;
				break;
			}
		}
		HessOneStep(pf_HessMatraix,r,N);
	}
	for(i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			if (fabs(pf_HessMatraix[i*N+j]) <= ee)
			{
				pf_HessMatraix[i*N+j] = 0;
			}
		}
	}		
	return pf_HessMatraix;
}

void QRDDStep4(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count);
void QRDDStep567(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count);
void QRDDStep8910(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count);
void QRDDStep11(int& m);
void QROneStep(double * pf_Matraix,double * pf_MatraixM,int r,int N);

//step(3)
void QRDDStep3(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
	if (fabs(pf_Matraix[(m-1)*N+m-2]) <= ee)	
	{
		pf_LmudaRe[count] = pf_Matraix[(m-1)*N+m-1];
		count++;
		m = m-1;
		QRDDStep4(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);

	}
	else
	{
		QRDDStep567(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
	}
	
}

//step(4)
void QRDDStep4(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
	if (m == 1)	
	{
		pf_LmudaRe[count] = pf_Matraix[0];
		count++;
		QRDDStep11(m);	
		
	}
	if(m == 0)
	{
		QRDDStep11(m);
	}
	
	if (m > 1)
	{
		QRDDStep3(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
	}
	
}

void GetTwoLumda(double * pf_LmudaRe,double * pf_LmudaIm,double b,double c,int& count)
{
	double Delta;	//delta = b*b - 4ac
	
	Delta = b*b-4*c;
		
	if ( Delta >= 0)
	{
		pf_LmudaRe[count] = (-b + sqrt(Delta))/2;
		count++;		
		pf_LmudaRe[count] = (-b - sqrt(Delta))/2;
		count++;		
	}
	else
	{
		pf_LmudaRe[count] = -b/2;
		pf_LmudaIm[count] = sqrt(-Delta)/2;
		count++;		
		
		pf_LmudaRe[count] = -b/2;
		pf_LmudaIm[count] = -sqrt(-Delta)/2;
		count++;		
		
	}
}

//step(5)&&step(6)&&step(7)
void QRDDStep567(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
	double s1,s2;
	int i;

	double * pf_MatraixD;

	pf_MatraixD = new double[4];

	pf_MatraixD[0] = pf_Matraix[(m-2)*N+m-2];
	pf_MatraixD[1] = pf_Matraix[(m-2)*N+m-1];
	pf_MatraixD[2] = pf_Matraix[(m-1)*N+m-2];
	pf_MatraixD[3] = pf_Matraix[(m-1)*N+m-1];	

	double b,c;	//二阶方程的系数b,c

	b = -(pf_MatraixD[0]+pf_MatraixD[3]);
	c = ((pf_MatraixD[0] * pf_MatraixD[3]) - (pf_MatraixD[1] * pf_MatraixD[2]));

	if (m == 2)	
	{
		GetTwoLumda(pf_LmudaRe,pf_LmudaIm,b,c,count);
		QRDDStep11(m);
	}
	else
	{
		if (fabs(pf_Matraix[(m-2)*N+m-3]) <= ee)
		{
			GetTwoLumda(pf_LmudaRe,pf_LmudaIm,b,c,count);			
			m = m - 2;
			QRDDStep4(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
		}
		else
		{
			QRDDStep8910(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
		}				
	}
}

//step(8)&&step(9)&&step(10)
void QRDDStep8910(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
	int i,j,r;
	
	double * pf_MatraixA = new double[m*m];
	double * pf_MatraixM = new double[m*m];
	
	for(i = 0; i < m; i++)
	{
		for(j = 0; j < m; j++)
		{
			pf_MatraixA[i*m+j] = pf_Matraix[i*N+j];
		}
	}

	double s,t;

	s = pf_Matraix[(m-2)*N+m-2] + pf_Matraix[(m-1)*N+m-1];
	t = (pf_Matraix[(m-2)*N+m-2] * pf_Matraix[(m-1)*N+m-1])
		 - (pf_Matraix[(m-1)*N+m-2] * pf_Matraix[(m-2)*N+m-1]);

	double * pf_MatraixI= new double[m*m];		//单位矩阵	
	
	//单位矩阵赋值
	for(i = 0; i < m; i++)
	{
		for(j = 0; j < m; j++)
		{
			if (i == j)
			{
				pf_MatraixI[i*m+j] = 1;
			}
			else
			{
				pf_MatraixI[i*m+j] = 0;
			}
			
		}
	}
	

	
	MatrixSub(Multiply(pf_MatraixA,pf_MatraixA,m),MultiplyNumM(pf_MatraixA,s,m),pf_MatraixM,m);	//M1 = A^2-sA
	MatrixSub(pf_MatraixM,MultiplyNumM(pf_MatraixI,-t,m),pf_MatraixM,m);						//M = M1-(-tI)

	BOOL flag;
	
	for(r = 0; r < m - 1; r++)
	{
		flag = FALSE;
		for(i = r+1; i < m; i++)
		{
			if (fabs(pf_MatraixM[i*m+r]) > ee)
			{
				flag = TRUE;
				break;
			}
		}
		QROneStep(pf_MatraixM,pf_MatraixA,r,m);
	}

	//矩阵A的前m*m元素赋值为QR分解后的矩阵
	for(i = 0; i < m; i++)
	{
		for(j = 0; j < m; j++)
		{
			pf_Matraix[i*N+j] = pf_MatraixA[i*m+j];
		}
	}
	
	delete  pf_MatraixM;
	delete  pf_MatraixI;
	delete  pf_MatraixA;
	
	
}

//step(11)
void QRDDStep11(int& m)
{
	m = 0;
}

//用具体简化算法进行QR分解
void QROneStep(double * pf_MatraixB,double * pf_MatraixC,int r,int m)
{
	int i;
	double dr,cr,hr,temp;
	
	temp = 0;
	for(i = r; i < m; i++)
	{
		temp += pf_MatraixB[i*m+r]*pf_MatraixB[i*m+r];
	}	
	dr = sqrt(temp);
	
	if (pf_MatraixB[r*m+r] <= 0)
	{
		cr = dr;
	}
	else
	{
		cr = -dr;
	}
	
	hr = cr*cr - cr*pf_MatraixB[r*m+r];
	
	double * ur = new double [m];
	double * pr;
	double * qr;
	double * vr;
	double tr = 0;
	double * wr;

	for(i = 0; i < m; i++)
	{
		if (i < r) 
		{
			ur[i] = 0;
		}
		
		if (i == r)
		{
			ur[i] = pf_MatraixB[i*m+r] - cr;
		}
		if( i > r)
		{
			ur[i] = pf_MatraixB[i*m+r];
		}
	}

	vr = VectorMultiNum(MultiplyMV(MatrixT(pf_MatraixB,m),ur,m),1/hr,m);
	MatrixSub(pf_MatraixB,MultiplyCR(ur,vr,m),pf_MatraixB,m);	//求B(r+1)
	pr = VectorMultiNum(MultiplyMV(MatrixT(pf_MatraixC,m),ur,m),1/hr,m);
	qr = VectorMultiNum(MultiplyMV(pf_MatraixC,ur,m),1/hr,m);
	tr = MultiplyRC(pr,ur,m) / hr;
	wr = VectorSub(qr,VectorMultiNum(ur,tr,m),m);

	MatrixSub(pf_MatraixC,MultiplyCR(wr,ur,m),pf_MatraixC,m);	//求C(r+1)
	MatrixSub(pf_MatraixC,MultiplyCR(ur,pr,m),pf_MatraixC,m);	//求C(r+1)

	delete ur;
	delete pr;
	delete qr;
	delete wr;
	delete vr;
}

//对拟上三角化的矩阵进行QR分解求特征值
BOOL QRDoubleDispace(double * pf_Matraix,int N, double * pf_LmudaRe,double * pf_LmudaIm)
{
	int i,k,m,count;

	count = 0;
	m = N;
	
	while (m != 0)
	{
		QRDDStep3(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
		
	}	
	for(i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			if (fabs(pf_Matraix[i*N+j]) <= ee)
			{
				pf_Matraix[i*N+j] = 0;
			}
		}
	}		
	return TRUE;
}

//对矩阵进行选主元素高斯消去
BOOL GaussRe(double * P_Matrix,double * P_MatrixG, int N)
{
	double m,temp;
	
	int i,j,k,l;

	int flag = 0;
	
	if (P_Matrix[0] == 0)
	{
		return FALSE;
	}

	//复制原始矩阵和向量
	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			P_MatrixG[i*N+j] = P_Matrix[i*N+j];
		}
	}
	//开始消去法
	for(k = 0; k < N-1; k++)
	{
		temp = P_MatrixG[k];
		flag = k;
		//找出主元素
		for(l = k; l < N; l++ )
		{
			if (fabs(P_MatrixG[l*N+k]) > fabs(temp)) 
			{
				temp = P_MatrixG[l*N+k];
				flag = l;
			}
		}
		//交换行
		if (flag != k)
		{
			for(l = k; l < N; l++)//交换系数矩阵元素
			{				
				temp = P_MatrixG[flag*N+l];
				P_MatrixG[flag*N+l] = P_MatrixG[k*N+l];
				P_MatrixG[k*N+l] = temp;
			}		
		}				
		for(i = k+1; i < N; i++)
		{			
			m = P_MatrixG[i*N+k]/P_MatrixG[k*N+k];
			P_MatrixG[i*N+k] = 0;	//消去的元素为0			
			for(j = k+1; j < N; j++)
			{
				P_MatrixG[i*N+j] = P_MatrixG[i*N+j]-m*P_MatrixG[k*N+j];
			}
		}
	}
	return TRUE;
}

//回代求解
BOOL Calculate(double * P_MatrixG, double * P_VectorX, int N)
{
	int j,k;
	double temp;

	P_VectorX[N-1] = 1;	//因为是齐次线性方程组，将最后一个x直接赋值为1
	for(k = N-2; k >= 0; k--)
	{
		temp = 0;
		for(j = k+1; j < N; j++)
			{
				temp += P_MatrixG[k*N+j]*P_VectorX[j];
			}
		P_VectorX[k] =  -temp/P_MatrixG[k*N+k];
	}	
	return TRUE;	
}

//求出每个实特征值的特征向量
void CalVectors(double * pf_Matrix,int N, double * pf_LmudaRe,double * pf_LmudaIm)
{
	int i,j,k;
	double * pf_MatrixASubLmuda = new double[N*N];
	double * pf_VectorX = new double[N];

	for (k = 0; k < N; k++)
	{
		if (pf_LmudaIm[k] != 0)	//若特征值不是实数，不进行计算
		{
			continue;
		}

		cout<<endl<<"This is the vector of "<<k+1<<" (";
		cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[k]<<"):";
		for(i = 0; i < N; i++)
		{
			for(j = 0; j < N; j++)
			{
				if (i == j)
				{
					pf_MatrixASubLmuda[i*N+j] = pf_Matrix[i*N+j] - pf_LmudaRe[k];
				}
				else
				{
					pf_MatrixASubLmuda[i*N+j] = pf_Matrix[i*N+j];
				}
			}
		}		
		double * pf_MatrixG = new double[N*N];		
		GaussRe(pf_MatrixASubLmuda,pf_MatrixG,N);			
		Calculate(pf_MatrixG,pf_VectorX,N);		
		ShowVector(pf_VectorX,N);		
	}	
}

void main()
{
	int N;
	int i,j;

	N = 10;
	
	double * pf_MatrixA = new double[N*N];
 	double * pf_Result; 
	double * pf_LmudaIm = new double[N]; 
	double * pf_LmudaRe = new double[N]; 
	
	//Lumda虚部为0
	for(i = 0; i < N; i++)
	{
		pf_LmudaIm[i] = 0;
	}

	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			if (i == j)
			{
				pf_MatrixA[i*N+j] = 1.5*cos(i+1+1.2*(j+1));
			}
			else
			{
				pf_MatrixA[i*N+j] = sin(0.5*(i+1)+0.2*(j+1));
			}
		}
	}
		pf_Result = Hessenberg(pf_MatrixA,N);
		cout<<"The hessenberg matrix is:"<<endl;
		ShowMatrix(pf_Result,N);
		
 		QRDoubleDispace(pf_Result,N,pf_LmudaRe,pf_LmudaIm);
		cout<<endl<<"The matrix after QR is :";		
		ShowMatrix(pf_Result,N);

		cout<<endl<<"The Lmudas are :";		
		ShowVectorLumda(pf_LmudaRe,pf_LmudaIm,N);
		
		cout<<endl<<"The X vectors are :";
		CalVectors(pf_MatrixA,N,pf_LmudaRe,pf_LmudaIm);	
}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -