TITLE filtered white noise membrane channel with IIR filter

COMMENT
This mod file produces white noise current, which is filtered through
a IIR filter. The Filter coeeficient are provided in the array 'a' and 'b'

  i(n) = b[0]*z(t) + b[1]*z(t-dt) + ... + b[nb]*z(t-nb*dt)
                   - a[1]*i(t-dt) - ... - a[na]*i(t-na*dt)
 
Where na = nb = bufsz. Note that a[0] is assumed to be one,
and skipped in the summation.

The random generator is implemnted in RNG.mod and uses is gasdev of numerical rec. 
using the rejection method to produce gaussian RV given a uniform distribution. 
This is emulated using ran1.

Important Notes: 
1. The variable dummy is added such that Neuron will not use VECTORIZED
   method, in this way we can use VERBATIM functions. (Otherwise some of 
   the variables get random values).

2. The variable Am is a scaling factor, in order to match the desired Power 
   Spectral Density of the current. Since the noise is independet for each channel 
   the variance is scaled linearly with the area of the patch, and thus we
   need to divide by the sqrt of the area, to get the true scaling factor (SF)
   The units doesn't seems to make sense but it is o.k.
ENDCOMMENT


UNITS {
	(mV) = (millivolt)
	(mA) = (milliamp)
}

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}


DEFINE SIZE 5

NEURON {
	SUFFIX iirsyn
	NONSPECIFIC_CURRENT i
	RANGE f,a,b,n,o,bufsz,Am,SF,ptr
	GLOBAL dummy
}

ASSIGNED { 
     dummy
     i (mA/cm2)
     f (mA/cm2)
     a[SIZE]          : output  coefficient
     b[SIZE]          : input coefficient
     o[SIZE] (mA/cm2) : output buffer
     n[SIZE] (mA/cm2) : noise  buffer
}

PARAMETER {
	bufsz = SIZE
        Am   = 1e-4 (mA/cm) 
        SF   = 1.0  (mA/cm2)
        ptr = 0

}

INITIAL {
  SF = Am/sqrt(area)
  fillbuf()    
}

BREAKPOINT {
        SOLVE filter
        i = -f         
}

PROCEDURE fillbuf()
{
LOCAL j
VERBATIM
  int jk;
  extern double gasdev();

  for (jk = 0;jk<(int)bufsz;jk++){
      n[jk] = gasdev();
      o[jk] = gasdev();
  }
  ENDVERBATIM
  j = 0
  WHILE (j<bufsz*100) {
      filter()
      j = j + 1
  }
  VERBATIM
    return 0;
  ENDVERBATIM
}
 
PROCEDURE filter()
{
VERBATIM
   int pr;                      /* index pointer to circular arrays */
   int kj,ki;                               
   extern double gasdev();

   pr = (int)ptr;
   n[pr] = gasdev();
   f = b[0]*n[pr];
   for (kj =1;kj<(int)bufsz;kj++) {        /* loop through coefficients */
      ki = ((int)pr+kj)%(int)bufsz;        /* index to circular array */
      f += b[kj]*n[ki] - a[kj]*o[ki];      /* convolve the noise & former output with filter */
   }  
   o[pr] = f;                              /* store new output in array */
   f = f*SF;                               /* scale the current to have the right */
                                           /* amplitude in the approp units. */
   pr = ((int)bufsz+(pr-1))%((int)bufsz);  /* update position pointer  */
   ptr = pr;

ENDVERBATIM
}