fats: Frequency Analysis of Time Series

#' Function plotting the filter frequency response.
#'
#' The function calculates and plots the frequency response of the given list of filters.
#'
#' This function takes a list of filters either developed via the signal package (i.e. fir1() or butter()),
#' or custom developed filters such as moving average, which need to have the class Ma added manually.
#' Several arguments allow customizing the plot generated by the function.
#' The main purpose is to provide the right scaling on the x-Axis which is often a problem for beginners.
#' Also, the filter response of Arma models can be calculated and plotted.
#'
#' @param filterList List of filters of class Ma or Arma; FIR or IIR filters as derived from package signal
#' @param timeStep Numeric value (default: 1); time between two samples.
#' @param resolution Numeric value (default: 512); defining frequency resolution of filter response.
#' @param xAxis Character string (default: "period"); Either "frequency"
#' or "period". Defining if on the x-axis the linear frequency or the reciprocal
#' oscillation period is displayed.
#' @param xAxisLimits Numeric vector of length 2 (default: c(0, 1));
#' Defines the range of the spectrum being displayed. 0 -> 0Hz/Inf Period, 1 -> Nyquist Frequency/Nyquist Period (twice the time step).
#' The range is the same as for band edges during filter design in package signal (i.e. fir1())
#' @param returnPlottingData Logical value; Should plot data be returned as data frame.
#'
#' @return Plots the filter response of all filters in argument filterList.
#' Returns a data frame in long format with all plotting information.
#'
#' @author Daniel Beiter, \email{daniel.beiter@@gfz-potsdam.de}
#' @keywords filter response
#'
#' @import signal
#' @import dplyr
#' @import tidyr
#' @import ggplot2
#' @importFrom magrittr %>%
#' @export
#'
#' @examples
#'
#' ### plot FIR filters with various cutoff frequencies
#' cutoffFrequencies <- seq(0.1, 0.9, 0.1)
#' lowPassFIR <- lapply(cutoffFrequencies, function(x) signal::fir1(n = 10, w = x, type = "low"))
#' names(lowPassFIR) <- paste0("Cutoff frequency=", cutoffFrequencies)
#' plotFilterResponse(lowPassFIR)
#'
#' ### plot filter repsonse of moving average with various window lengths
#' windowLength <- seq(2, 5)
#' movingAverage <- lapply(windowLength, function(x) {averageFilter <- rep(1/x, x); class(averageFilter) <- "Ma"; averageFilter})
#' names(movingAverage) <- c(paste0("Moving Average n=", windowLength))
#' plotFilterResponse(movingAverage)
#'
#' ### comparing moving average with FIR filter
#' movingAverage <- rep(1/10, 10)
#' class(movingAverage) <- "Ma"
#' filterList <- list(signal::fir1(n = 15, w = 0.01, type = "low"), movingAverage)
#' names(filterList) <- c("FIR_n=15_w=0.01", "MA_n=10")
#' plotFilterResponse(filterList)
plotFilterResponse <- function(filterList, timeStep = 1, resolution = 512, xAxis = "period", xAxisLimits = c(0, 1), returnPlottingData = FALSE) {

  ### testing arguments before proceeding
  if (class(filterList) != "list")
    stop("Argument \"filterList\" must be a list! Value used was: \"", class(filterList),  "\"!")
  if (is.null(names(filterList))) {
    warning("Argument \"filterList\" should have names and not be NULL! Names added.")
    names(filterList) <- paste0("Filter_", seq(filterList))
  }
  if (!all(sapply(filterList, class) %in% c("Ma", "Arma")))
    stop("Argument \"filterList\" must only contain elements of class \"Ma\" or \"Arma\"!")
  if (!(is.numeric(timeStep) & length(timeStep) == 1))
    stop("Argument \"timeStep\" must be single value!")
  if (!xAxis %in% c("frequency", "period"))
    stop("Use either \"frequency\" or \"period\" for xAxis! Value used was: \"", xAxis, "\"!")

  ### pre allocate list for use in for loop
  nFilters <- length(filterList)
  ggFilterSpectrum <- vector("list", nFilters)

  ### using for loop instead of lapply because list names are needed during creation of data frame
  ### only few loop runs necessary; advantage of lapply not significant
  for (filterIdx in 1:nFilters) {
    currentFilter <- filterList[[filterIdx]]

    ### actual calculation of frequency response
    filterResponse <- signal::freqz(filt = currentFilter, n = resolution, Fs = 1 / timeStep)

    ### set-up data frame containing filter spectrum
    filterSpectrum <- data.frame(frequency = filterResponse$f,
                                 period = 1 / filterResponse$f,
                                 complex = filterResponse$h,
                                 amplitude = Mod(filterResponse$h),
                                 phase = Arg(filterResponse$h))

    ### group delay must be calculated differently for FIR and IIR
    ### rounding because of numeric problems when using single sin signals
    if (class(currentFilter) == "Ma")
      filterSpectrum$groupDelay <- round(signal::grpdelay(filt = as.numeric(currentFilter))$gd, 1)
    if (class(currentFilter) == "Arma")
      filterSpectrum$groupDelay <- round(signal::grpdelay(filt = currentFilter)$gd, 1)

    ### instead of plotting phase, group delay is plotted
    ggFilterSpectrum[[filterIdx]] <- filterSpectrum %>%
      dplyr::select(frequency, period, amplitude, groupDelay) %>%
      tidyr::gather(key = "variable", value = "value", -c("frequency", "period"), factor_key = TRUE) %>%
      dplyr::mutate(filterName = names(filterList)[filterIdx])
  }
  ### converting list of data frames into single data frame
  ggFilterSpectrum <- ggFilterSpectrum %>%
    bind_rows() %>%
    mutate(filterName = factor(filterName))

  ### create ggplot with amplitude and group delay spectrum
  ### add horizontal lines in amplitude spectrum at 0 and 1
  pSpectrum <- ggplot2::ggplot() +
    ggplot2::geom_line(ggplot2::aes(x = frequency, y = value, colour = filterName), data = ggFilterSpectrum, size = 1.2) +
    ggplot2::geom_hline(ggplot2::aes(yintercept = yIntercept), data = data.frame(yIntercept = c(0, 1), variable = "amplitude"), size = 1.2, linetype = "dashed") +
    ggplot2::labs(title = "Filter response", x = paste0(toupper(substring(xAxis, 1, 1)), substring(xAxis, 2))) +
    ggplot2::coord_cartesian(xlim = xAxisLimits / (2 * timeStep)) +
    ggplot2::facet_grid(variable ~ ., scales = "free_y", switch = "y", labeller = labeller(variable = c(amplitude = "Amplitude [-]", groupDelay = "Group Delay [samples]"))) +
    ggplot2::theme_bw() +
    ggplot2::theme(axis.title.y = element_blank(), plot.title = element_text(hjust = 0.5, size = 18, face = "bold")) +
    ggplot2::guides(colour = guide_legend(title = "Filters"))

  if (xAxis == "frequency") {
    pSpectrum <- pSpectrum +
      ggplot2::scale_x_continuous(labels = function(x) sprintf("%.2f", x))
  }
  if (xAxis == "period") {
    pSpectrum <- pSpectrum +
      ggplot2::scale_x_continuous(labels = function(x) sprintf("%.2f", 1 / x))
  }
  plot(pSpectrum)

  if (returnPlottingData)
    return(ggFilterSpectrum)
}