matrix.h

实现复数矩阵的基本运算

			for ( int j = 0; j < mCol; j++)
				sum =sum + pow(x(0, j), i+1);
			
			X1(0, i) = sum;
		}

		X2(0, 0) = mCol;
		for( int i = 1; i < n+1; i++)
			X2(0, i) = X1(0, i-1);
		
		for ( int i = 0; i < n; i++)
			for( int j = 0; j < n+1; j++)
				X3(i, j) = X1(0, j+i);

		for ( int i = 0;i < n+1; i++)
			XSUM(0, i) = X2(0, i);

		for ( int i = 1;i < n+1; i++)
			for (int j = 0; j < n+1; j++)
				XSUM(i, j) = X3(i-1, j);

		for( int k = 0; k < n; k++)
		{
			_Ty sum = 0;
			_MatTy xTemp(x);

			for ( int i = 0; i < mCol; i++)
				xTemp(0, i) = pow(x(0, i), k+1);

			xTemp = PM(xTemp, y);

			for ( int j = 0; j < mCol; j++)
				sum = sum + xTemp(0, j);

			YSUM(0, k+1) = sum;
		}

		YSUM(0, 0) = SumNum;

		

		T(YSUM);
		Inv(XSUM);
		coef = XSUM * YSUM;
		T(coef);

		for ( int i = 0,m = n; i < (n+1)/2; i++,n--)
		{
			_Ty Temp = coef(0, i);
			coef(0, i) = coef(0, m);
			coef(0, m) = Temp;
		}
		
		int nc = coef.GetColNum();
		if (nc > 0)
		{
			for ( int i = 0; i < mCol; i++)
				yFit(0, i) = coef(0, 0);
		}	
			

		for ( int i = 1; i < nc; i++)
		{
			yFit = PM(x, yFit) + coef(0, i);
		}
				
		return ;
	}

	//////////////////////////
	// 矩阵运算辅助函数//
	/////////////////////////

	friend _MatTy& Multiply(_MatTy& mOut, const _MatTy& lhs, const _MatTy& rhs)
	{	//断定左边矩阵的列数与右边矩阵的行数相等
		assert(lhs.GetColNum() == rhs.GetRowNum());
		//生成矩阵新对象，用lhs的行作为新阵的行数，用rhs的列数作为新阵的列数
		_MatTy mTmp(lhs.GetRowNum(), rhs.GetColNum());

		for(size_t i = 0; i < mTmp.GetRowNum(); i ++)
		{
			for(size_t j = 0; j < mTmp.GetColNum(); j ++)
			{
				mTmp(i, j) = _Ty(0);		//赋初值0
				for(size_t k = 0; k < lhs.GetColNum(); k ++)
				{
					mTmp(i, j) += lhs(i, k) * rhs(k, j);
				}
			}
		}

		mOut = mTmp;		//将最后结果转放入mOut矩阵中
		return mOut;		//返回结果矩阵mOut
	}


private:
	//矩阵元素
	valarray<_Ty> m_Datas;	

	size_t m_stRow;		//矩阵行数变量
	size_t m_stCol;		//矩阵列数变量
};

//复数共轭
COMPLEX_DOUBLE conjg(COMPLEX_DOUBLE& cc)
{
	double creal,cimag;
	
	creal=cc.real();
	cimag=cc.imag();
	
	return COMPLEX_DOUBLE(creal,-cimag);
	
}

//矩阵共轭
template <class _Ty>
void Conj(CMatrix<_Ty>& mIn)
{
	size_t sR, sC;
	sR = mIn.GetRowNum();			//取原矩阵行数
	sC = mIn.GetColNum();			//取原矩阵列数
	
	//CMatrix<_Ty> mTemp(sC, sR);		//生成一新阵，行数与列数与原阵互换

	for(size_t stC=0; stC<sC; stC++)
		for(size_t stR=0; stR<sR; stR++)
		{	
			double cc_real = mIn(stR, stC).real();
			double cc_imag = -mIn(stR, stC).imag();
			mIn(stR, stC) = COMPLEX_DOUBLE(cc_real, cc_imag);
		}
	return ;
}

