In bartekkroczek/detrendeR: Remove linear or polynomial trend from signal

knitr::opts_chunk$set(
     collapse = TRUE,
     comment = "#>",
     fig.width = 7.1,
     fig.height = 5
)

This package provides a small set of functions that detect and remove linear or polynomial trends in the time-series signal. The main functions are 'getltrend' and 'gettrend'. This package does not depend on any other package or function outside of the 'stat' namespace.

Installation

library(devtools)
devtools::install_github("bartekkroczek/detrendeR")

# or build with documentation
devtools::install_github("bartekkroczek/detrendeR", build_vignettes = TRUE)

Rationals

There is often to see a trend in the time series. It can be understood as a systematic pattern (e.g. rising or falling) that disturbs the stationarity of the signal. Removing this trend is important for many reasons, including to emphasize the oscillatory nature of the signal.

n <- 48 # no. of data points
time <- 1 # how many sec artifical signal lasts
freq <- 6 # frequency of miningful oscilaltion
t <- seq(0, time, length.out = n)
set.seed(123) # make runif getting this same result every time
noise <- runif(n)
trend <- seq(0.01, 3, length.out = n)
sine <- sin(2 * pi * freq * t) / 2
signal <- sine + trend + noise

f <- seq_along(t) / time
y <- fft(signal)
mag <- sqrt(Re(y)^2 + Im(y)^2) * 2 / n

layout(matrix(c(1, 2), 2, 1, byrow = TRUE))
plot(t, signal,
     type = "l", xlab = "Time [s]", ylab = "sine + trend + noise",
     main = "Example signal containing 6 Hz oscillation, trend and noise"
)
plot(f[1:length(f) / 2], mag[1:length(f) / 2],
     type = "l", xlab = "Frequency [Hz]", ylim = c(0, 0.5),
     ylab = "Amplitude", main = "Frequency spectra, oscillation is hard to spot"
)

Although the oscillation in the signal is quite distinct, the trend makes it barely noticeable in the spectrum.

A popular method of getting rid of a trend is to "filter" the signal with a moving average. However, this can lead to frequency artifacts, as shown in (@placeholder).

roll_mean <- (function(x, n = 5) filter(x, rep(1 / n, n), sides = 2))
roll_mean_ord  <- 5 # how many points averaged

mov_avg <- na.omit(signal - roll_mean(signal, roll_mean_ord))

f <- seq_along(mov_avg) / time
y <- fft(mov_avg)
mag <- sqrt(Re(y)^2 + Im(y)^2) * 2 / n

layout(matrix(c(1, 2), 2, 1, byrow = TRUE))
plot(seq_along(mov_avg) + floor(roll_mean_ord / 2), mov_avg,
     type = "l", xlab = "Time [s]", ylab = "sine + noise + artifacts?",
     main = "Effect of moving avg trend removal method"
)
plot(f[1:length(f) / 2], mag[1:length(f) / 2],
     type = "l", xlab = "Frequency [Hz]", ylim = c(0, 0.5), ylab = "Amplitude",
     main = "The desired effect is clearer, but there are also potential artifacts"
)
arrows(x0 = 15.1, y0 = 0.37, x1 = 9.1, y1 = 0.17, lwd = 0.8)
arrows(x0 = 15.4, y0 = 0.37, x1 = 12.05, y1 = 0.185, lwd = 0.8)
arrows(x0 = 15.7, y0 = 0.37, x1 = 16, y1 = 0.165, lwd = 0.8)
arrows(x0 = 16.0, y0 = 0.37, x1 = 17.95, y1 = 0.13, lwd = 0.8)
text(x = 15.6, y = 0.43, label = "Artifacts?")
arrows(x0 = 3.9, y0 = 0.38, x1 = 6.8, y1 = 0.22)
text(x = 3.75, y = 0.43, label = "Desired effect")

As an alternative, we propose a small set of functions that allows to detect and remove a linear or polynomial trend in the time series without contaminating signal.

library(detrendeR)

detrend <- detrendeR::getltrend(signal)

f <- seq_along(t) / time
y <- fft(detrended(detrend))
mag <- sqrt(Re(y)^2 + Im(y)^2) * 2 / n

layout(matrix(c(1, 2), 2, 1, byrow = TRUE))
plot(t, detrended(detrend),
     type = "l", xlab = "Time [s]", ylab = "sine + noise",
     main = "Signal detrended with detrendeR::getltrend"
)
plot(f[1:length(f) / 2], mag[1:length(f) / 2],
     type = "l", xlab = "Frequency [Hz]", ylim = c(0, 0.5),
     ylab = "Amplitude", main = "The desired effect is sound, artifacts are neglected"
)

Usage

Random data

n <- 48 # no. of data points
time <- 1 # how many sec artifical signal lasts
freq <- 6 # frequency of miningful oscilaltion
t <- seq(0, time, length.out = n)
noise <- runif(n)
trend <- seq(0.01, 3, length.out = n) + (sin(2 * pi * 2 * t) / 2)
sine <- sin(2 * pi * freq * t) / 2
signal <- sine + trend + noise

Remove and plot linear trend

detrend <- detrendeR::getltrend(signal)
plot(detrend)

Remove polynomial trend by setting a degree

detrend <- detrendeR::gettrend(signal, degree = 3)
plot(detrend)

... or set degree automatically

detrend <- detrendeR::gettrend(signal, maxfreq = 3, time = 1)
plot(detrend)

Get detrended signal

detrend <- detrendeR::gettrend(signal, degree = 5)

time  <- time(detrend) # same as time  <- t
orig  <- orisignal(detrend) # same as orig  <- signal
res <- detrended(detrend)

layout(matrix(c(1, 2), 2, 1, byrow = TRUE))
plot(time, orig, type = "l", xlab = "Time [s]", ylab = "Val", main = "Original signal")
plot(time, res, type = "l", xlab = "Time [s]", ylab = "Val", main = "Detrended signal")