From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Fri 21 Oct 2005 - 19:33:03 EST

>>>>> "tom" == tom wright <tom@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=================

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