R/wavelet-transform-250.R

In Teresa00/hfmAnnotation: HFM Internal Package for Annotations

# Wavelet transform calculation
# using Quadratic spline wavelet filterbank (Sampling Frequency 250Hz)
# Teresa Feb 21st 2018
#' @importFrom assertthat assert_that
#' @importFrom magrittr %>%
#' @importFrom stats sd
#'
library("magrittr")
library("assertthat")

length_synchronization <- function(length) {
  floor((length - 1) / 2)
}


# Normalize signal to scale [a, b]
normalization <- function(signal, a, b) {
  if(min(signal) == max(signal)){
    #((b - a) * (signal - min(signal))) + a
    rep(a, length(signal))
  }else{
    ((b - a) * (signal - min(signal))) / (max(signal) - min(signal)) + a
  }
}


generalized_scale_ <- function(signal){
  sigma_ <- sd(signal)
  sigma_[sigma_ == 0] <- 1

  (signal - mean(signal)) / sigma_
}


ecg_normalization <- function(signal) {
  assert_that(
    all(is.numeric(signal)),
    all(is.finite(signal))
  )

  # scale is the build-in function in R
  # standarize the signal to zero mean and unit variance
  scaled_signal <-
    generalized_scale_(signal) %>%
    as.vector()

  normalized_signal <-
    normalization(scaled_signal, -1, 1)

}


trim_vector_length <- function(q_i) {
  q <-
    q_i %>%
    round(10)

  q[1:max(which(q != 0))]
}


wavelet_transform_250 <- function(ecg) {
  assert_that(
    all(is.numeric(ecg)),
    all(is.finite(ecg))
  )

  # Initialization
  h <- 0.125 * c(1, 3, 3, 1)
  g <- 2 * c(1, -1)

  h2 <-
    h %>%
    kronecker(c(1, 0))

  h4 <-
    h %>%
    kronecker(c(1, rep(0, 3)))

  h8 <-
    h %>%
    kronecker(c(1, rep(0, 7)))

  g2 <-
    g %>%
    kronecker(c(1, 0))

  g4 <-
    g %>%
    kronecker(c(1, rep(0, 3)))

  g8 <-
    g %>%
    kronecker(c(1, rep(0, 7)))

  g16 <-
    g %>%
    kronecker(c(1, rep(0, 15)))

  q1 <-
    g %>%
    trim_vector_length()

  q2 <-
    h %>%
    signal::conv(g2) %>%
    trim_vector_length()

  q3 <-
    h %>%
    signal::conv(h2) %>%
    signal::conv(g4) %>%
    trim_vector_length()

  q4 <-
    h %>%
    signal::conv(h2) %>%
    signal::conv(h4) %>%
    signal::conv(g8) %>%
    trim_vector_length()

  q5 <-
    h %>%
    signal::conv(h2) %>%
    signal::conv(h4) %>%
    signal::conv(h8) %>%
    signal::conv(g16) %>%
    trim_vector_length()

  wt <- matrix(0, nrow = length(ecg), ncol = 5)

  wt[, 1] <- signal::filter(q1, 1, ecg)
  wt[, 2] <- signal::filter(q2, 1, ecg)
  wt[, 3] <- signal::filter(q3, 1, ecg)
  wt[, 4] <- signal::filter(q4, 1, ecg)
  wt[, 5] <- signal::filter(q5, 1, ecg)

  wavelets <- list(
    q1 = q1,
    q2 = q2,
    q3 = q3,
    q4 = q4,
    q5 = q5,
    wt = wt
  )
}


#' Preprocess ECG signal before detection and delineation.
#' Return processed ecg signal, wavelet transform matrix,
#' thresholds for maxima modulus search and "alignments" vector for index alignments
#' caused by wavelet transform filter bank
#'
#'
#' @param ecg ECG signal.
#' @param fs_initial The sampling frequency of original ECG signal.
#'
#' @export
#'


wavelet_preprocessing <- function(ecg, fs_initial) {
  fs_ratio <-
    250 / fs_initial

  ecg <-
    if(fs_ratio == 1){
      ecg %>% ecg_normalization()
    }else{
      ecg %>%
        signal::resample(fs_ratio, d = 1) %>%
        ecg_normalization()
    }

  # Wavelet transform calculation and synchronization
  wavelets_list <-
    wavelet_transform_250(ecg)

  d1 <-
    length(wavelets_list$q1) %>%
    length_synchronization()

  d2 <-
    length(wavelets_list$q2) %>%
    length_synchronization()

  d3 <-
    length(wavelets_list$q3) %>%
    length_synchronization()

  d4 <-
    length(wavelets_list$q4) %>%
    length_synchronization()

  l5 <-
    length(wavelets_list$q5)

  d5 <-
    length_synchronization(l5)

  wt <- wavelets_list$wt[-c(1:(l5 - 1)), ]

  wavelet_transforms <- matrix(0, nrow(wt) - d5, 5)

  wavelet_transforms[, 5] <- wt[(d5 + 1):nrow(wt), 5]
  wavelet_transforms[, 4] <- wt[(d4 + 1):(nrow(wt) + d4 - d5), 4]
  wavelet_transforms[, 3] <- wt[(d3 + 1):(nrow(wt) + d3 - d5), 3]
  wavelet_transforms[, 2] <- wt[(d2 + 1):(nrow(wt) + d2 - d5), 2]
  wavelet_transforms[, 1] <- wt[(d1 + 1):(nrow(wt) + d1 - d5), 1]

  thresholds <-
    0.5 * sqrt(apply(wavelet_transforms^2, 2, mean))

  thresholds[4] <-
    thresholds[4] * 2

  list(
    processed_ecg = ecg,
    wavelets = wavelet_transforms,
    thresholds = thresholds,
    alignments = l5:(length(ecg) - d5)
  )
}