generate hermitian positive definite matrix

[naufal]

Hi,

I would like to generate a hermitian positive definite matrix Z based on random rayleigh fading channel matrix H. The rayleigh fading channel with i.i.d, zero-mean, and unit-variance complex Gaussian random variables. This matrix is used for MIMO system. I successfully generated the Rayleigh fading channel matrix H and matrix Z in MATLAB as below:

K=4;
M=16;
H=(1/sqrt(2))*(randn(K,M) + 1i*(randn(K,M));
Z=H*H';

Unfortunately, I don't have any idea how to generate this matrix in c++.
Hope to hear from you guys.

Thank you.

[dutch]

First you need to choose a library. Both the C and C++ libraries listed here will work with C++.
https://en.wikipedia.org/wiki/List_of_numerical_libraries#C

GSL (C) is very popular: https://www.gnu.org/software/gsl/

Eigen (C++) is supposed to be pretty cool: http://eigen.tuxfamily.org/index.php?title=Main_Page

[lastchance]

#include <iostream>
#include <iomanip>
#include <vector>
#include <complex>
#include <random>
#include <ctime>
using namespace std;

using vec    = vector< complex<double> >;
using matrix = vector< vec >;


mt19937 gen( time( 0 ) );
normal_distribution<double> dist( 0.0, 1.0 );


//======================================================================


void print( const matrix &M, int width )
{
   for ( auto &row : M )
   {
      for ( auto z : row ) cout << setw( width ) << z << "  ";
      cout << '\n';
   }
}


//======================================================================


matrix genNormal( int rows, int cols )
{
   const double sq = sqrt( 2.0 );
   matrix result( rows, vec( cols ) );

   for ( int i = 0; i < rows; i++ )
   {
      for ( int j = 0; j < cols; j++ ) result[i][j] = complex<double>( dist( gen ), dist( gen ) ) / sq;
   }
   return result;
}


//======================================================================


matrix Hermitian( const matrix &A )
{
   int N = A.size();
   int nk = A[0].size();
   matrix result( N, vec( N ) );
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ )
      {
         result[i][j] = 0.0;
         for ( int k = 0; k < nk; k++ ) result[i][j] += A[i][k] * conj( A[j][k] );
      }
   }
   return result;
}


//======================================================================


int main()
{
   const int K = 4, M = 16;

   matrix H = genNormal( K, M );
   matrix Z = Hermitian( H );
   print( Z, 20 );
}

         (16.6123,0)    (2.96544,-4.94904)   (5.62544,0.0848946)  (-0.0605223,1.80707)  
   (2.96544,4.94904)           (16.2707,0)    (-2.08518,7.43584)     (-2.8041,1.47581)  
(5.62544,-0.0848946)   (-2.08518,-7.43584)           (19.6652,0)      (2.7578,0.95839)  
(-0.0605223,-1.80707)    (-2.8041,-1.47581)     (2.7578,-0.95839)           (13.3219,0)  

[naufal]

Thank you.