This article serves the purpose of illustrating that signal processing with R is possible – thanks to the signal package – and to keep a reference of some of the stuff that I learned at my last edX course. Anyway, I am by no means an expert on signal processing so I’d prefer to let the pictures and the code speak for themselves. But to give you the idea – I show case the creation and application of an FIR band pass filter (Chebyshev Type 1 in this case) and of an FIR filter created using the Parks-McClellan method with the Remez exchange algorithm. The code snippets are taken from a larger R script which you can find on GitHub. I aim to focus on the essential parts. You’re welcome to share your knowledge and corrections by leaving a comment.

(The expression 50 in 1000 or 50/1000 is supposed to refer to 50 repetitions / periods within a sample of length 1000.)

Chebyshev Type 1 Band Pass Filter

The signal is an additive combination of four sinusoids with frequencies 300, 500, 700 and 1100 in 10000. The goal is to filter out all components except for 500/10000.

# length of signal
n <- 10000

# the frequency to be kept
F <- 500
F0 <- 2 * F / n

# input signal
sig <- add_sin_sig(n, c(300,500,700,1100))

# filter specification
filter <- cheby1(6,2,c(F0-F0*.1,F0+F0*.1),type="pass")

# filtered signal
sig2 <- signal::filter(filter, x=sig)

FIR Band Stop Filter with Parks-McClellan Method

In this case I use a single sinus function whose frequency increases linearly from 1 to 10’000 in 100’000. The band stop filter I apply here is created using Parks-McClellan method which is provided by the signal package implementing the Remez exchange algorithm.

pm_filter <- function(n_sig, n_fil, freq, type, d1=.05, d2=.07) {
  c0 <- 2 * freq / n_sig
  
  if(type == "stop") {
    v <- c(1,1,0,0,1,1)
  } else if(type == "pass") {
    v <- c(0,0,1,1,0,0)
  }
  
  fir <- remez(n = n_fil,f=c(0,c0-d2,c0-d1,c0+d1,c0+d2,1),a = v)
  freq <- freqz(fir, n=n_fil)
  return(list(fir = fir, freq = freq))
}

apply_filter_to_sig <- function(fir, sig) {
  y <- signal::filter(as.vector(fir), 1, x = sig)
}

# signal length
n <- 100000

# filter length
n_fil <- 800

# the frequency to be filtered out
F <- 5000
F0 <- 2 * F / n

sig <- sin_sig_incr_freq(n, from = 0, to = 10000)

filter <- pm_filter(
  n_sig = n, n_fil = n_fil, freq = F, type = "stop", d1 = F0*.05, d2 = F0*.1)

sig2 <- apply_filter_to_sig(filter$fir, sig)

Approximating the z-Plane for the Band Stop Filter

I am not sure if it is possible to calculate the zeros and poles algebraically. So I simply used a numerical method. As the absolute values practically explode towards infinity towards the center and implode towards zero beyond the unit circle I combined two  scales – one for lower than 1 and one for larger than 1. Note that I take the 10 Mio.th logarithm.

G <- expand.grid(a=-100:100/100,b=-100:100/100)

# http://en.wikipedia.org/wiki/Z-transform
for(i in 1:nrow(G)){ 
  G[i,"v"] <- abs(sum(filter$freq$h*((G[i,"a"]+G[i,"b"]*1i)^-(0:799))))
}

levelplot(
  log(matrix(G[,"v"], ncol=201, byrow=TRUE))/log(10000000), 
  col.regions=colorRampPalette(
    c("blue","red", "green"), space = "Lab")(100),
  at=(
    c(seq(-.65,0,length.out=50)^7,
      seq(0,70,length.out=51)[2:51]))
)

(original article published on www.joyofdata.de)