fftilt

Accession:62685
Low-pass filtering using NEURON's fft() routine
Tool Information (Click on a link to find other Tools with that property)
Tool Type: Analysis; Convolution; Utility;
Simulation Environment: NEURON;
\
fftilt
readme.txt
hanning.gif
blackman.gif
reflected.gif
test.gif
fftilt.hoc
                            
// $Id: fftilt.hoc,v 1.3 2006/01/21 15:20:56 billl Exp $

objref vec,vec0,g[10]
vec=new Vector()
vec0=vec.c
g=new Graph()
pstep=0.05 // time step (dt) of the sample

//* fftilt(src,dest[,cutoff])
// use fft to filter noise out of a vector
obfunc fftilt () { local sz,sz2,min,max,loc,cutoff localobj v1,v2
  v1=new Vector() v2=v1.c
  v1.copy($o1)
  min=v1.min max=v1.max
  if (numarg()==3) cutoff=$3 else cutoff=50
  loc=round(cutoff/2/PI/pstep) // translate freq to location on 
  sz=v1.size
  sz2=2^round(log(sz)/log(2))      // nearest power of 2 length
  // sz2=2^int((log(sz)/log(2))+1) // alternative: alway enlarge
  v1.resize(sz2)
  if (sz2>sz) v1.fill(0,sz,sz2-1) // pad with zeros
  v2.copy(v1)
  v1.reverse()
  v1.append(v2) // reflection across y-axis places beginning of trace in center
  v1.mul(hanning(v1.size)) // size will now be sz2*2
  v2.fft(v1)
  v2.mul(blackman(v1.size,loc))
  v1.fft(v2,-1)
  v1.resize(v1.size/2,v1.size-1) // get rid of the reflection
  v1.reverse
  if (sz2>sz) v1.resize(sz) // get rid of any zero padding
  v1.scale(min,max)  // restore the scale
  $o2.copy(v1)
  return $o2
}

// ffilt_test(cutoff[,blackman_size]) 
// mix freq's of 10 and 50 Hz; filter out above 'cutoff'; display
// eg ffilt_test(20,32)
proc ffilt_test () { local f,coff
  coff=$1
  if (numarg()==2) blackman_size=$2
  pstep=0.05
  vec.resize(15000)
  vec.fill(0)
  for case(&f,10,50) vec.add(vec.c.sin(f,0,pstep))
  fftilt(vec,vec0,coff) // cutoff at $1
  g.erase_all()
  vec.line(g,pstep,1,2) vec0.line(g,pstep,2,2)
  g.exec_menu("View = plot")
}

//* hanning(vec_size)
// Hanning window drops off to zero at edges so as to avoid discontinuity artifact
// note that this means that beginning and end of signal are unrepresented
// (see 'reflection' note in fftilt())
nsz=0
func hann () { return 0.5-(0.5*cos(2*PI*$1/nsz)) }
obfunc hanning () { localobj v1
  nsz=$1
  v1=new Vector(nsz)
  v1.indgen(1)
  v1.apply("hann")
  return v1
}

//* blackman(vec_size,fall_off_loc)
// Blackman window gives gradual fall-off so don't get ringing from sudden cutoff of fft
// eg aa=blackman(16384,512,30)
blackman_size=128
func blkmn () { 
  return 1-0.42+(0.5*cos(2*PI*$1/blackman_size/2))-(.08*cos(4*PI*$1/blackman_size/2)) 
}
obfunc blackman () { local a,sz,loc localobj v1,vr
  sz=$1 loc=$2
  if (blackman_size>loc) { 
    printf("blackman() ERRA: window of blackman_size %d won't fit before cutoff at index %d\n",\
           blackman_size,loc)
    err() } // generate an error break by calling nonexistent function
  if (loc>sz) { 
    printf("blackman() ERRB: cutoff too high; index %d > sample size %d\n",loc,sz) err() }
  vr=new Vector(blackman_size) v1=vr.c
  vr.resize(loc-blackman_size) vr.fill(1)
  v1.resize(blackman_size) v1.indgen(1)
  v1.apply("blkmn")
  vr.append(v1)
  v1.resize(sz-vr.size) v1.fill(0)
  vr.append(v1)  
  return vr
}

// round() round off to nearest integer
func round () { local ii
  if (argtype(1)==1) {
    if ($o1.size==0) return 1e9
    for ii=0,$o1.size-1 {
      if ($o1.x[ii]>0) $o1.x[ii]=int($o1.x[ii]+0.5) else $o1.x[ii]=int($o1.x[ii]-0.5) 
    }
    return($o1.x[0])
  } else {
    if ($1>0) return int($1+0.5) else return int($1-0.5) 
  }
} 

ffilt_test(20,32)