mimo_channel.cpp

该函数可以实现任意行列数double型矩阵的张量乘积 用数组实现

#include "iostream.h"
#include "math.h"
#include <malloc.h>
#include "complex.h"

#define pi 3.1415

class  CMatrix             //矩阵类
{
public:

 
CMatrix();
CMatrix(const CMatrix &M);
CMatrix(int m, int n);
CMatrix(int m, int n, double * array);
void SetMatrix(double* array);
~CMatrix();
void SetMatrix(int m, int n);
int  GetM()const;
int  GetN()const;
CMatrix Transpose();
CMatrix & operator=(const CMatrix & H);
CMatrix operator+( const CMatrix& H);
CMatrix operator-(const CMatrix& H);
 CMatrix operator*( const CMatrix& H);
CMatrix operator*(double k);
double * & operator[](int i)const;
void display();


//double Surplus() ;
//CMatrix  MatrixOpp();

double matrix_det();
double matrix_alg_cofactor(int ai,int aj);
CMatrix matrix_inverse();
CMatrix cholesky ();



private:

  int m;     //行数
  int n;     //列数
  double ** p;     //二维数组指针

};



CMatrix::CMatrix()//无参构造函数
{
   m = 0; n  = 0;
   p = NULL;
}




 CMatrix::CMatrix(int m, int n)//有参构造函数,实现全零矩阵

{
    int i, j;
    p = new double*[m];
    for (i = 0; i < m; i++)
 {
      p[i] = new double[n];

       for (j = 0; j < n; j++)

       {
           p[i][j]=0;
       }
    }
    this->m = m; this->n = n;
}

 CMatrix::CMatrix(int m, int n, double* array)//在构造矩阵对象时候就直接将一个数组赋值到对象中
{
    int i, j;
    p = new double*[m];
    for (i = 0; i < m; i++)
    {
       p[i] = new double[n];
       for (j = 0; j < n; j++)
       {
          p[i][j] = array[i*n+j];//若传递int**array二维指针的话，元素可以用*(*(array+i)+j)表示
       }
   }
    this->m = m; this->n = n;
}

void CMatrix::SetMatrix(int m, int n)//设置矩阵的行列值，无参构造函数可以调用设置行列值

{
    int i, j;
    p = new double*[m];
    for (i = 0; i < m; i++)
 {
      p[i] = new double[n];

       for (j = 0; j < n; j++)

       {
           p[i][j]=0;
       }
    }
    this->m = m; this->n = n;
}
void CMatrix::SetMatrix(double* array) //在构造有参对象之后，将一个数组赋值进去
{ 
	int i, j; 
 	if (p) 
	{ 
		for (i = 0; i < this->m; i++) 
		{ 
			delete[] p[i]; 
		} 
		delete[] p; 
	} 


	 p = new double*[m];
    for (i = 0; i < m; i++)
      p[i] = new double[n];
	
    
 	for (i = 0; i < m; i++) 
	{ 
		for (j = 0; j < n; j++) 
		{ 
			p[i][j] = array[i*n+j]; 
		} 
	} 
} 

CMatrix::CMatrix(const CMatrix &M)//拷贝构造函数*****

{
	m=M.GetM();
    n=M.GetN();
    //赋值前先分配空间!
     p=new double *[m];
     for(int t=0; t<m; t++)
       p[t]=new double [n];
       
  for(int i=0;i<m;i++)
    for(int j=0;j<n;j++)
      p[i][j]=M.p[i][j];

}

 CMatrix::~CMatrix()//析构函数
 {
	int i;
    if (p)
    {
       for (i = 0; i < m; i++)
      {
           delete[] p[i];
      }
       delete[] p;
    }
 }


int  CMatrix::GetM()const
{
    return m;
}


int  CMatrix::GetN()const
{
    return n;
}



 double * & CMatrix::operator[](int i) const//重载运算符[]
{
    return p[i];
}

CMatrix CMatrix::Transpose() //实现矩阵的转置
{
    int i, j;
   
	CMatrix t(n,m);
 
      for (i = 0; i < m; i++)
       {
          for (j = 0; j < n; j++)
           {
              t[j][i] = p[i][j];
		  }
       }

     return t ;
 }







CMatrix & CMatrix::operator=(const CMatrix &H)//重载运算符=
{
    int i, j;
     if(this != &H)
{
     for ( i = 0 ;i < m;i++ )
      if(p[i])
      delete[] p[i];
      if(p)
   delete[] p;
 }

    m = H.GetM();
	n = H.GetN();
   p = new double*[m];
    for (i = 0; i < m; i++)
      p[i] = new double[n];

    for (i = 0; i < m; i++)
    	for (j = 0; j < n; j++)
       {
           p[i][j] = H[i][j];
       }
    

     return (*this);
}


