## Copyright (C) 2003 David Bateman
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File}  {} qaskenco (@var{m})
## @deftypefnx {Function File} {} qaskenco (@var{msg},@var{m})
## @deftypefnx {Function File} {@var{y} = } qaskenco (@var{...})
## @deftypefnx {Function File} {[@var{inphase}, @var{quadr}] =} qaskenco (@var{...})
##
## Map a digital signal using a square QASK constellation. The argument
## @var{m} must be a positive integer power of 2. With two input arguments
## the variable @var{msg} represents the message to be encoded. The values
## of @var{msg} must be between 0 and @code{@var{m}-1}. In all cases
## @code{qaskenco(@var{M})} is equivalent to @code{qaskenco(1:@var{m},@var{m})}
##
## Three types of outputs can be created depending on the number of output 
## arguments. That is
##
## @table @asis
## @item No output arguments
## In this case @dfn{qaskenco} plots the constellation. Only the
## points in @var{msg} are plotted, which in the case of a single input
## argument is all constellation points.
## @item A single output argument
## The returned variable is a complex variable representing the in-phase 
## and quadrature components of the mapped  message @var{msg}. With, a 
## single input argument this effectively gives the mapping from symbols
## to constellation points
## @item Two output arguments
## This is the same as two ouput arguments, expect that the in-phase
## and quadrature components are returned explicitly. That is
##
## @example
## octave> c = qaskenco(msg, m);
## octave> [a, b] = qaskenco(msg, m);
## octave> all(c == a + 1i*b)
## ans = 1
## @end example
## @end table
##
## If @code{sqrt(@var{m})} is an integer, then @dfn{qaskenco} uses a Gray
## mapping. Otherwise, an attempt is made to create a nearly square mapping 
## with a minimum Hamming distance between adjacent constellation points.
## @end deftypefn
## @seealso{qaskdeco}

## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com>
##     * modified for new gnuplot interface (octave > 2.9.0)

function [a, b] = qaskenco(msg, M)

  if (nargin == 1)
    M = msg;
  elseif (nargin == 2)
    if ((min(msg(:)) < 0) || (max(msg(:)) > M-1))
      error ("qaskenco: message invalid");
    endif
  else
    error ("qaskenco: incorrect number of arguments");
  endif

  if (!isscalar(M) || (M != ceil(M)) || (M < 2))
    error ("qaskenco: order of modulation must be a positive integer greater than 2");
  endif

  if (log2(M) != ceil(log2(M)))
    error ("qaskenco: the order must be a power of two");
  endif
  
  if (M == 2)
    inphase = [-1,  1];
    quadr =   [ 0,  0];
  elseif (M == 4)
    inphase = [-1, -1,  1,  1];
    quadr =   [-1,  1, -1,  1];
  elseif (M == 8)
    inphase = [-1, -1,  1,  1, -3, -3,  3,  3];
    quadr =   [-1,  1, -1,  1, -1,  1, -1,  1];
  else
    NC =2^floor(log2(sqrt(M)));
    MM = NC * NC;
    Gray = [0, 1];
    for i=2:ceil(log2(NC))
      Gray = [Gray 2^(i-1) + fliplr(Gray)];
    end
    Gray = fliplr(de2bi(shift(Gray,length(Gray)/2 - 1)));
    Gray2 = zeros(MM,log2(MM));
    Gray2(:,1:2:log2(MM)) = repmat(Gray,NC,1);
    for i=1:NC
      Gray2(i:NC:MM,2:2:log2(MM)) = Gray;
    end
    layout = reshape(bi2de(fliplr(Gray2)),NC,NC);

    if (MM != M)
      ## Not sure this is the best that can be done for these mappings. If 
      ## anyone wants to improve this, go ahead, but do it in qaskdeco too.
      OFF = sqrt(M/32);
      NR = NC + 2*OFF;
      layout2 = NaN * ones(NR);
      layout2(1+OFF:OFF+NC,1+OFF:OFF+NC) = layout;
      
      layout2(1:OFF,1+OFF:OFF+NC) = MM + layout(OFF:-1:1,:);
      layout2(NR-OFF+1:NR,1+OFF:OFF+NC) = MM + layout(NC:-1:NC-OFF+1,:);

      layout2(1+2*OFF:NC,1:OFF) = MM + layout(OFF+1:NC-OFF,OFF:-1:1);
      layout2(1+2*OFF:NC,NR-OFF+1:NR) = MM + ...
	  layout(OFF+1:NC-OFF,NC:-1:NC-OFF+1);

      layout2(1+OFF:2*OFF,1:OFF) = MM + ...
	  layout(NC/2:-1:NC/2-OFF+1,NC/2:-1:OFF+1);
      layout2(NC+1:OFF+NC,1:OFF) = MM + ...
	  layout(NC-OFF:-1:NC/2+1,NC/2:-1:OFF+1);

      layout2(1+OFF:2*OFF,NR-OFF+1:NR) = MM + ...
	  layout(NC/2:-1:NC/2-OFF+1,NC-OFF:-1:NC/2+1);
      layout2(NC+1:OFF+NC,NR-OFF+1:NR) = MM + ...
	  layout(NC-OFF:-1:NC/2+1,NC-OFF:-1:NC/2+1);
      NC = NR;
      layout = layout2;
    endif

    inphase = repmat([0:NC-1]*2 - NC + 1,1,NC);
    for i=1:NC
      quadr(i:NC:NC*NC) = [0:NC-1]*2 - NC + 1;
    end
    [dummy, indx] = sort(layout(:));
    indx = indx(1:M);			## Get rid of remaining NaN's
    inphase = inphase(indx);
    quadr = quadr(indx);
  endif

  if (nargin == 2)
    inphase = inphase(msg+1);
    quadr = quadr(msg+1);
    ## Fix up indexing if using column vector
    if (size(msg,2) == 1)
      inphase = inphase';
      quadr = quadr';
    endif
  endif

  if (nargout == 0)
    inphase = inphase(:);
    quadr = quadr(:);
    plot (inphase, quadr, "r+");
    title("QASK Constellation");
    xlabel("In-phase");
    ylabel("Quadrature");
    axis([min(inphase)-1, max(inphase)+1, min(quadr)-1, max(quadr)+1]);
    xd = 0.02 * max(inphase);
    if (nargin == 2)
      msg = msg(:);
      for i=1:length(inphase)
	text(inphase(i)+xd,quadr(i),num2str(msg(i)));
      end
    else
      for i=1:length(inphase)
	text(inphase(i)+xd,quadr(i),num2str(i-1));
      end
    endif
  elseif (nargout == 1)
    a = inphase + 1i * quadr;
  else
    a = inphase;
    b = quadr;
  endif

endfunction