[R] Filter design in R?

[Martin Maechler]

>>>>> "tom" == tom wright <tom at maladmin.com>
>>>>>     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" http://www.dspguide.com

    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
https://svn.R-project.org/R/trunk/src/library/stats/R/kernel.R

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> filt.total<-sum(as.numeric(f))
    tom> f<-as.numeric(f)/filt.total
    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> filt.total<-sum(as.numeric(f.low))
    tom> f.low<-as.numeric(f.low)/filt.total

    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> filt.total<-sum(as.numeric(f.high))
    tom> f.high<-as.numeric(f.high)/filt.total

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