inst/doc/gsignal.R

In gsignal: Signal Processing

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(gsignal)

## ----pulstran, fig.height=7, fig.width=7--------------------------------------

op <- par(mfrow = c(2, 2))

## periodic rectangular pulse
t <- seq(0, 60, 1/1e3)
d <- cbind(seq(0, 60, 2), sin(2 * pi * 0.05 * seq(0, 60, 2)))
y <- pulstran(t, d, 'rectpuls')
plot(t, y, type = "l", xlab = "", ylab = "",
     main = "Periodic rectangular pulse")

## assymetric sawtooth waveform
fs <- 1e3
t <- seq(0, 1, 1/fs)
d <- seq(0, 1, 1/3)
x <- tripuls(t, 0.2, -1)
y <- pulstran(t, d, x, fs)
plot(t, y, type = "l", xlab = "", ylab = "",
     main = "Asymmetric sawtooth ")

## Periodic Gaussian waveform
fs <- 1e7
tc <- 0.00025
t <- seq(-tc, tc, 1/fs)
x <- gauspuls(t, 10e3, 0.5)
ts <- seq(0, 0.025, 1/50e3)
d <- cbind(seq(0, 0.025, 1/1e3), sin(2 * pi * 0.1 * (0:25)))
y <- pulstran(ts, d, x, fs)
plot(ts, y, type = "l", xlab = "", ylab = "",
     main = "Gaussian pulse")

## Custom pulse trains
fnx <- function(x, fn) sin(2 * pi * fn * x) * exp(-fn * abs(x))
ffs <- 1000
tp <- seq(0, 1, 1 / ffs)
pp <- fnx(tp, 30)
fs <- 2e3
t <- seq(0, 1.2, 1 / fs)
d <- seq(0, 1, 1/3)
dd <- cbind(d, 4^-d)
z <- pulstran(t, dd, pp, ffs)
plot(t, z, type = "l", xlab = "", ylab = "",
     main = "Custom pulse")

par(op)


## ----findpeaks, fig.height=4, fig.width=7-------------------------------------

t <- 2 * pi * seq(0, 1,length = 1024)
y <- sin(3.14 * t) + 0.5 * cos(6.09 * t) +
     0.1 * sin(10.11 * t + 1 / 6) + 0.1 * sin(15.3 * t + 1 / 3)
data1 <- abs(y) # Positive values
peaks1 <- findpeaks(data1)
data2 <- y # Double-sided
peaks2 <- findpeaks(data2, DoubleSided = TRUE)
peaks3 <- findpeaks (data2, DoubleSided = TRUE, MinPeakHeight = 0.5)

op <- par(mfrow=c(1,2))
plot(t, data1, type="l", xlab="", ylab="")
points(t[peaks1$loc], peaks1$pks, col = "red", pch = 1)
plot(t, data2, type = "l", xlab = "", ylab = "")
points(t[peaks2$loc], peaks2$pks, col = "red", pch = 1)
points(t[peaks3$loc], peaks3$pks, col = "red", pch = 4)
par (op)
title("Finding the peaks of smooth data is not a big deal")

t <- 2 * pi * seq(0, 1, length = 1024)
y <- sin(3.14 * t) + 0.5 * cos(6.09 * t) + 0.1 *
     sin(10.11 * t + 1 / 6) + 0.1 * sin(15.3 * t + 1 / 3)
data <- abs(y + 0.1*rnorm(length(y),1))   # Positive values + noise
peaks1 <- findpeaks(data, MinPeakHeight=1)
dt <- t[2]-t[1]
peaks2 <- findpeaks(data, MinPeakHeight=1, MinPeakDistance=round(0.5/dt))
op <- par(mfrow=c(1,2))
plot(t, data, type="l", xlab="", ylab="")
points (t[peaks1$loc],peaks1$pks,col="red", pch=1)
plot(t, data, type="l", xlab="", ylab="")
points (t[peaks2$loc],peaks2$pks,col="red", pch=1)
par (op)
title(paste("Noisy data may need tuning of the parameters.\n",
            "In the 2nd example, MinPeakDistance is used\n",
            "as a smoother of the peaks"))


