// $Id: filtutils.hoc,v 1.6 2010/10/11 18:34:24 samn Exp $ //* mkgauss(vector,average,standard-dev) proc mkgauss () { local i,x,sz,mu,sd localobj vin vin=$o1 sz=vin.size mu=$2 sd=$3 for vtr(&x,vin,&i) vin.x(i) = exp( -(x-mu)^2 / (2*sd*sd) ) vin.mul( 1 / (sd*sqrt(2*PI)) ) } //* mktriangwin(vec,size - should be odd,[skip the wraparound]) proc mktriangwin () { local i,j,sz localobj vin vin=$o1 vin.resize($2) vin.x(int($2/2))=1 j=1 sz=1/(vin.size/2) for (i=int($2/2)-1;i>=0;i-=1) { vin.x(i)=j j-=sz } j=1 for i=int($2/2)+1,vin.size-1 { vin.x(i)=j j-=sz } vin.div(vin.sum) if(numarg()>2) return vin.wraparound(vin.size) } //* mkgaussfilt(vec,stdev[,vx]) proc mkgaussfilt () { local sd,minx,maxx,dx,a localobj vin,vx vin=$o1 sd=$2 if(numarg()>2) {vx=$o3 vin.resize(0) vin.copy(vx)} mkgauss(vin,0,sd) vin.wraparound(vin.size) vin.div(vin.sum) dealloc(a) } //* dofilt(vsignal,vwindow) - filters with convlv proc dofilt () { local a,i localobj vsig,vwin,vtmp a=allocvecs(vtmp) vsig=$o1 vwin=$o2 sz=vsig.size vtmp.convlv(vsig,vwin) vsig.copy(vtmp) vsig.resize(sz) // make sure size doesn't change dealloc(a) } //* triangfilt(vin,filtsize) - run a triangle filter proc triangfilt () { local a localobj vin,vwin vin=$o1 a=allocvecs(vwin) mktriangwin(vwin,$2) // make the window dofilt(vin,vwin) // do the filtering dealloc(a) } //* boxfilt(vin,filtsize) - run a box(moving average) filter proc boxfilt () { local a localobj vin,vwin vin=$o1 a=allocvecs(vwin) {vwin.resize($2) vwin.fill(1) vwin.div(vwin.size)} // make the window dofilt(vin,vwin) // do the filtering dealloc(a) } //* gaussfilt(vin,stdev,vx) - run a gaussian filter - vx is x-values used to make gaussian proc gaussfilt () { local a,sd localobj vin,vwin,vx vin=$o1 sd=$2 vx=$o3 a=allocvecs(vwin) mkgaussfilt(vwin,sd,vx) // make the window dofilt(vin,vwin) // do the filtering dealloc(a) } //* myfilt(code,vec) - code:0=gauss,1=triangle,2=box proc myfilt () { local a localobj vx if($1==0) { a=allocvecs(vx) vx.indgen(-3,3,.03) gaussfilt($o2,stdg,vx) } else if($1==1) { triangfilt($o2,winsz) } else if($1==2) { boxfilt($o2,winsz) } } //* resample(vec,new size) - resample a vec to new size using linear interpolation proc resample(){ local newsz,idxdest,idxsrc,val,fctr,frac localobj vtmp {vtmp=new Vector($2) fctr=$o1.size/$2 vtmp.x(0)=$o1.x(0) idxsrc=fctr} for(idxdest=1;idxdest<$2-1;idxdest+=1){ idxsrc = idxdest * fctr frac = idxsrc - int(idxsrc) idxsrc = int(idxsrc) if(idxsrc+1>=$o1.size){ vtmp.x(idxdest) = $o1.x(idxsrc) continue } val = (1-frac) * $o1.x(idxsrc) + frac * $o1.x(idxsrc+1) vtmp.x(idxdest) = val } {vtmp.x($2-1)=$o1.x($o1.size-1) $o1.resize($2) $o1.copy(vtmp)} }