How can I model a Rayleigh channel gains

I want to generate a 4*4 [a11,a12 ; a21,a22]

path gains which are modeled by independent complex Gaussian random variables. The real and imaginary parts of the path gains at each time slot are i.i.d Gaussian random variables which have a zero mean and a variance equal to 0.5. The envelope of the channel path gains|aij| has a Rayleigh distribution. Also,|aij|^2 is a chi-square random variable with two degrees of freedom and the average channel energy is E[|aij|^2].