## ----fir, fig.height=4, fig.width=7-------------------------------------------

## FIR filter design by windowing

# low-pass filter 10 Hz
fs = 256
h <- fir1(40, 10/ (fs / 2), "low")
freqz(h, fs = fs)
# observe the effect of filter length
h <- fir1(80, 10/ (fs / 2), "low")
freqz(h, fs = fs)

# fir2 allows specifying arbitrary frequency responses
f <- c(0, 0.3, 0.3, 0.6, 0.6, 1)
m <- c(0, 0, 1, 1/2, 0, 0)
fh <- freqz(fir2(100, f, m))
op <- par(mfrow = c(1, 2))
plot(f, m, type = "b", ylab = "magnitude", xlab = "Frequency")
lines(fh$w / pi, abs(fh$h), col = "blue")
legend("topright", legend = c("specified", "actual"), lty = 1,
       pch = c(1, NA), col = c("black", "blue"))
# plot in dB:
plot(f, 20*log10(m+1e-5), type = "b", ylab = "dB", xlab = "Frequency")
lines(fh$w / pi, 20*log10(abs(fh$h)), col = "blue")
par(op)
title("specify arbitrary frequency responses with fir2")

## 50 Hz notch filter with remez
fs <- 200
nyquist <- fs / 2
f <- c(0, 48.5 / nyquist, 49.5 / nyquist, 50.5 / nyquist, 51.5 / nyquist, 1)
a <- c(1, 1, 0, 0, 1, 1)
h <- remez(200, f, a)
freqz(h, fs = fs)


## ----FIR_delay, fig.height=7, fig.width=7-------------------------------------
op <- par(mfrow = c(2, 1))
# design the filter
fs = 256
h <- fir1(40, 30/ (fs / 2), "low")

# group delay is constant at N/2
gd <- grpdelay(h)
plot(gd, ylim = c(0, 40),
     main = paste("(a) Group delay for FIR filters is constant\n",
                  "(here 40 / 2 = 20)"))

# filter electrocardiogram data with added noise
data(signals)
npts <- nrow(signals)
ecg <- signals$ecg + 1000 * runif(npts)
time <- seq(0, 10, length.out = npts)
plot(time, ecg, type = "l", main = "(b) Example ECG signal",
     xlab = "Time", ylab = "", xlim = c(0,2))
title(ylab = expression(paste("Amplitude (", mu, "V)")), line = 2)
f1 <- filter(h, ecg)
lines(time, f1, col = "red", lwd = 2)
delay <- mean(gd$gd)
f2 <- c(f1[(delay + 1):npts], rep(NA, delay + 1))
lines(time, f2, col = "blue", lwd = 2)
legend("topright", legend = c("Original", "Filtered", "Corrected"),
       lty = 1, lwd = c(1, 2, 2), col = c("black", "red", "blue"))

par(op)

## ----iir, fig.height=7, fig.width=7-------------------------------------------

op <- par(mfrow = c(3,1))

# compare Butterworth and Chebyshev filters.
bfr <- freqz(butter(5, 0.1))
cfr <- freqz(cheby1(5, .5, 0.1))
plot(bfr$w, 20 * log10(abs(bfr$h)), type = "l", ylim = c(-80, 0),
     xlab = "Frequency (rad)", ylab = c("dB"),
     main = "(a) Butterworth and Chebyshev")
lines(cfr$w, 20 * log10(abs(cfr$h)), col = "red")
legend("topright", legend = c("5th order Butterworth", "5th order Chebyshev"),
       lty = 1, col = c("black", "red"))

# compare Butterworth and elliptic filters.
efr <- freqz(ellip(5, 3, 40, 0.1))
plot(bfr$w, 20 * log10(abs(bfr$h)), type = "l", ylim = c(-80, 0),
    xlab = "Frequency (rad)", ylab = c("dB"),
    main = "(b) Butterworth and Elliptic")