//矩阵转置
template <class _Ty>
void T(CMatrix<_Ty>& mIn)
{
	size_t sR, sC;
	sR = mIn.GetRowNum();			//取原矩阵行数
	sC = mIn.GetColNum();			//取原矩阵列数
	
	CMatrix<_Ty> mTemp(sC, sR);		//生成一新阵，行数与列数与原阵互换
	
	for(size_t stC=0; stC<sC; stC++)
		for(size_t stR=0; stR<sR; stR++)
			mTemp(stC, stR) = mIn(stR, stC);	//对新阵赋值
		mIn = mTemp;
		return ;					//返回新的转置阵
}

//矩阵共轭转置
template <class _Ty>
void H(CMatrix<_Ty>& mIn)
{
	size_t sR, sC;
	sR = mIn.GetRowNum();	//取原矩阵行数
	sC = mIn.GetColNum();	//取原矩阵列数
	
	CMatrix<COMPLEX_DOUBLE> mTemp(sC, sR);
	
	for ( size_t i = 0; i < sR; i++ )
		for ( size_t j = 0; j < sC; j++)
			mTemp(i, j) = conjg(mIn(i,j));
		
		T(mTemp);
		mIn = mTemp;
		
		return ;
}

//全选主元高斯-约当(Gauss-Jordan)法求矩阵逆
template <class _Ty>
int Inv(CMatrix<_Ty >& rhs)
{
	size_t stRank = rhs.GetColNum();	// 矩阵阶数
	if(stRank != rhs.GetRowNum())
		return int(-1);					// rhs不是方阵
	
	valarray<size_t> is(stRank);		// 行交换信息
	valarray<size_t> js(stRank);		// 列交换信息
	
	CMatrix<_Ty> m(rhs);					// 生成一matrix对象
	
	for(size_t k = 0; k < stRank; k++)
	{	// 全选主元
		long double MaxValue(0);
		for(size_t i = k; i < stRank; i ++)
		{
			for(size_t j = k; j < stRank; j ++)
			{
				long double tmp(Abs(m(i, j)));	//求m(i,j)绝对值
				if(tmp > MaxValue)				//主元不在对角线上
				{
					MaxValue = tmp;
					is[k] = i;					//记录主元行数
					js[k] = j;					//记录主元列数
				}
			}
		}
		
		if(FloatEqual(MaxValue, 0)) 
			return int(0);						//主元为0，矩阵奇异
		
		if(is[k] != k)							//主元不在当前行
		{
			for(size_t j = 0; j < stRank; j ++)	//交换行元素
				swap(m(k, j), m(is[k], j));
		}
		
		if(js[k] != k)							//主元不在当前列
		{
			for(size_t i = 0; i < stRank; i ++)	//交换列元素
				swap(m(i, k), m(i, js[k]));
		}
		
		m(k, k) = 1.0 / m(k, k);				//主元倒数
		for(size_t j = 0; j < stRank; j ++)
			if(j != k)
				m(k, j) *= m(k, k);
			
			for(int i = 0; i < stRank; i ++)
				if(i != k)
					for(size_t j = 0; j < stRank; j ++)	
						if(j != k)
							m(i, j) = m(i, j) - m(i, k) * m(k, j);
						
						for(int i = 0; i < stRank; i ++)
							if(i != k)
								m(i, k) = -m(i, k) * m(k, k);
	}
	
	for(int r = stRank - 1; r >= 0; r--)
	{
		if(js[r] != r)
			for(size_t j = 0; j < stRank; j ++)
				swap(m(r, j), m(js[r], j));
			if(is[r] != r)
				for(size_t i = 0; i < stRank; i ++)
					swap(m(i, r), m(i, is[r]));
	}
	
	rhs = m;
	
	return int(1);
}


///////////////////////////////////////////////////////////
	//fft的函数
	// n:  输入的样本点数；k:满足n=(2)^k
	// l:  l＝0，表示求FFT变换；l＝1，表示求FFT反变化
	// il: il＝0，表示不求模和幅角；il＝1，表示要求模和幅角
	// 默认：l=0; il=1
	// matIn 表示输入的数据，返回时其实部存放模值，虚部存放幅角
	// matOut 表示FFT后输出的数据
