function  [x]=qam_gen(M,fl,t,d)
% [x]=qam_gen(M,fl,t,d);
%  example:
% [x]=qam_gen(16,1000,0.5,1);
% 'M'--qam modulation type 16/32/64
% 'fl'-- signal rate
% 't'-- time last
% 'd'-- min distance between symbols
error(nargchk(1,4,nargin));
if M<=0
   error('Wrong parameter, the value of M must be greater than 0')
end
if mod(M,4)~=0
   error('Wrong parameter, the value of M must be 4,16,32,64')
end
if nargin==1
   N=1000;d=1;
elseif nargin==3
   d=1;
end
N=fix(fl*t);
%generation of the data source follows
% 16QAM rand('state',308); 1198.3Hz   for del_fc=1200Hz
% 32QAM rand('state',118); 1200.4Hz
% 64QAM rand('state',171); 1197.7Hz
for i=1:N
   temp=rand;  % a uniform R.V. between 0 and 1
   dsource(i)=1+floor(M*temp);  % a number between 1 and 16, uniform
end;
% mapping to the signal constellation follow
map16 = [ -3*d 3*d; -d 3*d; d 3*d; 3*d 3*d;
          -3*d   d;   -d d; d   d; 3*d   d;
           -3*d -d;  -d -d;  d -d; 3*d  -d;
         -3*d -3*d;-d -3*d;d -3*d; 3*d -3*d];
map32=[      -3*d  5*d;-d  5*d; d  5*d;3*d  5*d;
   -5*d 3*d; -3*d  3*d;-d  3*d; d  3*d;3*d  3*d; 5*d 3*d;
   -5*d 3*d; -3*d    d;-d    d; d    d;3*d    d; 5*d d;
   -5*d  -d; -3*d   -d;-d   -d; d   -d;3*d   -d; 5*d -d;
   -5*d -3*d;-3*d -3*d;-d -3*d; d -3*d;3*d -3*d; 5*d -3*d;
   -3*d -5*d;-d -5*d; d -5*d;3*d -5*d];
map64=[-7*d  7*d;-5*d  7*d;-3*d  7*d;-d  7*d;d  7*d;3*d  7*d;5*d  7*d;7*d 7*d;
       -7*d  5*d;-5*d  5*d;-3*d  5*d;-d  5*d;d  5*d;3*d  5*d;5*d  5*d;7*d 5*d;
       -7*d  3*d;-5*d  3*d;-3*d  3*d;-d  3*d;d  3*d;3*d  3*d;5*d  3*d;7*d 3*d;
       -7*d    d;-5*d    d;-3*d    d;-d    d;d    d;3*d    d;5*d    d;7*d   d;
       -7*d   -d;-5*d   -d;-3*d   -d;-d   -d;d   -d;3*d   -d;5*d   -d;7*d  -d;
       -7*d -3*d;-5*d -3*d;-3*d -3*d;-d -3*d;d -3*d;3*d -3*d;5*d -3*d;7*d -3*d
       -7*d -5*d;-5*d -5*d;-3*d -5*d;-d -5*d;d -5*d;3*d -5*d;5*d -5*d;7*d -5*d;
       -7*d -7*d;-5*d -7*d;-3*d -7*d;-d -7*d;d -7*d;3*d -7*d;5*d -7*d;7*d -7*d];

 for i=1:N
    if M==16
       qam_sig(i,:)=map16(dsource(i),:);
    elseif M==32
       qam_sig(i,:)=map32(dsource(i),:);
    elseif M==64
       qam_sig(i,:)=map64(dsource(i),:);
    end
 end;
 x1=qam_sig(:,1).';x2=qam_sig(:,2).';
 x=x1+j*x2;
 return