lines(efr$w, 20 * log10(abs(efr$h)), col = "red")
legend ("topright", legend = c("5th order Butterworh", "5th order Elliptic"),
       lty = 1, col = c("black", "red"))

# compare type I and type II Chebyshev filters.
c1fr <- freqz(cheby1(5, .5, 0.1))
c2fr <- freqz(cheby2(5, 20, 0.1))
plot(c1fr$w, 20 * log10(abs(c1fr$h)), type = "l", ylim = c(-80, 0),
     xlab = "Frequency (rad)", ylab = c("dB"),
     main = "(c) Type I and II Chebyshev")
lines(c2fr$w, 20 * log10(abs(c2fr$h)), col = "red")
legend ("topright", legend = c("5th order Type I", "5th order Type II"),
       lty = 1, col = c("black", "red"))

par(op)

## ----zplane, fig.height=7, fig.width=7----------------------------------------

op <- par(no.readonly = TRUE)
n <- layout(matrix(c(1, 2, 3, 3), nrow = 2, byrow = TRUE))
stable <- butter(3, 0.2, "low", output = "Zpg")

# artificially adapt pole
instable <- stable
instable$p[2] <- instable$p[2] - 2

zplane(stable, main = "Stable")
zplane(instable, main = "Instable")

t <- seq(0, 1, len = 100)
x <- sin(2* pi * t * 2.3) + 0.5 * rnorm(length(t))
z1 <- filter(stable, x)
z2 <- filter(instable, x)
plot(t, x, type = "l", xlab = "", ylab = "")
lines(t, z1, col = "green", lwd = 2)
lines(t, z2, col = "red")
legend("bottomleft", legend = c("Original", "Stable", "Instable"),
       lty = 1, col = c("black", "green", "red"), ncol = 3)
par(op)


## ----IIR_delay, fig.height=7, fig.width=7-------------------------------------
op <- par(mfrow = c(2, 1))
ell <- ellip(5, 0.1, 60, 30/(fs/2), "low")
ellf <- freqz(ell, fs = fs)
argh <- Arg(ellf$h)
argh[which(is.na(argh))] <- 0
phase <- unwrap(argh)
plot(ellf$w, phase, type = "l", xlab = "Frequency (Hz)", ylab = "Phase",
     main = paste("30 Hz 5th order elliptical low-pass IIR filter\n",
                  "phase response is not linear"))
gd <- grpdelay(ell, fs = fs)
plot(gd, main = paste("group delay depends on frequency\n",
                      "mean:", round(mean(gd$gd), 1), "samples"))
par(op)


## ----IIR_filtfilt, fig.height=4, fig.width=7----------------------------------

f <- filter(ell, ecg)
ff <- filtfilt(ell, ecg)
plot(time, ecg, type = "l", xlab = "Time", ylab = "", xlim = c(0,2))
title(ylab = expression(paste("Amplitude (", mu, "V)")), line = 2)
lines(time, f, col = "red", lwd = 2)
lines(time, ff, col = "blue", lwd = 2)
legend("topright", legend = c("Original", "filter()", "filtfilt()"),
       lty = 1, lwd = c(1, 2, 2), col = c("black", "red", "blue"))

## ----conv1--------------------------------------------------------------------

u <- rep(1, 3)
v <- c(1, 1, 0, 0, 0, 1, 1)
conv(u, v, "full")
conv(u, v, "same")
conv(u, v, "valid")
conv(v, u, "valid")



## ----conv2, fig.height=7, fig.width=7-----------------------------------------

  short <- runif(20L)
  long <- runif(1000L)

  # convolve two long series
  ll <- microbenchmark::microbenchmark(conv(long, long),
                                       cconv(long, long),
                                       fftconv(long, long))
  plot1 <- ggplot2::autoplot(ll)

  # convolve a short and a long series
  sl <- microbenchmark::microbenchmark(conv(short, long),
                                       cconv(short, long),
                                       fftconv(short, long))
  plot2 <- ggplot2::autoplot(sl)
 
  gridExtra::grid.arrange(plot1, plot2, nrow = 2, ncol = 1)
  

