octave-communications-common 1.2.1-2 / usr / share / octave / packages / communications-1.2.1 / qaskdeco.m

## 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} {@var{msg} =} qaskdeco (@var{c}, @var{m})
## @deftypefnx {Function File} {@var{msg} =} qaskdeco (@var{inphase}, @var{quadr}, @var{m})
## @deftypefnx {Function File} {@var{msg} =} qaskdeco (@dots{}, @var{mnmx})
##
## Demaps an analog signal using a square QASK constellation. The input signal
## maybe either a complex variable @var{c}, or as two real variables
## @var{inphase} and @var{quadr} representing the in-phase and quadrature
## components of the signal.
##
## The argument @var{m} must be a positive integer power of 2. By default the
## same constellation as created in @code{qaskenco} is used by @code{qaskdeco}.
## If is possible to change the values of the minimum and maximum of the
## in-phase and quadrature components of the constellation to account for
## linear changes in the signal values in the received signal. The variable
## @var{mnmx} is a 2-by-2 matrix of the following form
##
## @multitable @columnfractions 0.125 0.05 0.25 0.05 0.25 0.05
## @item @tab | @tab min in-phase   @tab , @tab max in-phase   @tab |
## @item @tab | @tab min quadrature @tab , @tab max quadrature @tab |
## @end multitable
##
## If @code{sqrt (@var{m})} is an integer, then @code{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.
## @seealso{qaskenco}
## @end deftypefn

function a = qaskdeco (varargin)

  have_mnmx = 0;
  if (nargin == 2)
    c = varargin{1};
    inphase = real (c);
    quadr = imag (c);
    M = varargin{2};
  elseif (nargin == 3)
    if (all (size (varargin{3}) == [2 2]))
      c = varargin{1};
      inphase = real (c);
      quadr = imag (c);
      M = varargin{2};
      mnmx = varargin{3};
      have_mnmx = 1;
    else
      inphase = varargin{1};
      quadr = varargin{2};
      M = varargin{3};
    endif
  elseif (nargin == 4)
    inphase = varargin{1};
    quadr = varargin{2};
    M = varargin{3};
    mnmx = varargin{4};
    have_mnmx = 1;
  else
    print_usage ();
  endif

  if (iscomplex (inphase) || iscomplex (quadr))
    error ("qaskdeco: INPHASE and QUADR must be real");
  endif

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

  if (log2 (M) != ceil (log2 (M)))
    error ("qaskdeco: M must be a power of 2");
  endif

  if (have_mnmx)
    if (any (size (mnmx) != [2 2]))
      error ("qaskdeco: MNMX must be a 2-by-2 matrix of min and max values");
    endif
  else
    if ((M == 2) || (M == 4))
      mnmx = [-1, 1; -1, 1];
    elseif (M == 8)
      mnmx = [-3, 3; -1, 1];
    elseif (sqrt (M) == fix (sqrt (M)))
      NC = 2^floor (log2 (sqrt (M)));
      mnmx = [-NC+1, NC-1; -NC+1, NC-1];
    else
      NC = 2^floor (log2 (sqrt (M))) + 2*sqrt (M/32);
      mnmx = [-NC+1, NC-1; -NC+1, NC-1];
    endif
  endif

  if (M == 2)
    layout = [0, 1]';
  elseif (M == 4)
    layout = [0, 1; 2, 3];
  elseif (M == 8)
    layout = [4, 5; 0, 1; 2, 3; 6, 7];
  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)];
    endfor
    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;
    endfor
    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 qaskenco 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);
      layout = layout2;
    endif
  endif

  ix = 1 + (inphase - mnmx(1,1)) / (mnmx(1,2)-mnmx(1,1)) * (size (layout, 1) - 1);
  qx = 1 + (quadr - mnmx(2,1)) / (mnmx(2,2)-mnmx(2,1)) * (size (layout, 2) - 1);

  try    wfi = warning ("off", "Octave:fortran-indexing");
  catch  wfi = 0;
  end_try_catch

  unwind_protect
    a = layout(size (layout, 1) * (max (min (round (qx), size (layout, 2)), 1) - 1) + ...
               max (min (round (ix), size (layout, 1)), 1));
    ## FIXME: Why is this necessary??
    if ((M == 2) && (size (inphase, 1) == 1))
      a = a';
    endif

    if (any (isnan (a(:))))
      ## We have a non-square constellation, with some invalid points.
      ## Map to nearest valid constellation points...
      indx = find (isnan (a(:)));
      ix = ix(indx);
      qx = qx(indx);
      ang = atan2 (quadr(indx), inphase(indx));

      qx(find (ang > 3*pi/4)) = NR-OFF;
      ix(find ((ang <= 3*pi/4) & (ang > pi/2))) = OFF+1;
      ix(find ((ang <= pi/2) & (ang > pi/4))) = NR - OFF;
      qx(find ((ang <= pi/4) & (ang > 0))) = NR - OFF;
      qx(find ((ang <= 0) & (ang > -pi/4))) = OFF+1;
      ix(find ((ang <= -pi/4) & (ang > -pi/2))) = NR - OFF;
      ix(find ((ang <= -pi/2) & (ang > -3*pi/4))) = OFF+1;
      qx(find (ang <= -3*pi/4)) = OFF+1;

      a(indx) = layout(size (layout, 1) * (max (min (round (qx), ...
                       size (layout, 2)), 1)-1) + max (min (round (ix), size (layout, 1)), 1));
    endif
  unwind_protect_cleanup
    warning (wfi);
  end_unwind_protect

endfunction

%!function dec = __fntestqask1__ (msg, m)
%! [inp, qudr] = qaskenco (msg, m);
%! dec = qaskdeco (inp, qudr, m);

%!function __fntestqask2__ (m, dims)
%! msg = floor (rand (dims) * m);
%! assert (__fntestqask1__ (msg, m), msg);

%!test __fntestqask2__ (2, [100, 100])
%!test __fntestqask2__ (4, [100, 100])
%!test __fntestqask2__ (8, [100, 100])
%!test __fntestqask2__ (16, [100, 100])
%!test __fntestqask2__ (32, [100, 100])
%!test __fntestqask2__ (64, [100, 100])

%!test __fntestqask2__ (2, [100, 100, 3])
%!test __fntestqask2__ (4, [100, 100, 3])
%!test __fntestqask2__ (8, [100, 100, 3])
%!test __fntestqask2__ (16, [100, 100, 3])
%!test __fntestqask2__ (32, [100, 100, 3])
%!test __fntestqask2__ (64, [100, 100, 3])

%% Test input validation
%!error qaskdeco ()
%!error qaskdeco (1)
%!error qaskdeco (1, 2, 3, 4, 5)
%!error qaskdeco (1, 1)
%!error qaskdeco (1, 5)