CMatrix CMatrix::operator+(const CMatrix & H)//重在运算符+
 {
   int i, j;
   
   CMatrix t(m, n);
   
      for (i = 0; i < m; i++)
      {
        for (j = 0; j <n; j++)
		{
           t[i][j] = p[i][j] + H[i][j];
		}
	
      }
   
   return t;
}

CMatrix CMatrix::operator-(const CMatrix& H)//重在运算符-
{
   int i, j;
  CMatrix t(m, n);
   if (m == H.GetM() && n == H.GetN())    
   {
      for (i = 0; i < m; i++)
      {
        for (j = 0; j < n; j++)
        {
           t[i][j] = p[i][j] - H[i][j];
        }
      }
   }
   return t;
}


CMatrix CMatrix::operator*( const CMatrix& H) //重载运算符*
{
  int i, j, k;
  
    
  
	  CMatrix t(m, H.GetN());
      for (i = 0; i < m; i++)
      {
        for (j = 0; j < H.GetN(); j ++)
        {
           for (k = 0; k < n; k++)
           {
            
              t[i][j] =  t[i][j]+p[i][k] * H[k][j];
           }
        }
      } 
	  
   
 return t;
  

}

CMatrix CMatrix::operator*(double k)//重载运算符*，实现实数和矩阵相乘
{
   int i, j;
   CMatrix t(m, n);
   for (i = 0; i < m; i++)
   {
      for (j = 0; j < n; j++)
      {
        t[i][j] =p[i][j] * k  ;
      }
   }
  return t;
}

void CMatrix::display()
{
	for (int i=0;i<m;i++)
	{
		for(int j=0;j<n;j++)
	{
		cout<<p[i][j];
	}
   	cout<<endl;
	}
}

 




double CMatrix::matrix_det()//求矩阵行列式，可以作为工具函数使用
{
	if( m!=n) 
	{cout<<"此矩阵不是方阵,没有逆矩阵，故不存在行列式!"<<endl;return 0;}

  double det = 0;
    
  if (m == 1&&n==1)
 {
  det = p[0][0];
 }
 else if (m==2&&n == 2) 
 {
  det = p[0][0] * p[1][1]- p[0][1] * p[1][0];
 }
  else 
  	  
 {
 
 
   CMatrix c(this->m-1,this->n-1);
 // double * t = new double [(m-1)*(n-1)];
  for (int k = 1;k <= n;k++)
  { 
	 
    int temp=1;
//	int temp1=1;
   for (int i = 1;i <=this-> m;i++)
    for (int j = 1;j <= this->n;j++)
    {
     if (i != 1 && j != k) 
     {
		 
		
     //*(*(c.p+i-2)+temp-1)=*((*((this->p)+i-1)+j-1));
   *(*(c.p+(temp-1)/(m-1))+(temp-1)%(this->n-1))=*((*((this->p)+i-1)+j-1));//重点考虑如何实现规律
                                                                           //如果类中数组指针为一维时，比较好实现规律
                                                                           //(c.p+temp-1) temp++ 再对实际的二维数组做强制转换    
	//	t[temp-1] =*((*((this->p)+i-1)+j-1));

       temp++;
	  // temp1++;
     }
    }
//	c.SetMatrix(t);
//	delete []t;
   if (k%2 == 1)
   
    det += p[0][k-1] * c.matrix_det();
   
   else
   
    det -= p[0][k-1] * c.matrix_det();
   
 
 }
 
 }
  return det;
}


//矩阵的代数余子式
double CMatrix::matrix_alg_cofactor(int ai,int aj)//可以作为工具函数使用
{
 int i,j;
 double ac;
// double * t = new double [(m-1)*(n-1)];
 CMatrix c(m-1,n-1);
 //ai = (k - 1) / n + 1;
// aj = (k - 1) % n + 1;
// m = 1;
 int temp=1;
 //int  temp1=1;
 for (i = 1;i <= m;i++)
  for (j = 1;j <= n;j++)
  { 
   if (i != ai && j != aj) 
   {
    // t[temp-1] = *((*((p+i-1))+j-1)); 
	  //  *(*(c.p+temp-1)+temp-1)=*((*((p+i-1))+j-1));//重点考虑如何实现规律
	   *(*(c.p+(temp-1)/(m-1))+(temp-1)%(this->n-1))=*((*((p+i-1))+j-1));
        temp++;
	 //  temp1++;
   // m++;
   }
  }

  //c.SetMatrix(t);
 // delete []t;

 ac = c.matrix_det();
 

 if ((ai + aj) %2 == 0)
 {
  return ac;
 }
 else
 {
  return 0 - ac;
 }
}