## ----chunks-------------------------------------------------------------------

N <- 10000L
long <- runif(N)
part1 <- long[1:(N / 2)]
part2 <- long[((N / 2) + 1):N]

b <- c(2, 3)
a <- c(1, 0.2)

y <- filter(b, a, long)
y1 <- filter(b, a, part1, 'zf')
y2 <- filter(b, a, part2, y1$zf)
yy <- c(y1$y, y2$y)
all.equal(y, yy)



## ----zi, fig.height=5, fig.width=7--------------------------------------------

t <- seq(-1, 1, length.out =  201)
x <- (sin(2 * pi * 0.75 * t * (1 - t) + 2.1)
     + 0.1 * sin(2 * pi * 1.25 * t + 1)
     + 0.18 * cos(2 * pi * 3.85 * t))
h <- butter(3, 0.05)
zi <- filter_zi(h)
## alternatively, use:
## lab <- max(length(h$b), length(h$a)) - 1
## zi <- filtic(h, rep(1, lab), rep(1, lab))
z1 <- filter(h, x)
z2 <- filter(h, x, zi * x[1])
plot(t, x, type = "l", xlab ="", ylab = "")
lines(t, z1, col = "red")
lines(t, z2$y, col = "green")
legend("bottomright", legend = c("Original signal",
       "Filtered without initial conditions",
       "Filtered with initial conditions"),
      lty = 1, col = c("black", "red", "green"))



## ----fft_ar, fig.height=7, fig.width=7----------------------------------------

op <- par(mfrow = c(3, 1))
fs <- 200
nsecs <- 10
lx <- fs * nsecs
t <- seq(0, nsecs, length.out =  lx)
# signal of 5 Hz + 12 Hz + noise
x <- (sin(2 * pi * 5 * t)
     + sin(2 * pi * 12 * t)
     + runif(lx))
plot(t, x, type = "l", xlab = "Time (s)", ylab = "", main = "Original signal")
pw <- pwelch(x, window = lx, fs = fs, detrend = "none")
plot(pw, xlim = c(0, 20), main = "PSD estimate using FFT")

py <- pyulear(x, 30, fs = fs)
plot(py, xlim = c(0, 20), main = "PSD estimate using Yule-Walker")
par(op)


## ----welch, fig.height=7, fig.width=7-----------------------------------------

op <- par(mfrow = c(3, 1))
fs <- 200
nsecs <- 100
lx <- fs * nsecs
t <- seq(0, nsecs, length.out =  lx)
# sine and cosine of signal of 5 Hz noise
x1 <- cos(2 * pi * 5 * t) + runif(lx)
x2 <- sin(2 * pi * 5 * t) + runif(lx)
x <- cbind(x1, x2)
pw <- pwelch(x, fs = fs)
plot(pw, plot.type = "spectrum", yscale = "dB", xlim = c(0, 50),
    main = "A sine and a cosine of 5 Hz have the same PSD")
legend("topright", legend = c("Cosine", "Sine"), lty = 1:2, col = 1:2)
rect(3, -35, 7, -4, border = "red", lwd = 3)
plot(pw, plot.type = "phase", xlim = c(0, 50),
    main = expression(bold(paste("but differ ", pi/2, " radians in phase at 5 Hz"))))
rect(3, -pi, 7, pi, border = "red", lwd = 3)
plot(pw, plot.type = "coherence", xlim = c(0, 50),
    main = "leading to coherence ~ 1 at 5 Hz")
rect(3, 0, 7, 1, border = "red", lwd = 3)
par(op)


## ----specgram, fig.height=7, fig.width=7--------------------------------------

op <- par(mfrow = c(2, 1))

jet <- grDevices::colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F",
                                     "yellow", "#FF7F00", "red", "#7F0000"))
sp <- specgram(chirp(seq(0, 5, by = 1/8000), 200, 2, 500, "logarithmic"), fs = 8000)
plot(sp, col = jet(20))

c2w <- grDevices::colorRampPalette(colors = c("red", "white", "blue"))
plot(sp, col = c2w(50))

par(op)