Re: [R] Filter design in R?

From: Martin Maechler <>
Date: Fri 21 Oct 2005 - 19:33:03 EST

>>>>> "tom" == tom wright <> >>>>> on Thu, 20 Oct 2005 09:26:05 -0400 writes:

    tom> On Wed, 2005-19-10 at 17:49 -0400, Israel Christie wrote:
>> Dr. Williams,
>> I ran across your inquiry on one of the R-help mailing lists regarding
>> digital filter design and implementation. I found no response to your
>> email in the archives and was wondering if you were able to find anything.
>> Thanks,
>> Israel

    tom> I'm not Dr Williams, but I've been doing some work on filter design
    tom> recently. I'm also no expert in this area but I found some very useful
    tom> resources, primarily the online book "The scientist and engineers guide
    tom> to digital signal processing"

    tom> I came up with some code for generating simple highpass; low pass and     tom> bandpass filters, the filters can be applied using the filter() function

    tom> Since I'm no expert here I'd really appreciate any comment from people     tom> who know more than me about these techniques.

I'm neither Dr. Williams nor an expert in the ("engineer point of view on") filter design.

Note however that R (in 'stats') , additionally to filter() has a function kernel() which is used to build some common (non-recursive) filters and produce objects of (S3) class "tskernel".
These were primarily designed for use in spectrum(), but should be more generally useful, and could serve as an initial role model for extention.

The "real" source is in

Martin Maechler

    tom> Regards
    tom> Tom

>> ================== BEGIN USAGE CODE==================
>> t<-c(1:1000)/1000 #timeline 1KHz
>> s1<-sin(2*pi*t*3) #3Hz waveform
>> s2<-sin(2*pi*t*5) #5Hz waveform
>> s3<-sin(2*pi*t*10) #10Hz waveform
>> stot<-s1+s2+s3 #complex waveform
>> plot(stot,type='l')
>> #create the filter, the longer it is the better cutoff
>> #length must be an even number
>> f<-calcbpfilt(length=900,samplerate=1000,lowfreq=7,highfreq=4)
>> sfilt<-filter(stot,f,circular=TRUE) #apply the filter
>> lines(sfilt,type='l',col='red')
>> #only the 5Hz freq should be let through
>> ================== END USAGE CODE==================

    tom> ==============BEGIN CODE=================
    tom> calclpfilt<-function(length,fc){     
    tom> t<-c(0:length+1)
    tom> for (val in t){
    tom> if(val-length/2==0){
    tom> f[val]<-as.numeric(2*pi*fc)
    tom> }else{

    tom> f[val]<-as.numeric(sin(2*pi*fc*(val-length/2))/(val-length/2))
    tom> }
    tom> f[val]=as.numeric(f[val])*(0.54-0.46*cos(2*pi*val/length))
    tom> }

    tom> #f<-convolve(f,f)
    tom> #normalise filter
    tom> f<-as.numeric(f)/
    tom> }

    tom> calcbpfilt<-function(length,samplerate,lowfreq,highfreq){
    tom> f.low<-list()
    tom> f.high<-list()

    tom> fc.low<-1/(samplerate/lowfreq)
    tom> fc.high<-1/(samplerate/highfreq)

    tom> t<-c(0:length+1)

    tom> #calculate the lowpass filter
    tom> for (val in t){
    tom> if(val-length/2==0){
    tom> f.low[val]<-as.numeric(2*pi*fc.low)
    tom> }else{

    tom> f.low[val]<-as.numeric(sin(2*pi*fc.low*(val-length/2))/(val-length/2))     tom> }

    tom> f.low[val]=as.numeric(f.low[val])*(0.54-0.46*cos(2*pi*val/length))
    tom> }
    tom> #f<-convolve(f,f)
    tom> #normalise filter

    tom><-sum(as.numeric(f.low))     tom> f.low<-as.numeric(f.low)/

    tom> #calculate the second filter
    tom> for (val in t){
    tom> if(val-length/2==0){
    tom> f.high[val]<-as.numeric(2*pi*fc.high)
    tom> }else{

    tom> f.high[val]<-as.numeric(sin(2*pi*fc.high*(val-length/2))/(val-length/2))     tom> }

    tom> f.high[val]=as.numeric(f.high[val])*(0.54-0.46*cos(2*pi*val/length))
    tom> }
    tom> #f<-convolve(f,f)
    tom> #normalise filter

    tom><-sum(as.numeric(f.high))     tom> f.high<-as.numeric(f.high)/

    tom> #invert the high filter to make it high pass

    tom> f.high<-0-f.high
    tom> f.high[length/2]<-f.high[length/2]+1

    tom> #add lowpass filterkernal and highpass filter kernel
    tom> #makes band reject filter
    tom> f.bandreject<-f.low+f.high

    tom> #make band pass by spectral inversion
    tom> f.bandpass<-0-f.bandreject
    tom> f.bandpass[length/2]<-f.bandpass[length/2]+1
    tom> f.bandpass

    tom> }
    tom> ==============END CODE=================

______________________________________________ mailing list PLEASE do read the posting guide! Received on Fri Oct 21 19:39:56 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 19:15:55 EST