## -*- texinfo -*-
## @deftypefn {} {y =} vco (@var{x}, @var{fc}, @var{fs})
## @deftypefnx {} {y =} vco (@var{x}, [@var{fmin}, @var{fmax}], @var{fs})
## @end deftypefnx

function [y,yn] = vco (x, fc = [], fs = [])

  if (nargin < 1 || nargin > 3)
    print_usage ();
  endif

  ## We can't just set default values in parameter list because
  ## default FC depends on value of FS and empty arguments may be
  ## used to request default values.

  if (isvector (x) && rows (x) == 1)
    x = x(:);
  endif

  if (any (x(:) > 1 | x(:) < -1))
    error ("vco: values in X must be in the range [-1, 1]");
  endif

  if (isempty  (fs))
    fs = 1;
  endif

  if (isempty (fc))
    fc = fs / 4;
  endif

  x_max = max (max (abs (x)));

  if (numel (fc) == 2)
    fmin = fc(1);
    fmax = fc(2);
    if (fmin > fs/2 || fmax > fs/2)
      error ("vco: FMIN and FMAX must be less than FS/2");
    endif
    if (fmin >= fmax)
      error ("vco: FMAX must be greater than FMIN");
    endif
    if (fmin < 0 || fmax < 0)
      error ("vco: FMIN and FMAX must be positive values");
    endif
    ## I don't know whether this is the best way, but I think it gives
    ## the proper value for fc and the proper deviation (+/-) from
    ## that value.
    fc = (fmin + fmax) / 2;
    kf = (fmax - fmin) / 2 / fs * 2 * pi / x_max;
  elseif (isscalar (fc))
    kf = (fc / fs) * 2 * pi / x_max;
  else
    error ("vco: FC must be a scalar");
  endif

  if (! (isscalar (fs)))
    error ("vco: FS must be a scalar");
  endif

  ## Note, no need to generate a matrix for T, we can rely on
  ## broadcasting instead.
  len = rows (x);
  t = (0 : (1/fs) : ((len-1)/fs))';

  ## Use modulate instead?
  y = cos (2*pi*fc*t + kf*cumsum (x));

endfunction
