#This file contains R code to read the ecg.csv #file from my website, perform some digital #signal processing, and produce a nice plot #of the raw and filtered signals. This script #may be run from R with the follwing command #source("http://biostatmatt.com/R/ecg.R") #Transform the real and imaginary portions of the #FFT into magnitude and phase. The argument #ff should be the output of the fft function amplitude <- function( x ) { sqrt(Re(x)^2+Im(x)^2) } phase <- function( x ) { atan(Im(x)/Re(x)) } #sinc function of frequency f sinc <- function( x, f ) { ifelse(x==0, 2*pi*f, sin(2*pi*f*x)/x) } #Blackman window from 0..m Blackman <- function( m ) { 0.42-0.5*cos(2*pi*(0:m)/m)+0.08*cos(4*pi*(0:m)/m) } #Hamming window from 0..m Hamming <- function( m ) { 0.54-0.46*cos(2*pi*(0:m)/m) } #simple low pass filter #y - vector to filter #t - time interval between measurements (s) #f - low pass frequency (Hz) lpf <- function( y, t, f ) { rc <- 1 / ( 2 * pi * f ) a <- t / ( t + rc ) n <- length( y ) yf <- y for( i in 2:length(y) ) { yf[i] <- a * y[i] + (1-a) * yf[i-1] } return( yf ) } #windowed sinc low pass filter #y - vector to filter #t - time interval between measurements (s) #f - low pass frequency (Hz) wlpf <- function( y, t, f ) { m <- min(floor(length(y)/2), 500) #generate the sinc kernel rk <- sinc(-m:m, f*t) #apply the Blackman window bk <- Blackman(2*m) * rk #pad the filter with zeros k <- c(bk, rep(0,length(y)-length(bk))) #convolve y with the filter kernel fy <- fft(fft(k)*fft(y), inverse=TRUE) return(Re(fy)) } dat <- scan("http://biostatmatt.com/csv/ecg.csv") dat <- ( dat - mean(dat) ) / sd(dat) #filter high frequency noise fdat <- wlpf(dat, 1/1000, 30) fdat <- (fdat-mean(fdat))/sd(fdat) #isolate respiration bias rdat <- wlpf(dat, 1/1000, 1) rdat <- (rdat-mean(rdat))/sd(rdat) #subtract respiration bias from ecg signal edat <- fdat - rdat require(lattice) xplot <- rep((0:(length(dat)-1))/1000,4) yplot <- c(dat, fdat, rdat, edat) gplot <- c(rep("Raw",length(dat)), rep("High Frequency Filter",length(dat)), rep("Low Frequency Filter",length(dat)), rep("ECG",length(dat))) tp <- xyplot(yplot~xplot|gplot,type="l",layout=c(1,4), xlab="Time", ylab="V") #uncomment the following to save an image #trellis.device(png, file="ecgfilter.png", height=750, width=750) print(tp) #dev.off()