>>>>> "tom" == tom wright <[EMAIL PROTECTED]> >>>>> 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 https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html