template <class _Ty>
CMatrix<_Ty>  fft( CMatrix<_Ty>& matIn,
		  int l,int il)
{
//////////////////////////////////////////////////////////////////////////////////////////
// 把矩阵类matIn转化成数组
// n:  输入的样本点数；k:满足n=(2)^k
	int m1, n;
	int k;
	n = matIn.GetRowNum();
	m1= matIn.GetColNum();
	//k =(int)_logb(n);

	CMatrix<_Ty> matOut(n, m1);
	double *pr = new double[n*m1];
	double *pi = new double[n*m1];
	double *fr = new double[n*m1];
	double *fi = new double[n*m1];

	if( n == 1)
	{
		n = m1;
		k = static_cast<int>( (log10((double) n))/(log10(2.0)) );
		for ( int j = 0; j < n; j++)
		{
			pr[j] = matIn(0, j).real();
			pi[j] = matIn(0, j).imag();

		}

		//////////////////////////////////////////////////////////////////////////////////////////
		int it,m,is,i,j,nv,l0;
		double p,q,s,vr,vi,poddr,poddi;
		for (it=0;it<=n-1;it++)
		{
			m=it;is=0;
			for (i=0;i<=k-1;i++)
			{
				j=m/2;
				is=2*is+(m-2*j);
				m=j;
			}
			fr[it]=pr[is];fi[it]=pi[is];
		}
		pr[0]=1.0;pi[0]=0.0;
		p=6.283185306/(1.0*n);
		pr[1]=cos(p);pi[1]=-sin(p);
		if (l!=0) pi[1]=-pi[1];
		for (i=2;i<=n-1;i++)
		{
			p=pr[i-1]*pr[1];q=pi[i-1]*pi[1];
			s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
			pr[i]=p-q;pi[i]=s-p-q;
		}
		for (it=0;it<=n-2;it=it+2)
		{
			vr=fr[it];vi=fi[it];
			fr[it]=vr+fr[it+1];fi[it]=vi+fi[it+1];
			fr[it+1]=vr-fr[it+1];fi[it+1]=vi-fi[it+1];
		}
		m=n/2;nv=2;
		for (l0=k-2;l0>=0;l0--)
		{
			m=m/2;nv=2*nv;
			for (it=0;it<=(m-1)*nv;it=it+nv)
				for (j=0;j<=(nv/2)-1;j++)
				{
					p=pr[m*j]*fr[it+j+nv/2];
					q=pi[m*j]*fi[it+j+nv/2];
					s=pr[m*j]+pi[m*j];
					s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
					poddr=p-q;poddi=s-p-q;
					fr[it+j+nv/2]=fr[it+j]-poddr;
					fi[it+j+nv/2]=fi[it+j]-poddi;
					fr[it+j]=fr[it+j]+poddr;
					fi[it+j]=fi[it+j]+poddi;
				}
		}
		if(l!=0)
			for (i=0;i<=n-1;i++)
			{
				fr[i]=fr[i]/(1.0*n);
				fi[i]=fi[i]/(1.0*n);
			}
			if(il!=0)
				for (i=0;i<=n-1;i++)
				{
					pr[i]=(sqrt(fr[i]*fr[i]+fi[i]*fi[i]));
					pi[i]=atan2(fi[i],fr[i])*360.0/6.283185306;
				}
				/////////////////////////////////////////////////////////////////////////////////////////////////

				for ( int j = 0; j < n; j++ )
				{
					matOut(0, j) = COMPLEX_DOUBLE(fr[j ], fi[j]);
				}

				for ( int j = 0; j < n; j++ )
				{
					matIn(0, j) = COMPLEX_DOUBLE(pr[j], pi[j]);
				}
	}
	else
	{
		k = static_cast<int>( (log10((double) n))/(log10(2.0)) );
		for ( int nLength = 0; nLength < m1; nLength++)
		{

			for ( int j = 0; j < n; j++)
			{
				pr[j] = matIn(j, nLength).real();
				pi[j] = matIn(j, nLength).imag();

			}

			//////////////////////////////////////////////////////////////////////////////////////////
			int it,m,is,i,j,nv,l0;
			double p,q,s,vr,vi,poddr,poddi;
			for (it=0;it<=n-1;it++)
			{
				m=it;is=0;
				for (i=0;i<=k-1;i++)
				{
					j=m/2;
					is=2*is+(m-2*j);
					m=j;
				}
				fr[it]=pr[is];fi[it]=pi[is];
			}
			pr[0]=1.0;pi[0]=0.0;
			p=6.283185306/(1.0*n);
			pr[1]=cos(p);pi[1]=-sin(p);
			if (l!=0) pi[1]=-pi[1];
			for (i=2;i<=n-1;i++)
			{
				p=pr[i-1]*pr[1];q=pi[i-1]*pi[1];
				s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
				pr[i]=p-q;pi[i]=s-p-q;
			}
			for (it=0;it<=n-2;it=it+2)
			{
				vr=fr[it];vi=fi[it];
				fr[it]=vr+fr[it+1];fi[it]=vi+fi[it+1];
				fr[it+1]=vr-fr[it+1];fi[it+1]=vi-fi[it+1];
			}
			m=n/2;nv=2;
			for (l0=k-2;l0>=0;l0--)
			{
				m=m/2;nv=2*nv;
				for (it=0;it<=(m-1)*nv;it=it+nv)
					for (j=0;j<=(nv/2)-1;j++)
					{
						p=pr[m*j]*fr[it+j+nv/2];
						q=pi[m*j]*fi[it+j+nv/2];
						s=pr[m*j]+pi[m*j];
						s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
						poddr=p-q;poddi=s-p-q;
						fr[it+j+nv/2]=fr[it+j]-poddr;
						fi[it+j+nv/2]=fi[it+j]-poddi;
						fr[it+j]=fr[it+j]+poddr;
						fi[it+j]=fi[it+j]+poddi;
					}
			}
			if(l!=0)
				for (i=0;i<=n-1;i++)
				{
					fr[i]=fr[i]/(1.0*n);
					fi[i]=fi[i]/(1.0*n);
				}
				if(il!=0)
					for (i=0;i<=n-1;i++)
					{
						pr[i]=(sqrt(fr[i]*fr[i]+fi[i]*fi[i]));
						pi[i]=atan2(fi[i],fr[i])*360.0/6.283185306;
					}
					/////////////////////////////////////////////////////////////////////////////////////////////////

					for ( int j = 0; j < n; j++ )
					{
						matOut(j, nLength) = COMPLEX_DOUBLE(fr[j ], fi[j]);
					}

					for ( int j = 0; j < n; j++ )
					{
						matIn(j, nLength) = COMPLEX_DOUBLE(pr[j], pi[j]);
					}
		}
	}
	delete[] pr;
	delete[] pi;
	delete[] fr;
	delete[] fi;
				
	return matOut;
}

