R/filter_commands.R

In craddm/firfiltR: FIR filtering

#' FIR Hamming windowed-sinc filter
#'
#' Filters data using an Hamming windowed-sinc FIR filer. An implementation of
#' Andreas Widmann's firfilt plugin for EEGLAB.
#'
#' @author Andreas Widmann. Ported to R by Matt Craddock
#'   \email{matt@@mattcraddock.com}
#'
#' @param data Data to be filtered.
#' @param low_edge Low passband edge. Specify ONLY this for a high-pass filter.
#' @param high_edge High passband edge. Specify ONLY this for a low-pass filter.
#' @param filt_order Length of the filter in samples. If not specified, will be
#'   automatically estimated.
#' @param notch Convert to a notch filter. Defaults to FALSE (bandpass filter).
#'   Only relevant if two cut-offs specified.
#' @param minphase Minimum-phase filter. Not yet implemented.
#' @param srate Sampling rate.
#'
#' @export

simple_filt <- function(data,
                        low_edge = NULL,
                        high_edge = NULL,
                        filt_order = NULL,
                        notch = FALSE,
                        minphase = FALSE,
                        srate = NULL) {
  # ensure minphase is set false for now, until this is implemented.
  minphase <- FALSE

  ## Check parameters ----------------
  if (is.null(low_edge) & is.null(high_edge)) {
    stop("At least one cutoff must be specified.")
  } else if (is.null(low_edge) | is.null(high_edge) && notch == TRUE) {
    stop("Both edges must be specified for notch filtering.")
  }

  if (is.null(srate)) {
    stop("Sampling rate must be specified.")
  }

  trans_width_ratio <- 0.25 # transition width ratio
  f_nyquist <- srate / 2

  # this flips the filter around if only a low cutoff is given.
  if (is.null(high_edge)) {
    high_edge <- low_edge
    low_edge <- NULL
    if (!notch) {
      notch <- TRUE
    }
  }

  edge_array <- sort(c(low_edge, high_edge))

  if (any(edge_array < 0) | any(edge_array >= f_nyquist)) {
    stop('Cutoff frequency out of range.')
  }

  if (!is.null(filt_order) &&
      (filt_order < 2 || filt_order %% 2 != 0)) {
    stop('Filter order must be a real, even, positive integer.')
  }

  # Max stop-band width
  max_TBW_array <- edge_array # Band-/highpass

  if (!notch) {
    # Band-/lowpass
    max_TBW_array[length(max_TBW_array)] <- f_nyquist - edge_array[length(edge_array)]

  } else if (length(edge_array) == 2) {
    # Bandstop
    max_TBW_array <- diff(edge_array) / 2
  }

  max_Df <- min(max_TBW_array)

  # Transition band width and filter order
  if (is.null(filt_order)) {
    # Default filter order heuristic

    if (notch) {

      # Highpass and bandstop
      df <- min(c(max(c(max_Df * trans_width_ratio, 2)), max_Df))

    } else {

      #% Lowpass and bandpass
      df <- min(c(max(c(edge_array[1] * trans_width_ratio, 2)), max_Df))
    }

    filt_order <- 3.3 / (df / srate)
    filt_order <- ceiling(filt_order / 2) * 2 # Filter order must be even.

  } else {

    df <- 3.3 / filt_order * srate # Hamming window
    filtorderMin <- ceiling(3.3 / ((max_Df * 2) / srate) / 2) * 2
    filtorderOpt <- ceiling(3.3 / (max_Df / srate) / 2) * 2

    if (filt_order < filtorderMin) {
      warning(
        'Filter order too low. Minimum required filter order is %d. For better results a minimum filter order of %d is recommended.',
        filtorderMin,
        filtorderOpt
      )
    } else if (filt_order < filtorderOpt) {
      warning(
        'firfilt:filterOrderLow',
        'Transition band is wider than maximum stop-band width. For better results a minimum filter order of %d is recommended. Reported might deviate from effective -6dB cutoff frequency.',
        filtorderOpt
      )
    }
  }

  if (notch) {
    if (is.null(low_edge)) {
      filter_type <- 'highpass'
      cutoff_array <- edge_array - (df / 2)
    } else {
      filter_type <- 'bandstop'
      cutoff_array <- c(edge_array[1] + df / 2, edge_array[2] - df / 2)
    }
  } else if (is.null(low_edge)) {
    filter_type <- 'lowpass'
    cutoff_array <- edge_array + df
  } else {
    filter_type <- 'bandpass'
    cutoff_array <- c(edge_array[1] - df / 2, edge_array[2] + df / 2)
  }

  message(
    paste("Performing", filt_order + 1, filter_type, "filtering.")
  )

  message(sprintf('simple_filt - transition band width: %.4g Hz\n', df))
  message(sprintf('simple_filt - passband edge(s): %s Hz\n', edge_array))
  message(sprintf(
    'simple_filt - cutoff frequency(ies) (-6 dB): %s Hz\n',
    cutoff_array
  ))

  if (notch) {
    b <- firws(filt_order,
               cutoff_array,
               filter_type,
               "hamming",
               srate = srate)

  } else {
    b <- firws(filt_order,
               cutoff_array,
               "hamming",
               srate = srate)
  }

  # if (minphase) {
  #   message('Converting filter to minimum-phase (non-linear!)')
  #   b <- minphaserceps(b)
  #   causal <- 1
  #   dir <- '(causal)'
  # } else {
  #   causal <- 0
  #   dir <- '(zero-phase, non-causal)'
  # }

  # Filter
  # if (minphase) {
    # EEG <- firfiltsplit(data, b, causal) # Not yet implemented
  # } else {
    EEG_new <- firfilt(data, b)
  # }
  return(EEG_new)
}