CMatrix CMatrix::matrix_inverse()//矩阵的逆函数
{

 double det = this->matrix_det();
 if (det==0)
 {
 CMatrix a(m,n);
 cout<<"不可逆"<<endl;
 return a;
 }
  int i,j;
 // double * temp = new double[n*n];
  CMatrix c(m,n);
  CMatrix d;
 // for (i = 1;i <= n * n;i++)
  for( i=1;i<=m;i++)
	  for( j=1;j<=n;j++)
		  c[i-1][j-1]=(this->matrix_alg_cofactor(i,j))/det;
           
	d= c.Transpose();
	  return (d);
 }



//int cholesky (double * G,int n,double * L,double * D)

CMatrix CMatrix::cholesky ()// 对称正定矩阵的Cholesky分解
{
 
if(m!=n)//可以先判断是不是方阵

{
	cout<<"不是方阵，不可cholesky分解"<<endl;
	return *this;
}
 int i,j,r;
 double temp;
 //double * C = new double [n * n];
 CMatrix C(n,n);
 CMatrix L(n,n);
 CMatrix D(n,n);
 for (i = 1;i <= n * n;i++)
 {
  if ((i-1)%n == (i-1)/n) //注意下标的转换关系
  // L[i-1] = 1;
  L[(i-1)/(n)][(i-1)%(n)]=1;

  else
   L[(i-1)/(n)][(i-1)%(n)]= 0; 
  D[(i-1)/(n)][(i-1)%(n)] = 0;

 } 
 
 for(j = 1;j <= n;j++)
 {
  temp = 0;
  for(r=1;r<=j-1;r++)
  // temp += C[(j-1)*n+r-1] * C[(j-1)*n+r-1] / D[r-1];
  temp += C[j-1][r-1] * C[j-1][r-1] / D[(r-1)/(n)][(r-1)%(n)];
  //D[j-1] = G[(j-1)*n+j-1] - temp;
  D[(j-1)/(n)][(j-1)%(n)] = p[j-1][j-1] - temp;
  if ( D[(j-1)/(n)][(j-1)%(n)] <= 0) 
  {
    cout<<"该方阵不可cholesky分解"<<endl;
	return *this;break;
   
  }
  if (j == n)
  {
   for (r = 2;r <= n;r++)
   {
    //D[(r-1)*n+r-1] = D[r-1];
    D[r-1][r-1] = D[(r-1)/(n)][(r-1)%(n)];
    //D[r-1] = 0;
    D[(r-1)/(n)][(r-1)%(n)]=0;
   }
   
   return L;
  }
  for (i = j + 1;i <= n;i++)
  {
   temp = 0;
   for(r=1;r<=j-1;r++)
    //temp += C[(j-1)*n+r-1] * C[(i-1)*n+r-1] / D[r-1];
	temp += C[j-1][r-1] * C[i-1][r-1] /  D[(r-1)/(n)][(r-1)%(n)];
  // C[(i-1)*n+j-1] = G[(i-1)*n+j-1] - temp;
   C[i-1][j-1] = p[i-1][j-1] - temp;
   L[i-1][j-1] = C[i-1][j-1] / D[(j-1)/(n)][(j-1)%(n)];
  }
 }
 
 cout<<"出现错误"<<endl;
 return *this;
}









/*****************************************************************************************

ostream& operator<<(ostream&output,matrix&c)//重载输出函数，只能作为友元函数
{
	for (int i=0;i<row;i++)
	{
		for(int j=0;j<colum;j++)
	{
		//	output<<c.c[i][j]<<" ";
			c.c[i][j].display();
	}
	//	output<<endl;
	}
	return output;
 
}

istream& operator>>(istream &intput,matrix &c)//重载输入函数，只能作为友元函数
{

	cout<<"please input the matrix:"<<endl;
	   for (int i=0;i<row;i++)
			for(int j=0;j<colum;j++)
			//	intput>>c.c[i][j];
			c.c[i][j].setinitial(1,1);
			return intput;
}

*********************************************************/


/*******************************************
          channel function
  *****************************************/

