## Copyright (C) 2018 P Sudeepam
##
## 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; see the file COPYING. If not, see
## <https://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{y} =} modulate (@var{x}, @var{fc}, @var{fs})
## @deftypefnx {Function File} {@var{y} =} modulate (@var{x}, @var{fc}, @var{fs}, "method")
## @deftypefnx {Function File} {@var{y} =} modulate (@var{x}, @var{fc}, @var{fs}, "method", @var{opt})
## @deftypefnx {Function File} {[@var{y}, @var{t_out}] =} modulate (@var{x}, @var{fc}, @var{fs}, "method", @var{opt})
##
## Modulate the message signal {@var{x}}, sampled at a sampling frequency
## {@var{fs}}, with a carrier frequency {@var{fc}}, using the modulation
## scheme specified by 'method'. {@var{opt}} is an extra parameter whose
## purpose depends on the modulation scheme choosen.
##
## If X is a matrix, its columns are modulated.
##
## [@var{y}, @var{t_out}] = MODULATE(...) returns the modulated signal and a time
## vector of the same length as Y. The time vector contains the time points
## where a sample of the signal has been taken.
##
## fs must satisfy fs > 2*fc
##
## @seealso{demod}
## @end deftypefn


function [y, t_out] = modulate (x, fc, fs, method, opt)

  ## Input argument validation
  if (nargin < 3 || nargin > 5)
    print_usage ();
  endif

  ## Checks for an Aliased carrier
  if (fc >= fs/2)
  	error ("modulate: Aliased carrier wave");
  endif

  ## Sets the default modulation scheme.
  if (nargin < 4)
      method = 'am';
  endif

  ## r is the number of rows, and c is the number of columns in X
  [r, c] = size (x);

  ## Returns an empty matrix for a empty signal input
  if (r*c == 0)
      y = [];
      return
  endif

  ## Convert row vector to column vector
  if (r == 1)
      x = x(:);
      len = c;
  else
      len = r;
  endif

  ## Modulation according to the METHOD specified

  switch (method)

    case {"am" "amdsb-sc"}
      ### Calculates the total number of time samples in {@var{x}}.
      t = (0: 1/fs: ((len-1)/fs))';
      ### Replecates the time samples for each column for the column wise
      ### modulation of X when it is of a matrix type.
      t = t(:, ones (1, size (x, 2)));
      ### Modulates the signal with the amdsb-sc modulation scheme.
      y = x .* cos (2*pi*fc*t);

    case {"amdsb-tc"}
      t = (0: 1/fs: ((len-1)/fs))';
      t = t(:, ones (1, size (x, 2)));
      if (nargin == 5)
          offset = opt;
      else
          offset = min (min (x)); ### Default offset value in MATLAB
      endif
      y = (x - offset) .* cos (2*pi*fc*t);

    case {"amssb"}
      t = (0: 1/fs: ((len-1)/fs))';
      t = t(:, ones (1, size (x, 2)));
      y = x .* cos (2*pi*fc*t) + imag (hilbert (x)) .* sin (2*pi*fc*t);

    case {"fm"}
      t = (0: 1/fs: ((len-1)/fs))';
      t = t(:, ones (1, size (x, 2)));
      if (nargin == 5)
          kf = opt;
      else
          x_max = max (max (abs (x)));
          ### If opt is not passed, by default, MATLAB chooses a non-negative
          ### kf value equal to (fc/fs)*2*pi/(max(max(x)))
          if (x_max > 0)
              kf = (fc / fs) * 2 * pi / x_max;
          else
              kf = 0;
          endif
      endif
      y = cos (2*pi*fc*t + kf*cumsum (x)); # Approximated rectangular integral

    case {"pm"}
      t = (0: 1/fs: ((len-1)/fs))';
      t = t(:, ones (1, size (x, 2)));
      if (nargin == 5)
          kp = opt;
      else
          kp = pi / max (max (abs (x))); # MATLAB default
      endif
      y = cos (2*pi*fc*t + kp*x);

      case {"qam"}
        if (nargin < 5)
          error ("modulate: two message signals are required")
        endif
        x2 = opt;
        r2 = size (x2, 1);
        if (r2 == 1)
          x2 = x2';
        endif
        if (any (size (x2) != size (x)))
          error ("modulate: message signals, dimensions mismatch")
        endif
        t = (0: 1/fs: ((len-1)/fs))';
        t = t(:, ones (1, size (x, 2)));
        y = x .* cos (2*pi*fc*t) + x2 .* sin (2*pi*fc*t);

    case {"pwm"}
      ### For PWM, all the elements of X must be between 0 and 1,
      ### specifying the width of each pulse in fractions of a period.
      x_max = max (max (x));
      x_min = min (min (x));
      if ( x_max > 1 || x_min < 0 )
        error ("modulate: Invalid Range, X must lie in [0,1]")
      endif

      ### https://www.mathworks.com/help/physmod/sps/ref/pwmgenerator.html

      ### fc/fs = number of cycles per sample
      ### fs/fc = number of samples per cycle

      ### 't' holds an index for each time sample, for the 1st time sample, it
      ### has a value of 0, for the 2nd time sample, a value of 1 and so on...
      t = (0: (fs/fc*len)-1)';
      if (nargin < 5)
          opt = 'left';
      endif

      tval =  rem (t, (fs/fc)) / (fs/fc); ### contains the required time value for each time sample.
      ind = floor ((t * (fc/fs)) + 1); ### contains the values of x repeated for each time sample as required

      if ( strcmpi ('centered', opt) )
          x(len + 1, :) = zeros (1, size (x, 2));  # Adds a zero at end of input
          y = zeros (length (t), size (x, 2));
          for i = 1 : size (x, 2)   # each iteration pulse width modulates a column
              y(:,i) = (tval < x(ind, i) / 2)  +  ((1 - tval) <= x(ind + 1, i) / 2);
          endfor
      elseif ( strcmpi ('left', opt) )
          y = zeros (length (t), size (x, 2));
          for i = 1 : size (x, 2)   # each iteration pulse width modulates a column
              y(:, i) = tval < x(ind, i);
          endfor
      else
          error ('modulate: PWM wave must be left or center aligned')
      endif
      t = t/fs;

    case {"ptm" "ppm"}
      x_max = max (max (x));
      x_min = min (min (x));

      if ( x_max > 1 || x_min < 0 )
        error ("modulate: Invalid Range, X must lie in [0,1]")
      endif
      if ( nargin < 5 )
        opt = .1; ### default MATLAB width
      endif

      if ( opt >= 1 || opt <= 0 )
        error ("modulate: Invalid Range, OPT must lie in (0,1)")
      endif

      if (x_max > 1-opt)
        warning("modulate: Pulse Overlapped")
      endif

      t = (0: (fs/fc*len)-1)';
      t = t/fs;
      y = zeros (length (t) + ceil (opt*fs/fc), size (x,2));
      temp = 0 : (opt * fs/fc - 1);
      for i = 1 : size (x, 2)    # each iteration pulse position modulates a column
          ind = ceil ( (x(:, i) + (0 : len-1)') * fs/fc ) + 1;
          ind = ind(:, ones (1, floor (opt* fs/fc))) + temp(ones (1, length (ind)), :);
          y(ind(:), i) = ones (size (ind(:)));
      endfor

  endswitch

  ### Modulation done...


  t = t(:, 1);   ### The first column of time vector is sufficient

  if (r == 1)   ### converts the column vectors y and t to a row vector if the input is in the form of row vector
    y = y';
    t = t';
  endif

  if (nargout == 2)
    t_out = t;
  endif

endfunction
