FIR Filter Design and Digital Signal Processing in R

[This article was first published on joy of data » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

iconThis 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

Impulse Response
Frequency Response

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)

remez-zplaneApproximating 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 logarithm.

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

for(i in 1:nrow(G)){ 
  G[i,"v"] <- abs(sum(filter$freq$h*((G[i,"a"]+G[i,"b"]*1i)^-(0:799))))

  log(matrix(G[,"v"], ncol=201, byrow=TRUE))/log(10000000), 
    c("blue","red", "green"), space = "Lab")(100),

(original article published on

To leave a comment for the author, please follow the link and comment on their blog: joy of data » R. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)