/*
CMatrix MIMO_Channel(int Nr,int Nt,int t)
{
	int s=35;
	double fd=5.56;
	double h1,h2;
	double sita[35];
	CMatrix H(Nr,Nt);
	CMatrix HH(Nr,Nt);
	CMatrix CorrR;
	CMatrix CorrT;
	CMatrix CorrRT;
	for(int i=0;i<Nt;i++)
	 for(int j=0;i<Nr;j++)
	 {
		h1=0;
		h2=0;
		for(int k=1;k<=s-1;k++)
		{
		        sita[k]=2*pi*rand();
			h1+=sqrt(2)/sqrt(s-1/2)*sin(pi*k/(s-1))*cos(2*pi*fd*cos(pi*k/(2*s-1))*t+sita[k]);//产生独立分布高斯变量//添加时延
			h2+=sqrt(2)/sqrt(s-1/2)*cos(pi*k/(s-1))*cos(2*pi*fd*cos(pi*k/(2*s-1))*t+sita[k]);//产生独立分布高斯变量//添加时延

		}
	sita[35]=rand();
	h1+=1/(sqrt(2)*sqrt(s-1/2))*cos(2*pi*fd*t+sita[35]);//添加时延
        h2+=1/(sqrt(2)*sqrt(s-1/2))*cos(2*pi*fd*t+sita[35]);//添加时延
	complex c(h1,h2);
	H[j][i]=c;


	 }

	 CorrorT=MIM0_Corror(30,0,0.5,Nt);
	 CorrorR=MIMO_Corror(5,0,5,Nr);
	 CorrRT=Kron(CorrR,CorrT);
	 hr=Transpose(chol(CorrRT));//cholesky平方根求方阵的分解
//	 double f[16]={1,0.3409,0.5442,0.1855,0.3409,1,0.1855,0.5442,0.5442,0.1855,1,0.3409,0.1855,0.5442,0.3409,1};
//	 CMatrix hr(4,4,f);
	 H=hr*H;

	 for(int p=1;p<=Nr;p++)
		 for(int q=1;q<=Nt;q++)
		 {
			 HH[p][q]=H[p][q-1];
		 }

  return HH;

}

CMatrix MIMO_Corr(int anglespread,double angle,float d,int m)//相关矩阵是个方阵
{
	int L=1000;
	int anglespread1=720;
	int c=0;

	CMatrix p(1,L);
	CMatrix fai(1,L);
	CMatrix fai1(1,L);
	CMatrix FAI(1,L);
	CMatrix matrix1(m,1);
	CMatrix matrix2(1,m);
	CMatrix correlation1(m,m);
	CMatrix correlation2(m,m);
	CMatrix correlation(m,m);

	for (int i=0;i<L;i++)
	{
		fai1[0][i]=angle-anglespread1+2*anglespread1*i/L;
		fai[0][i]=2*pi*(angle-anglespread1+2*anglespread1*i/L)/360;
		FAI[0][i]=d*sin(fai[1][i]);
	}

	for(int i=0;i<L;i++)
		p[0][i]=1/(anglespread*sqrt(2))*exp(-sqrt(2)*abs(fai1[0][i]-angle)/anglespread)*2*anglespread1/L;

	for(int i=0;i<L;i++)
		c=p[0][i]+c;

	for(int j=0;j<L;j++)
	{
		for(int k=0;k<m;j++)
	 {matrix1[k][0]=exp(i*FAI[0][j]*2*pi*(k));}
		matrix2=matrix1.Transpose();
		correlation1=matrix1*matrix2*p[0][j];
		correlation2=correlation1+correlation2;

	}

   for(int i=0;i<m;i++)
	   for(int j=0;j<m;j++)
		   correlation[i][j]=abs(correlation2[i][j])/c;

   
   return correlation;
   
}


***********************/






void main()
{


double r=0,q=0,z=9,d=7,e=6,g=6,t=4,i=3,j=2,k=7,l=6,y=1;
double f[16]={1,0.3409,0.5442,0.1855,0.3409,1,0.1855,0.5442,0.5442,0.1855,1,0.3409,0.1855,0.5442,0.3409,1};
double h[9]={1,2,3,2,1,2,3,2,1};
double u[4]={0,8,6,0};
CMatrix a(4,4,f);
CMatrix b(3,3,h);
CMatrix c(2,2,u);
CMatrix v;
//a.display();
//b.display();
//c=a*b;
//c.display();


//v=b.matrix_inverse();

v=a.cholesky();
v.display();
//v=a.cholesky();











}