FIR Filter Design and Digital Signal Processing in R

May 15, 2014
By

[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
z-Plane

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)

cheb1-signal

FIR Band Stop Filter with Parks-McClellan Method

remez-imp-freq-res
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-signal
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 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)

To leave a comment for the author, please follow the link and comment on their blog: joy of data » R.

R-bloggers.com 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.



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

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)