#' Perform filtering with pre-designed FIR filter
#'
#' @param x Data. Should be in wide format, one column per electrode.
#' @param b Filter coefficients
#' @importFrom future.apply future_lapply
#' @export

firfilt <- function(x, b) {

  # Filter's group delay
  group_delay <- (length(b) - 1) / 2

  #create a DC array
  if (is.null(dim(x))) {
    dc_array <- c(1,length(x))
    } else {
      dc_array <- c(1,dim(x)[[1]])
      }

  for (iDc in 1:(length(dc_array) - 1)) {
    y <- rep(0,length(x))
    x <- rbind(x[rep(1, group_delay) * dc_array[iDc], , drop = FALSE],
                   x,
                   x[rep(dim(x)[[1]], group_delay) * dc_array[iDc], , drop = FALSE])

    if (!is.null(ncol(x)) && ncol(x) > 1) {
      y <- future.apply::future_lapply(x,
                                       function(ii) stats::filter(ii,
                                                                  b / 1,
                                                                  sides = 1))
      y <- do.call("rbind", y)
    } else {
      y <- stats::filter(x, b / 1, sides = 1)
    }

    y <- na.omit(t(y), b/1, sides = 1)
  }

  filtered_df <- as.data.frame(y)
  names(filtered_df) <- names(x)
  return(filtered_df)
}

#' Design a windowed-sinc type I linear phase FIR filter
#'
#' @author Andreas Widmann. Ported to R by Matt Craddock \email{matt@@mattcraddock.com}
#'
#' @param m Filter order. Must be even. Can be derived automatically with "auto".
#' @param f Normalized cutoff frequency(ies).
#' @param t Filter type. Options = "highpass", "bandstop", "lowpass", "bandpass"
#' @param w Window type. Defaults to Hamming window.
#' @param beta Required for use of Kaiser windows.
#' @param srate Sampling rate is required.
#' @param tbw Required when filter order is set to "auto".
#' @return b filter  coefficients
#' @export

firws <- function(m = "auto",
                  f,
                  t = NULL,
                  w = "hamming",
                  beta = NULL,
                  srate,
                  tbw = NULL) {

  #a <- 1

  if (m == "auto") {
    if (is.null(tbw)) {
      stop("Transition bandwidth parameter required.")
    }

    m <- est_filt_order(type = w,
                        tbw = tbw,
                        srate = srate)

  } else if (m %% 2 == 1) {
    stop('Must be even.')
  }

  f <- f / 2

  w <- select_window(w, m)
  b <- fkernel(m,
               f[1],
               w)

  if (length(f) == 1 && t == "highpass") {
    b <- fspecinv(b)
  } else if (length(f) == 2) {
    b <- b + fspecinv(fkernel(m, f[2], w))
    if (is.null(t) || t != "bandstop") {
      b <- fspecinv(b)
    }
  }
  return(b)
}

#' Estimate filter order.
#'
#' @author Andreas Widmann. Ported to R by Matt Craddock \email{matt@@mattcraddock.com}
#' @param type Window type. e.g. Hamming, Kaiser
#' @param tbw Transition bandwidth.
#' @param pb_dev Max passband ripple/deviation. Note, only used for Kaiser windows.
#' @param srate Sampling rate.
#' @export
#'

est_filt_order <- function(type, tbw, pb_dev = .0002, srate) {

  df <- tbw / srate #normalize transition bandwidth
  if (type == "hamming") {
    filt_ord <- 3.3 / df
  } else if (type == "kaiser") {
    filt_ord <- 1 + ((-20 * log10(pb_dev) - 8) / (2.285 * 2 * pi * df))
  }
  filt_ord <- ceiling(filt_ord / 2) * 2
}


