R/filters.R

In soundgen: Sound Synthesis and Acoustic Analysis

#' Bandpass/stop filters
#'
#' Filtering in the frequency domain with FFT-iFFT: low-pass, high-pass,
#' bandpass, and bandstop filters with flexible filter depth (dB) and adjustable
#' bandwidths at filter boundaries. Similar to the "ffilter" function from the
#' seewave package, but here we use FFT instead of STFT - that is, the entire
#' sound is processed at once. This works best for relatively short sounds
#' (seconds), but gives us maximum precision (e.g., for precise notch filtering)
#' and doesn't affect the attack and decay. NAs are accepted and can be
#' interpolated or preserved in the output. Because we don't do STFT,
#' arbitrarily short vectors are also fine as input - for example, we can apply
#' a low-pass filter prior to decimation when changing the sampling rate without
#' aliasing. Note that, unlike \code{\link{pitchSmoothPraat}}, \code{bandpass}
#' by default applies an abrupt cutoff instead of a smooth gaussian filter, but
#' this behavior can be adjusted with the \code{bw} argument.
#'
#' Algorithm: fill in NAs with constant interpolation at the edges and linear
#' interpolation in the middle; perform FFT; set the frequency ranges to be
#' filtered out to 0; perform inverse FFT; set to the original scale; put the
#' NAs back in.
#'
#' @inheritParams .roxygen_defaults
#' @param lwr,upr cutoff frequencies, Hz. Specifying just lwr gives a high-pass
#'   filter, just upr low-pass filter with action = 'pass' (or vice versa with
#'   action = 'stop'). Specifying both lwr and upr a bandpass/bandstop filter,
#'   depending on 'action'
#' @param action "pass" = preserve the selected frequency range (bandpass),
#'   "stop" = remove the selected frequency range (bandstop)
#' @param dB a positive number giving the strength of effect in dB (defaults to
#'   Inf - complete removal of selected frequencies)
#' @param bw bandwidth of the filter cutoffs, Hz. Defaults to 0 (abrupt, step
#'   function - avoid!), a positive number corresponds to the standard deviation
#'   of a Gaussian curve, and two numbers set different bandwidths for the lower
#'   and upper cutoff points
#' @param na.rm if TRUE, NAs are interpolated; otherwise, they are preserved in
#'   the output
#' @param normalize if TRUE, resets the output to the original scale (otherwise
#'   filtering often reduces the amplitude)
#' @param ... other graphical parameters passed to \code{plot()} and
#'   \code{\link{meanSpectrum}}
#' @return Numeric vector of the same length as input.
#' @export
#' @examples
#' # Filter white noise
#' s1 = fade(c(rnorm(2000, 0, 1), NA), samplingRate = 16000)
#'
#' # low-pass
#' s1_lp = bandpass(s1, 16000, upr = 2000, bw = 100, plot = TRUE)
#' meanSpectrum(s1_lp, 16000)
#'
#' # high-pass by 40 dB
#' s1_hp = bandpass(s1, 16000, lwr = 2000, dB = 40, bw = 50,
#'   plot = TRUE, windowLength = 5)
#' # windowLength is passed to meanSpectrum for plotting (ignore the warnings)
#'
#' # bandstop
#' bandpass(s1, 16000, lwr = 1000, upr = 1800, bw = 50, action = 'stop', plot = TRUE)
#'
#' # bandpass (turns noise into a whistle)
#' s1_bp = bandpass(s1, 16000, lwr = 2000, upr = 2100, bw = 25, plot = TRUE)
#' # playme(rep(s1_bp, 5))
#' # spectrogram(s1_bp, 16000)
#'
#' # low-pass and interpolate a short vector with some NAs
#' x = rnorm(150, 10) + 3 * sin((1:50) / 5)
#' x[sample(seq_along(x), 50)] = NA
#' plot(x, type = 'l')
#' x_bandp = bandpass(x, samplingRate = 100, upr = 10, bw = 1)
#' points(x_bandp, type = 'l', col = 'blue')
#'
#' # boost or attenuate a particular frequency band by 20 dB
#' s3 = bandpass(s1, 16000, lwr = 1700, upr = 2100, bw = 200,
#'   dB = 20, plot = TRUE)
#' # spectrogram(s3, 16000)
#' s4 = bandpass(s1, 16000, lwr = 2000, upr = 4300, bw = c(100, 500),
#'   dB = 60, action = 'stop', plot = TRUE)
#' # spectrogram(s4, 16000)
#'
#' # precise notch filtering is possible, even in low frequencies
#' whiteNoise = rnorm(16000)
#' s5 = bandpass(whiteNoise, 16000, lwr = 30, upr = 40, normalize = TRUE,
#'               plot = TRUE, xlim = c(0, 0.5))
#' spectrum(s5, 16000, xlim = c(0, 0.1))
#' # playme(s5)
#' # spectrogram(s5, 16000, windowLength = 150, yScale = 'log')
#'
#' \dontrun{
#' # high-pass all audio files in a folder
#' bandpass('~/Downloads/temp', lwr = 2000, saveAudio = TRUE, savePlots = TRUE)
#' }
bandpass = function(
    x,
    samplingRate = NULL,
    lwr = NULL,
    upr = NULL,
    action = c('pass', 'stop'),
    dB = Inf,
    bw = 0,
    na.rm = TRUE,
    from = NULL,
    to = NULL,
    normalize = FALSE,
    reportEvery = NULL,
    cores = 1,
    saveAudio = FALSE,
    plot = FALSE,
    savePlots = FALSE,
    embed = FALSE,
    width = 900,
    height = 500,
    units = 'px',
    res = NA,
    ...
) {
  # match args
  action = match.arg(action)
  myPars = c(as.list(environment()), list(...))
  # exclude some args
  myPars = myPars[!names(myPars) %in% c(
    'x', 'samplingRate', 'reportEvery', 'cores', 'savePlots', 'saveAudio',
    'embed', 'width', 'height', 'units')]
  if (isTRUE(savePlots) && !saveAudio)
    message(paste('Saving plots, but not audio - the html notebook will play',
                  'the original audio; set saveAudio = TRUE to save the modified sounds'))
  pa = processAudio(x = x,
                    samplingRate = samplingRate,
                    from = from,
                    to = to,
                    saveAudio = saveAudio,
                    savePlots = savePlots,
                    funToCall = '.bandpass',
                    suffix = 'bandpass',
                    myPars = myPars,
                    reportEvery = reportEvery,
                    cores = cores)

  # htmlPlots
  if (isTRUE(savePlots) && pa$input$n > 1) {
    try(htmlPlots(pa$input, changesAudio = saveAudio,
                  width = paste0(width, units), embed = embed))
  }

  # prepare output
  if (pa$input$n == 1) {
    result = pa$result[[1]]
  } else {
    result = pa$result
  }
  invisible(result)
}


