function Z = rndgg(n,m,s,b,mu,cq)
% RNDGG  random sample generation from unit dispersion complex generalized Gaussian (GG) distribution  
%                 
% INPUT 
%    m, n  size of the sample n    
%    s     exponent parameter (positive real number, s >0)
%    b     scale parameter (either b > 0, or, b = 'unitvar'). In the latter case, scale b is 
%          assined the value that yields unit variance (var(z)=1) GG distribution  
%    mu    mean (complex number)
%    cq    circularity quotient (complex number) with property  |cq|<= 1. Since dispersion = 1, cq  
%          also equals the pseudo-dispersion parameter of the distribution
% OUTPUT 
%    Z     random sample of size  m x n  from  GG_{s,b}(mu,dispersion=1,pseudo-dispersion=cq) 
% 
% (c) Esa Ollila, contact: esollila AT wooster.hut.fi
%
% If you use this program, please cite 
%    E. Ollila, J. Eriksson, V. Koivunen, "Complex univariate elliptically symmetric distributions", 
%       submitted to IEEE Transactions on Signal Processing, Dec. 2009. 
%-------------------------------------------------------------------------------

if ~license('test','statistics_toolbox')
   error('Program requires statistics toolbox: depends on GAMRND function')
end

% CHECK OUT THE VALIDITY OF THE INPUT PARAMETERS 

if ( ~isreal(s) ||  s <= 0 ) 
     error('Input scale parameter should be real-valued positive');
end  

if (b == 'unitvar')
   b = (gamma(1/s)/gamma(2/s))^s;
end 


if ( ~isreal(b)  || b <= 0) 
     error('Input exponent parameter should be real-valued positive');
end 

if (abs(cq)>1)  
     error('Input circularity quotient should have modulus <= 1')
end 


% STEP 1
cc  = abs(cq);    % cc = |cq| is the circularity coefficient
ang = angle(cq);  % ang = Arg(om) is the circularity angle

% STEP 2
R  =  sqrt(gamrnd(1/s,b,m,n).^(1/s));
Z0 = R.*exp(j*unifrnd(-pi,pi,m,n)); % random sample from circular GG_{s,b}(mu=0,dispersion=1)

% STEP 3
a = sqrt(1+cc); 
b = sqrt(1-cc); 
Zp = (1/2)*(a+b)*Z0 + (1/2)*(a-b)*conj(Z0);

% STEP 4
Z = mu + exp(j*ang/2)*Zp; 