#' Custom filtering.
#'
#' This function supports use of different transition bandwidths for low- and
#' high-pass cutoffs respectively. Note that since it is intended to be used
#' carefully, it has few defaults.
#'
#' @author Matt Craddock \email{matt@@mattcraddock.com}
#' @param data Data to be filtered. Wide format expected.
#' @param f Cutoff frequency/ies
#' @param tbw Transition bandwidth(s).
#' @param win_type Window type.
#' @param beta Required if window type is Kaiser
#' @param pb_dev Required if window type is Kaiser
#' @param type Filtering type - "highpass", "lowpass", "bandpass", "bandstop"
#' @param srate Sampling rate (required if data is not of class "eeg_data").
#'

custom_filter <- function(data,
                          f,
                          tbw = NULL,
                          win_type = "hamming",
                          beta = NULL,
                          pb_dev = 2e-04,
                          type,
                          srate = NULL) {

  if (type == "bandpass" | type == "bandstop" && length(f) <2) {
    stop("Two cutoffs must be supplied to use bandpass or bandstop filters.")
    }

  f <- sort(f)
  tbw_ratio <- 0.25

  if (is.null(tbw)) {
    message("Estimating transition bandwidths.")
  #  df <- min(c(max(c(max_Df * tbw_ratio, 2)), max_Df))
  }

  # pb_dev is only required for Kaiser windows, but is ignored for other window types.
  first_order <- est_filt_order(win_type,
                                tbw = tbw[1],
                                pb_dev = pb_dev,
                                srate = srate)
  if (length(tbw == 2)) {
    second_order <- est_filt_order(win_type,
                                   tbw = tbw[2],
                                   pb_dev = pb_dev,
                                   srate = srate)
  }
}


#' Compute filter kernel.
#'
#' @author Andreas Widmann. Ported to R by Matt Craddock \email{matt@@mattcraddock.com}
#' @param m Filter order.
#' @param f Cutoff frequencies.
#' @param w Window and filter coefficients.
#' @keywords internal
fkernel <- function(m, f, w) {
  m <- (-m / 2):(m / 2)
  b <- rep(0,length(m))
  b[m == 0] <- 2 * pi * f # No division by zero
  b[m != 0] <- sin(2 * pi * f * m[m != 0]) / m[m != 0] # Sinc
  b <- b * w # Window
  b <- b / sum(b) # Normalization to unity gain at DC
  b
}

#' Spectral inversion.
#'
#' @author Andreas Widmann. Ported to R by Matt Craddock \email{matt@@mattcraddock.com}
#' @param b filter coefficients
#' @keywords internal
fspecinv <- function(b) {
  b <- -b
  b[(length(b) - 1) / 2 + 1] <- b[(length(b) - 1) / 2 + 1] + 1
  b
}



# rcepsminphase() - Convert FIR filter coefficient to minimum phase
#
# Usage:
#   >> b = minphaserceps(b);
#
# Inputs:
 #   b - FIR filter coefficients
#
# % Outputs:
#   %   bMinPhase - minimum phase FIR filter coefficients
# %
# % Author: Andreas Widmann, University of Leipzig, 2013
# %
# % References:
#   %   [1] Smith III, O. J. (2007). Introduction to Digital Filters with Audio
# %       Applications. W3K Publishing. Retrieved Nov 11 2013, from
# %       https://ccrma.stanford.edu/~jos/fp/Matlab_listing_mps_m.html
# %   [2] Vetter, K. (2013, Nov 11). Long FIR filters with low latency.
# %       Retrieved Nov 11 2013, from
# %       http://www.katjaas.nl/minimumphase/minimumphase.html

# #' Convert FIR filter coefficients to minimum-phase.
# #'
# #' @param b FIR filter coefficients
# #'

 #minphaserceps <- function(b) {
#
#   # Line vector
#   b <- b(:)'
#
#   n = length(b);
#   upsamplingFactor = 1e3; % Impulse response upsampling/zero padding to reduce time-aliasing
#   nFFT = 2^ceil(log2(n * upsamplingFactor)); % Power of 2
#   clipThresh = 1e-8; % -160 dB
#
#   % Spectrum
#   s = abs(fft(b, nFFT));
#   s(s < clipThresh) = clipThresh; % Clip spectrum to reduce time-aliasing
#
#   % Real cepstrum
#   c = real(ifft(log(s)));
#
#   % Fold
#   c = [c(1) [c(2:nFFT / 2) 0] + conj(c(nFFT:-1:nFFT / 2 + 1)) zeros(1, nFFT / 2 - 1)];
#
#   % Minimum phase
#   bMinPhase = real(ifft(exp(fft(c))));
#
#   % Remove zero-padding
#   bMinPhase = bMinPhase(1:n)
# }