gsignal: Signal Processing

Documented in pulstran

# pulstran.R
# Copyright (C) 2019 Geert van Boxtel <gjmvanboxtel@gmail.com>
# Matlab/Octave signal package:
# Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
#
# 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/>.
#
# 20191126 Geert van Boxtel          First version for v0.1.0
#------------------------------------------------------------------------------

#' Pulse train
#'
#' Generate a train of pulses based on samples of a continuous function.
#'
#' Generate the signal \code{y <- sum(func(t + d, ...))} for each \code{d}. If
#' \code{d} is a matrix of two columns, the first column is the delay \code{d}
#' and the second column is the amplitude \code{a}, and \code{y <- sum(a *
#' func(t + d))} for each \code{d, a}. Clearly, \code{func} must be a function
#' which accepts a vector of times. Any extra arguments needed for the function
#' must be tagged on the end.
#'
#' If instead of a function name you supply a pulse shape sampled at frequency
#' \code{fs} (default 1 Hz), an interpolated version of the pulse is added at
#' each delay \code{d}. The interpolation stays within the the time range of the
#' delayed pulse. The interpolation method defaults to linear, but it can be any
#' interpolation method accepted by the function \code{interp1}
#'
#' @param t Time values at which \code{func} is evaluated, specified as a
#'   vector.
#' @param d Offset removed from the values of the array \code{t}, specified as a
#'   real vector, matrix, or array. You can apply an optional gain factor to
#'   each delayed evaluation by specifying \code{d} as a two-column matrix, with
#'   offset defined in column 1 and associated gain in column 2. If you specify
#'   \code{d} as a vector, the values are interpreted as delays only.
#' @param func Continuous function used to generate a pulse train based on its
#'   samples, specified as 'rectpuls', 'gauspuls', 'tripuls', or a function
#'   handle. If you use \code{func} as a function handle, you can pass the
#'   function parameters as follows:\cr \code{y <- pulstran(t, d, 'gauspuls',
#'   10e3, bw = 0.5)}.\cr This creates a pulse train using a 10 kHz Gaussian
#'   pulse with 50\% bandwidth. Alternatively, \code{func} can be a prototype
#'   function, specified as a vector. The interval of the function \code{0} to
#'   \code{(length(p) - 1) / fs}, and its samples are identically zero outside
#'   this interval. By default, linear interpolation is used for generating
#'   delays.
#' @param fs Sample rate in Hz, specified as a real scalar.
#' @param method Interpolation method, specified as one of the following
#'   options:
#' \describe{
#'   \item{"linear" (default)}{Linear interpolation. The interpolated value at a
#'   query point is based on linear interpolation of the values at neighboring
#'   grid points in each respective dimension. This is the default interpolation
#'   method.}
#'   \item{"nearest"}{Nearest neighbor interpolation. The interpolated value at a
#'   query point is the value at the nearest sample grid point.}
#'   \item{"cubic"}{Shape-preserving piecewise cubic interpolation. The
#'   interpolated value at a query point is based on a shape-preserving
#'   piecewise cubic interpolation of the values at neighboring grid points.}
#'   \item{"spline"}{Spline interpolation using not-a-knot end conditions. The
#'   interpolated value at a query point is based on a cubic interpolation of
#'   the values at neighboring grid points in each respective dimension.}
#' }
#' Interpolation is performed by the function \code{'interp1'} function in the
#' library \code{'pracma'}, and any interpolation method accepted by the
#' function \code{'interp1'} can be specified here.
#' 
#' @param ... Further arguments passed to \code{func}.
#'
#' @return Pulse train generated by the function, returned as a vector.
#'
#' @examples
#'
#' ## periodic rectangular pulse
#' t <- seq(0, 60, 1/1e3)
#' d <- cbind(seq(0, 60, 2), sin(2 * pi * 0.05 * seq(0, 60, 2)))
#' y <- pulstran(t, d, 'rectpuls')
#' plot(t, y, type = "l", xlab = "Time (s)", ylab = "Waveform",
#'      main = "Periodic rectangular pulse")
#'
#' ## assymetric sawtooth waveform
#' fs <- 1e3
#' t <- seq(0, 1, 1/fs)
#' d <- seq(0, 1, 1/3)
#' x <- tripuls(t, 0.2, -1)
#' y <- pulstran(t, d, x, fs)
#' plot(t, y, type = "l", xlab = "Time (s)", ylab = "Waveform",
#'      main = "Asymmetric sawtooth waveform")
#'
#' ## Periodic Gaussian waveform
#' fs <- 1e7
#' tc <- 0.00025
#' t <- seq(-tc, tc, 1/fs)
#' x <- gauspuls(t, 10e3, 0.5)
#' plot(t, x, type="l", xlab = "Time (s)", ylab = "Waveform",
#'      main = "Gaussian pulse")
#' ts <- seq(0, 0.025, 1/50e3)
#' d <- cbind(seq(0, 0.025, 1/1e3), sin(2 * pi * 0.1 * (0:25)))
#' y <- pulstran(ts, d, x, fs)
#' plot(ts, y, type = "l", xlab = "Time (s)", ylab = "Waveform",
#'      main = "Gaussian pulse train")
#'
#' # Custom pulse trains
#' fnx <- function(x, fn) sin(2 * pi * fn * x) * exp(-fn * abs(x))
#' ffs <- 1000
#' tp <- seq(0, 1, 1/ffs)
#' pp <- fnx(tp, 30)
#' plot(tp, pp, type = "l",xlab = 'Time (s)', ylab = 'Waveform',
#'      main = "Custom pulse")
#' fs <- 2e3
#' t <- seq(0, 1.2, 1/fs)
#' d <- seq(0, 1, 1/3)
#' dd <- cbind(d, 4^-d)
#' z <- pulstran(t, dd, pp, ffs)
#' plot(t, z, type = "l", xlab = "Time (s)", ylab = "Waveform",
#'      main = "Custom pulse train")
#'
#' @author Sylvain Pelissier, \email{sylvain.pelissier@@gmail.com}.\cr
#' Conversion to R by Geert van Boxtel \email{G.J.M.vanBoxtel@@gmail.com}.
#'
#' @export

pulstran <- function(t, d, func, fs = 1,
                     method = c("linear", "nearest", "cubic", "spline"), ...) {

  if (missing(t) || length(t) <= 0)
    stop("t must be an array")

  y <- rep(0L, length(t))
  if (missing(d) || length(d) <= 0) return(y)
  if (is.vector(d) || (is.array(d) && length(dim(d)) == 1)) {
    a <- rep(1L, length(d))
  } else if ((is.matrix(d) || is.array(d)) && ncol(d) == 2) {
    a <- d[, 2]
    d <- d[, 1]
  } else stop("invalid value specified for d")

  if (is.character(func)) {
    for (i in seq_along(d)) {
      y <- y + a[i] * do.call(func, list(t - d[i], ...))
    }
  } else {
    method <- match.arg(method)
    span <- (length(func) - 1) / fs
    t_pulse <- (0:(length(func) - 1)) / fs
    for (i in seq_along(d)) {
      dt <- t - d[i]
      idx <- which(dt >= 0 & dt <= span)
      if (length(idx) > 0) {
        y[idx] <- y[idx] + a[i] *
          pracma::interp1(t_pulse, func, dt[idx], method)
      }
    }
  }
  y
}