function [mod_symbols,table, P] = tx_modulate(bits_in, modulation)

%**************************************************************************
%This program modulates the input binary data.The inputs are:
% bits_in -> the binary input bits
%modulation -> we choose one of 'BPSK','QPSK',8PSK',16QAM'
%The outputs are:
%mod_symbols -> modulated output
%table -> the complete set of symbols in a chosen constellation. This is
%required by the demodulator.
%P -> the number of points in a chosen constellation. This is required by
%the demodulator
%**************************************************************************

full_len = length(bits_in);

% BPSK modulation
if ~isempty(findstr(modulation, 'BPSK '))
   % Angle [0 -pi] corresponds to 
   % Gray code vector [1 0], respectively.
   %table=exp(j*[0 -pi])./sqrt(0.5);  % generates BPSK symbols
   table=exp(j*[0 -pi])./sqrt(0.5);  % generates BPSK symbols
   table=table([1 0]+1); % Gray code mapping pattern for BPSK symbols
   inp=bits_in;
   mod_symbols=table(inp+1);  % maps transmitted bits into BPSK symbols
   P=2;% 2 constellation points
   % QPSK modulation
elseif ~isempty(findstr(modulation, 'QPSK '))
   % Angle [3/4*pi -3/4*pi -1/4*pi 1/4*pi] corresponds to 
   % Gray code vector [00 01 11 10], respectively.
   table=exp(j*[3/4*pi -3/4*pi -1/4*pi 1/4*pi]);  % generates QPSK symbols
   table=table([0 1 3 2]+1); % Gray code mapping pattern for QPSK symbols
   inp=reshape(bits_in,2,full_len/2);
   mod_symbols=table([2 1]*inp+1);  % maps transmitted bits into QPSK symbols
   P=4;% 4 constellation points
elseif  ~isempty(findstr(modulation, '8QAM '))
    % generate constellations
       % Angle [0 pi/4 pi/2 3*pi/4 pi 5*pi/4 3*pi/2 7*pi/4] corresponds to 
   % Gray code vector [000 001 011 010 110 111 101 100], respectively.
   %table=exp(j*[0 pi/4 pi/2 3*pi/4 pi 5*pi/4 3*pi/2 7*pi/4]);  % generates 8PSK symbols
   m=1;
   for k=-3:2:3
      for l=-1:2:1
         table(m) = (k+j*l)/sqrt(3); % power normalization
         m=m+1;
      end;
   end;
   table=table([1 0 3 2 7 6 5 4]+1); % Gray code mapping pattern for 8QAM symbols
   inp=reshape(bits_in,3,full_len/3);
   mod_symbols=table([4 2 1]*inp+1);  % maps transmitted bits into 8QAM symbols
   P=8;%8 constellation points
   % 16-QAM modulation
elseif ~isempty(findstr(modulation, '16QAM'))
   % generates 16QAM symbols
%       m=1;
%    for k=-3:2:3
%       for l=-3:2:3
%          table(m) = (k+j*l)/sqrt(18); % power normalization
%          m=m+1;
%       end;
%    end;
%    table=table([3 2 0 1 7 6 4 5 15 14 12 13 11 10 8 9]+1); % Gray code mapping pattern for 16QAM symbols
%    inp=reshape(bits_in,4,full_len/4);
%    mod_symbols=table([8 4 2 1]*inp+1);  % maps transmitted bits into 16QAM symbols
%    P=16; %16 constellation points
   inp=reshape(bits_in,4,full_len/4);
   mod_symbols=qammod([8 4 2 1]*inp,16,[],'gray')./sqrt(5);  % maps transmitted bits into 16QAM symbols
   table=qammod(0:15,16,[],'gray')./sqrt(5);
   P=16; %16 constellation points
elseif ~isempty(findstr(modulation, '32QAM'))
   % generates 32QAM symbols
%       m=1;
%    for k=-5:2:5
%       for l=-5:2:5
%          table(m) = (k+j*l)/sqrt(34); % power normalization
%          m=m+1;
%       end;
%    end;
%    table=table([3 2 0 1 7 6 4 5 15 14 12 13 11 10 8 9]+1); % Gray code mapping pattern for 32QAM symbols
   inp=reshape(bits_in,5,full_len/5);
   mod_symbols=qammod([16 8 4 2 1]*inp,32,[],'gray')./sqrt(10);  % maps transmitted bits into 32QAM symbols
   table=qammod(0:31,32,[],'gray')./sqrt(10);
   P=32; %32 constellation points
elseif ~isempty(findstr(modulation, '64QAM'))
   % generates 64QAM symbols
%       m=1;
%    for k=-3:2:3
%       for l=-3:2:3
%          table(m) = (k+j*l)/sqrt(18); % power normalization
%          m=m+1;
%       end;
%    end;
%    table=table([3 2 0 1 7 6 4 5 15 14 12 13 11 10 8 9]+1); % Gray code mapping pattern for 64QAM symbols
   inp=reshape(bits_in,6,full_len/6);
   mod_symbols=qammod([32 16 8 4 2 1]*inp,64,[],'gray')./sqrt(21);  % maps transmitted bits into 64QAM symbols
   table=qammod(0:63,64,[],'gray')./sqrt(21);
   P=64; %64 constellation points


else
   error('Unimplemented modulation');
end