template <class _Ty>
CMatrix<_Ty>  fftshift( CMatrix<_Ty>& matIn)
{
	int m_Row, n_Col;
	m_Row = matIn.GetRowNum();
	n_Col = matIn.GetColNum();

	for( int mLength = 0; mLength < m_Row; mLength++)
	{
		
		for( int nLong = 0, nWidth = n_Col/2; nLong < n_Col/2; nLong++, nWidth++)
		{
			_Ty Temp;
			Temp = matIn(mLength, nLong);
			matIn(mLength, nLong) = matIn(mLength, nWidth);
			matIn(mLength, nWidth) = Temp;
		}
	}
	return matIn;
}
////////////////////////////////////////////////////////////////////////
    //SVD的函数
	//默认：MMAX=M,NMAX=N; P=0
	// 输入为matA（M*N），NU 左向量矩阵维数，NV右向量矩阵维数
	// matU左向量矩阵（NU*NU）,matV右向量矩阵(NV*NV)
template <class _Ty>
void SVD(CMatrix< _Ty> &matA, int MMAX, int NMAX, int M, int N, int P, int NU, int NV, double * S, CMatrix< _Ty> &matU, CMatrix< _Ty> &matV) //一般复矩阵的奇异值分解
{
	COMPLEX_DOUBLE* A = new COMPLEX_DOUBLE[(M+1)*(N+1)];
	COMPLEX_DOUBLE* U = new COMPLEX_DOUBLE[(M+1)*(N+1)];
	COMPLEX_DOUBLE* V = new COMPLEX_DOUBLE[(M+1)*(N+1)];
	
	for(int i = 0; i < (M+1)*(N+1); i++)
		A[i] = COMPLEX_DOUBLE(0.0E0,0.0E0);