Date: Thu 20 Oct 2005 - 23:26:05 EST

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

______________________________________________
```
R-help@stat.math.ethz.ch mailing list