[R] Filter design in R?

[tom wright]

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

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

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


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

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

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

    filt.total<-sum(as.numeric(f))
    f<-as.numeric(f)/filt.total
}

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

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

    t<-c(0:length+1)

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

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

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

    filt.total<-sum(as.numeric(f.low))
    f.low<-as.numeric(f.low)/filt.total

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

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

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

    filt.total<-sum(as.numeric(f.high))
    f.high<-as.numeric(f.high)/filt.total

#invert the high filter to make it high pass

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

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

#make band pass by spectral inversion
    f.bandpass<-0-f.bandreject
    f.bandpass[length/2]<-f.bandpass[length/2]+1
    f.bandpass
}
==============END CODE=================