#' Bandpass filter per sound
#' @noRd
.bandpass = function(audio,
                     lwr = NULL,
                     upr = NULL,
                     action = 'pass',
                     dB = Inf,
                     bw = 0,
                     na.rm = TRUE,
                     normalize = FALSE,
                     plot = FALSE,
                     width = 900,
                     height = 500,
                     units = 'px',
                     res = NA,
                     ...) {
  if ((is.null(lwr) && is.null(upr)) || (!is.numeric(lwr) && !is.numeric(upr)))
    stop('Nothing to do: specify lwr and/or upr')
  if ((is.numeric(lwr) && is.numeric(upr)) && lwr > upr) {
    temp = lwr; lwr = upr; upr = temp
    warning('Found lwr >= upr; swapping them')
  }
  if (any(bw < 0)) stop('bw must be non-negative')
  x = audio$sound
  len = length(x)
  if (len < 3) return(x)

  mean_x = mean(x, na.rm = TRUE)
  ran_x = range(x, na.rm = TRUE)
  x = x - mean_x  # center to avoid a weird fade-out
  if (len > 10000) {
    # for long files, pad with 0 to the nearest power of 2 - much faster fft
    pad_with = 0  # ifelse(prod(ran_x) < 0, 0, mean(x))
    target_len = 2 ^ ceiling(log2(len))
    n_zeros = target_len - len
    x = c(x, rep(pad_with, n_zeros))
    len = len + n_zeros
  } else {
    n_zeros = 0
  }
  half_len = len %/% 2
  bin_width = audio$samplingRate / len

  # interpolate NAs
  idx_na = which(is.na(x))
  n_na = length(idx_na)
  if (n_na > 0) x = interpolateNA(x, idx_na = idx_na)

  # get spectrum
  sp = stats::fft(x)

  # bandwidths for gradual Gaussian slopes instead of step functions at lwr & upr
  # a = seq(-5, 5, .01)
  # b = pnorm(a, mean = 0, sd = 1)
  # plot(a, b)
  if (is.null(bw) || !is.numeric(bw)) {
    bw_lwr = bw_upr = 0
  } else {
    if (length(bw) == 1) {
      bw_lwr = bw_upr = bw
    } else {
      bw_lwr = bw[1]
      bw_upr = bw[2]
    }
  }
  bw_lwr_bins = round(bw_lwr / bin_width)
  bw_upr_bins = round(bw_upr / bin_width)

  ## set up the filter
  if (is.null(lwr)) idx_lwr = 1 else
    idx_lwr = max(1, min(half_len, round(lwr / bin_width)))
  if (is.null(upr)) idx_upr = half_len else
    idx_upr = max(1, min(half_len, round(upr / bin_width)))

  # lower
  a = 1:(half_len + 1)
  if (!is.null(lwr)) {
    if (bw_lwr_bins > 0) {
      b_lwr = 1 - pnorm(a, mean = idx_lwr, sd = bw_lwr_bins)
    } else {
      b_lwr = c(rep(1, idx_lwr - 1), rep(0, half_len + 1 - idx_lwr + 1))
    }
  } else {
    b_lwr = rep(0, half_len + 1)
  }
  # plot(b_lwr)

  # upper
  if (!is.null(upr)) {
    if (bw_upr_bins > 0) {
      b_upr = pnorm(a, mean = idx_upr, sd = bw_upr_bins)
    } else {
      b_upr = c(rep(0, idx_upr), rep(1, half_len + 1 - idx_upr))
    }
  } else {
    b_upr = rep(0, half_len + 1)
  }
  # plot(b_upr)

  # both combined
  if (action == 'pass') {
    filter = pmax(0, 1 - b_lwr - b_upr)
  } else {
    filter = b_lwr + b_upr
  }
  max_filt = max(filter)
  if (max_filt == 0) return(rep(0, audio$ls))
  filter = filter / max_filt
  # plot(filter)

  # normalize the filter
  if (is.finite(dB) && abs(dB) > .Machine$double.eps) {
    m = 10 ^ (-abs(dB) / 20)
    # normalize to (m, 1) - need to change min from 0 to m, preserving max at 1;
    # solve a system of equations (1+addid)/denom = 1 and (0+addit)/denom = m
    denom = 1 / (1 - m)
    addit = denom - 1
    filter = (filter + addit) / denom
  }
  # range(filter)
  # plot(a, filter)

  # mirror image of filter to cover the whole spectrum
  even = len %% 2 == 0
 if (even) {
    filter_full = c(filter, rev(filter[2:half_len]))
  } else {
    filter_full = c(filter, rev(filter[2:(half_len + 1)]))
  }
  # plot(filter_full)

  # inverse fft
  x_new = Re(fft(sp * filter_full, inverse = TRUE)) / len

  # trim extra zeros and fade-out (10 ms)
  if (n_zeros > 0) {
    x_new = .fade(
      list(sound = x_new[1:(len - n_zeros)],
           samplingRate = audio$samplingRate),
      fadeIn_points = 0,
      fadeOut_points = min(half_len, audio$samplingRate * .01))
  }

  # PLOTTING
  if (isTRUE(audio$savePlots)) {
    plot = TRUE
    png(filename = file.path(audio$path_output, paste0(audio$filename_noExt, ".png")),
        width = width, height = height, units = units, res = res)
    on.exit(dev.off())
  }
  if (plot) {
    par(mfrow = c(1, 2))
    on.exit(par(mfrow = c(1, 1)), add = TRUE)
    if (is.finite(dB)) {
      # plot in dB (log-scale)
      filter_dB = 20 * log10(pmax(filter, 1e-6))
      plot(a * bin_width / 1000, filter_dB, type = 'l',
           xlab = 'Frequency, kHz', ylab = 'Amplitude, dB', main = 'Filter', ...)
    } else {
      # plot on linear scale
      plot(a * bin_width / 1000, filter, type = 'l',
           xlab = 'Frequency, kHz', ylab = 'Amplitude', main = 'Filter', ...)
    }

    sp_old = try(.meanSpectrum(list(sound = audio$sound, samplingRate = audio$samplingRate),
                               yScale = 'max0', plot = FALSE, ...))
    sp_new = try(.meanSpectrum(list(sound = x_new, samplingRate = audio$samplingRate),
                               yScale = 'max0', plot = FALSE, ...))
    if (!inherits(sp_old, 'try-error') && any(!is.na(sp_old[, 2]))) {
      try(plot(sp_old, type = 'l', xlab = 'Frequency, kHz', ylab = 'dB', main = 'Spectra', ...))
      if (!inherits(sp_new, 'try-error') && any(!is.na(sp_new[, 2])))
        try(points(sp_new, type = 'l', col = 'blue', ...))
    } else {
      if (!inherits(sp_new, 'try-error') && any(!is.na(sp_new[, 2])))
        try(plot(sp_new, type = 'l', col = 'blue', xlab = 'Frequency, kHz',
                 ylab = 'dB', main = 'Spectra', ...))
    }
  }

  # back to original scale
  if (normalize && max(abs(x_new)) > 0)
    x_new = x_new / max(abs(x_new)) * audio$scale_used
  x_new = x_new + mean_x

  # save audio
  if (isTRUE(audio$saveAudio)) {
    filename = file.path(audio$path_output, paste0(audio$filename_noExt, ".wav"))
    writeAudio(x_new, audio = audio, filename = filename)
  }

  # put NAs back in
  if (!na.rm && n_na > 0) x_new[idx_na] = NA
  invisible(x_new)
}


