Filter design in R?

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=================