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
-
- 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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
# 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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 |
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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
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)