#' Pitch smoothing as in Praat
#'
#' Smooths an intonation (pitch) contour with a low-pass filter, as in Praat
#' (http://www.fon.hum.uva.nl/praat/). Algorithm: interpolates missing values
#' (voiceless frames), performs FFT to obtain the spectrum, multiplies by a
#' Gaussian filter, performs an inverse FFT, and fills the missing values back
#' in. The \code{bandwidth} parameter controls the width of the low‑pass filter;
#' larger values retain more high‑frequency detail. Similar to a special case of
#' low-pass filtering with \code{\link{bandpass}}, but faster and can handle NAs.
#'
#' @seealso \code{\link{analyze}}
#'
#' @param pitch numeric vector of pitch values (NA = voiceless)
#' @param bandwidth the bandwidth of low-pass filter, Hz (high = less smoothing,
#'   close to zero = more smoothing)
#' @param samplingRate the number of pitch values per second
#' @param plot if TRUE, plots the original and smoothed pitch contours
#' @param preprocess if TRUE, makes input non-negative and handles NA; if FALSE,
#'   expects a non-negative vector without any missing values (faster)
#' @return Numeric vector of the same length as input.
#' @export
#' @examples
#' pitch = c(NA, NA, 405, 441, 459, 459, 460, 462, 462, 458, 458, 445, 458, 451,
#' 444, 444, 430, 416, 409, 403, 403, 389, 375, NA, NA, NA, NA, NA, NA, NA, NA,
#' NA, 183, 677, 677, 846, 883, 886, 924, 938, 883, 946, 846, 911, 826, 826,
#' 788, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 307,
#' 307, 368, 377, 383, 383, 383, 380, 377, 377, 377, 374, 374, 375, 375, 375,
#' 375, 368, 371, 374, 375, 361, 375, 389, 375, 375, 375, 375, 375, 314, 169,
#' NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 238, 285, 361, 374, 375, 375,
#' 375, 375, 375, 389, 403, 389, 389, 375, 375, 389, 375, 348, 361, 375, 348,
#' 348, 361, 348, 342, 361, 361, 361, 365, 365, 361, 966, 966, 966, 959, 959,
#' 946, 1021, 1021, 1026, 1086, 1131, 1131, 1146, 1130, 1172, 1240, 1172, 1117,
#' 1103, 1026, 1026, 966, 919, 946, 882, 832, NA, NA, NA, NA, NA, NA, NA, NA,
#' NA, NA)
#' pitchSmoothPraat(pitch, bandwidth = 10, samplingRate = 40, plot = TRUE)
#' pitchSmoothPraat(pitch, bandwidth = 2, samplingRate = 40, plot = TRUE)
pitchSmoothPraat = function(pitch,
                            bandwidth,
                            samplingRate,
                            plot = FALSE,
                            preprocess = TRUE) {
  len = length(pitch)
  if (len == 0) return(pitch)

  bin_width = samplingRate / len
  half_len = len %/% 2
  even = len %% 2 == 0

  if (preprocess) {
    orig = pitch

    # handle all-NA case early
    if (all(is.na(pitch))) return(pitch)

    # make non-negative
    m = min(pitch, na.rm = TRUE)
    pitch = pitch - m
    idx_unv = which(is.na(pitch))

    # interpolate NAs
    if (length(idx_unv) > 0) {
      pitch = interpolateNA(pitch, idx_na = idx_unv)
    }
  } else {
    idx_unv = integer(0)
  }

  # get spectrum
  sp = stats::fft(pitch)

  # gaussian filter
  # Frequencies for the first half of the spectrum (0 to Nyquist)
  freq = (0:half_len) * bin_width
  filter_half = exp(-(freq / bandwidth)^2)
  # plot(freq, filter_half / max(filter_half), type = 'l')

  # Mirror to cover the whole spectrum exactly
  if (even) {
    # For even length: len = 2 * half_len
    # We need indices 0 to half_len, then mirror half_len-1 down to 1
    filter_full = c(filter_half, rev(filter_half[2:half_len]))
  } else {
    # For odd length: len = 2 * half_len + 1
    # We need indices 0 to half_len, then mirror half_len down to 1
    filter_full = c(filter_half, rev(filter_half[2:(half_len + 1)]))
  }
  # plot((0:(len-1)) * bin_width, filter_full, type = 'l')

  # inverse fft
  pitch2 = Re(fft(sp * filter_full, inverse = TRUE)) / len
  if (preprocess) {
    if (length(idx_unv) > 0) pitch2[idx_unv] = NA
    pitch2 = pitch2 + m  # back to the original scale
  }

  if (plot) {
    if (preprocess) {
      plot(orig)
    } else {
      plot(pitch)
    }
    lines(pitch2, col = 'red')
  }

  pitch2
}