: $Id: place.mod,v 1.141 2011/07/06 15:24:40 samn Exp $
:* COMMENT
COMMENT
CODE modified from ~/nrniv/place/fenton/
SimulateTS/MakeGaussianField.c,SimulateTS.c
CorrelateSpikeTrains/CorrelateSpikeTrains.c,Kendall.c
ENDCOMMENT
NEURON {
SUFFIX PLACE
GLOBAL INSTALLED, MAKE_GAUSSIAN, REFRAC
GLOBAL DEBUG, PEAK, THRESHOLD
GLOBAL samps_per_sec,scale_y,first_time_stamp,scale_x
GLOBAL MaxTime,AnalysisInterval,BinSize,IntervalStart,VERBOSE
}
PARAMETER {
INSTALLED=0
VERBOSE=0
DEBUG=0
MAKE_GAUSSIAN=1
REFRAC=3
PEAK = 5.0
THRESHOLD= 1e-5
samps_per_sec=0
scale_y=0
first_time_stamp=0
scale_x=0
MaxTime=1200 : seconds
AnalysisInterval=1200 : seconds
IntervalStart=0 : seconds
BinSize=250
}
ASSIGNED { }
VERBATIM
#include "bpf.h"
#include "misc.h"
void KendallTauPairs(si1 *dx, si1 *dy, ui4 NumEvents, ui4 NumPairs, sf8 *tau, sf8 *z, sf8 *prob);
static float **rate_map;
static double storage[10];
#define MAP_YSZ 32
#define MAP_XSZ 32
#define YSZ 256
#define XSZ 256
#define MAX_RATE 100.0
#define MAX_SPKS_PER_SAMP 100
#define TIMESTAMP_UNITS_PER_SEC 10000
#define TRUE 1
#define FALSE 0
#define PARTS_PER_SAMP 2
//* from MakeGaussianField.c
static double gaussaf (mu, sd, x)
double mu, sd, x;
{ double exponent, num, denom, z;
double value, max_value;
z = (x - mu) / sd;
exponent = -0.5 * z * z;
max_value = (1.0 / (sd * (double)SQRT2PI));
value = (exp(exponent) / (sd * (double)SQRT2PI));
value /= max_value;
return(value);
}
// MakeGaussianField(int x, int y, double sd, int xsz, int ysz, float peak, float
// threshold, float **map)
static void MakeGaussianField (x, y, sd, xsz, ysz, peak, threshold, map)
int x, y, xsz, ysz;
double sd;
float peak, threshold, **map;
{ int i, j;
double hypot(), d;
for(i = 0; i < ysz; i++){
for(j = 0; j < xsz; j++){
d = hypot((double)(i - y), (double)(j - x));
map[i][j] = peak * gaussaf(0.0, sd, d);
if(map[i][j] < threshold) map[i][j] = 0.0;
if (DEBUG==2 && map[i][j]!=0) printf("MAP: %d %d %g\n", i,j,map[i][j]);
}
}
return;
}
//* From generate_spks.c
static int generate_spks(rate, max_rate, samps_per_sec)
float rate, max_rate;
double samps_per_sec;
{ register int i;
register double prob;
double drand48();
int reduction_factor = 1, parts_per_samp = PARTS_PER_SAMP;
int n_spks;
while(rate >= max_rate){
rate /= 2.0;
reduction_factor *= 2;
}
prob = (double)rate / (double)samps_per_sec / (double)parts_per_samp;
n_spks = 0;
for(i = parts_per_samp * reduction_factor; i--; )
if(prob > (-1.0 * log(drand48())))
n_spks++;
return n_spks;
}
//* original SimulateTC.c
// why not just use fgets?
static char *get_line(FILE* f,char* s,int iMaxLen)
{ int i=0;
while(fscanf(f,"%c",&s[i]) != EOF && i<iMaxLen-1){
if(s[i] == '\n'){
s[i+1] = '\0';
return s;
}
i++;
}
return(NULL);
}
static void instruct()
{
printf("Call with a time_series format file in TS_DIR.\n");
printf("Output to STS_DIR.\n");
printf("Simulates a location-specific spike train using the track in the TS file. Positional rates are derived from a 2-D Gaussian rate map\n");
printf("or a rate map stored as a file. See options.\n");
printf("Options:\n");
printf("\tF <filename> Don't use a Gaussian field for the simulation. Use the rate map stored in RATE_DIR/filename.\n");
printf("\t\tRequires options X and Y\n");
printf("\tP <cycles per minute> Modulate rate with a periodic function. Give the period (cycles per min).\n");
printf("\tS <Seed the random number generator with your own integer seed. Default is a random seed.\n");
printf("\tX <x resolution of rate map> The x dimension resolution of the rate map array. Necessary for option F. Default: %d\n", MAP_XSZ);
printf("\tY <y resolution of rate map> The y dimension resolution of the rate map array. Necessary for option F. Default: %d\n", MAP_YSZ);
hxe();
}
static int ascending_int(unsigned int *a, unsigned int *b)
{ if (*a > *b)
return(1);
if (*a < *b)
return(-1);
return(0);
}
ENDVERBATIM
:* v1.glob(L_BOUND,R_BOUND,cut) -- group values together and give bounds (indices)
VERBATIM
static double glob (void* vv) {
int i, j, nx, ny, nz;
double *x, *y, *z, last, cut, diff;
nx = vector_instance_px(vv, &x);
ny = vector_arg_px(1, &y);
nz = vector_arg_px(2, &z);
if (ny!=nz || ny<2) {printf("Place:glob ERR0: ny!=nz %d %d\n",ny,nz); hxe();}
cut=*getarg(3);
for (last=x[0],i=1,j=0,y[0]=0; i<nx; i++) {
diff=x[i]-last;
if (diff<0) {printf("Place:glob ERRA: not monotonic %g < gd\n",x[i],last); hxe();}
if (diff>cut) {
z[j]=(double)i-1;
j++;
if (j>=ny) { printf("Place:glob ERRB: OOR %d\n",ny); hxe(); }
y[j]=(double)i;
}
last=x[i];
}
z[j]=nx-1; // index at end of vector
vector_resize(vector_arg(1), j+1);
vector_resize(vector_arg(2), j+1);
return (double)j+1;
}
ENDVERBATIM
:* v1.cumul(out)
VERBATIM
static double cumul (void* vv) {
int i, j, nx, ny, nyn, fli;
unsigned int* scr;
double *x, *y, cut,fl;
void *vvy;
nx = vector_instance_px(vv, &x);
ny = vector_arg_px(1, &y); vvy=vector_arg(1);
cut=*getarg(2);
if (nx<2) { printf("cumul:ERRA nx too small %d\n",nx); hxe();}
vector_resize(vvy,ny=vector_buffer_size(vvy));
scr=scrset(nx);
for (i=0;i<nx;i++) scr[i]=i;
nrn_mlh_gsort(x, scr, nx, cmpdfn);
fl=floor(x[scr[nx-1]]); // this is the largest value
if ((fli=(int)fl)>=ny-1) { printf("cumul:ERRB fli>=ny %d %d\n",fli+1,ny); hxe();}
if (fl>=cut) {
nyn=fli+2;
for (i=cut;i<nyn-1;i++) y[i]=1;
for (i=0;i<cut;i++) y[i]=0;
y[fli+1]=0;
for (i=nx-2;i>=0 && x[scr[i]]>=cut;i--) {
if (x[scr[i]]<fl) fli=(int)(fl=floor(x[scr[i]]));
for (j=cut;j<=fli;j++) y[j]++;
}
} else nyn=0;
vector_resize(vvy, nyn);
return (double)nyn;
}
ENDVERBATIM
:* main from SimulateTS.c
VERBATIM
static double simts (void* vv) {
FILE *fopen(), *ifp, *ofp, *ofq, *rfp;
Object* ob;
unsigned char xys[3];
int i,j,k,s,mt,fflag;
int op, AngularModulation;
int x, y, nmb_header_lines;
double center_x, center_y, sd, value;
double time;
double max_rate;
char line[256], keyword[64], *get_line(), rate_file[256];
char in_file[256], out_file[256], file_name[256];
double modulate, peak, threshold;
double CyclesPerMin;
double map_x, map_y, collapse_x, collapse_y;
void instruct(), MakeGaussianField();
int found, n_samps;
double sample_time, timestamp_units_per_samp;
int SampsPerCycle, SampsPerHalfCycle, samps;
double *xx, *yy, jit, *loc[2];
int nx,ny,maxsz,ii, q, sflag, nloc[2];
void *vv1, *vloc[2];
nx = vector_instance_px(vv, &xx);
maxsz=vector_buffer_size(vv);
vector_resize(vv, maxsz);
sflag=0; jit=0.;
ob = *hoc_objgetarg(1);
if (strncmp(hoc_object_name(ob),"List",4)==0) {
fflag=0;
// pick up list of 2 vectors
for (i=0;i<2;i++) nloc[i]=list_vector_px2(ob, i, &loc[i], &vloc[i]);
if (nloc[0]!=nloc[1]) {printf("simtst_PLACE ERRA: %d %d\n",nloc[0],nloc[1]); hxe();}
} else {
fflag=1;
ifp=hoc_obj_file_arg(1);
}
if (ifarg(2)) {
if (hoc_is_double_arg(2)) {
jit= *getarg(2)*10.; // convert from ms to 100 mus
} else { // vector for storage
ny = vector_arg_px(2, &yy);
vv1=vector_arg(2);
if ((ii=vector_buffer_size(vv1))<maxsz) {
fprintf(stderr, "2nd vec maxsz (%d) should be at least %d\n",ii,maxsz); hxe(); }
vector_resize(vv1, maxsz);
sflag=1;
}
}
q=found=0;
AngularModulation = FALSE;
max_rate = MAX_RATE;
map_x = MAP_XSZ; map_y = MAP_YSZ;
while (0) switch (op) {
case 'F': ;// Gaussian_flag -> replaced with MAKE_GAUSSIAN_PLACE
// strcpy(rate_file, optarg);
break;
case 'P': AngularModulation = TRUE;
// CyclesPerMin = atof(optarg);
break;
case 'S': ; // Seed = atoi(optarg);
break;
case 'X': map_x = MAP_XSZ; //atof(optarg);
break;
case 'Y': map_y = MAP_YSZ; // atof(optarg);
break;
default:
instruct();
}
if (rate_map==NULL) {
if ((rate_map = (float **)calloc((size_t)YSZ, sizeof(float *))) == NULL) {
fprintf(stderr, "Can't allocate\n"); hxe(); }
for(i=0; i < YSZ; i++) if ((rate_map[i]=(float *)calloc(XSZ, sizeof(float))) == NULL){
fprintf(stderr, "Can't allocate\n"); hxe(); }
}
if (MAKE_GAUSSIAN_PLACE==1) {
// pick the center of the field
center_x = center_y = 0.0;
while(hypot(center_x - (double)(XSZ/2.0), center_y - (double)(YSZ/2.0)) > (XSZ/2.0)){
center_x = (double)(XSZ * drand48());
center_y = (double)(YSZ * drand48());
}
peak=PEAK; // AP/s
threshold = THRESHOLD; // minimum rate AP/s
sd = 0.2 * XSZ/2.0; // sd set to 20% of radius
MakeGaussianField((int)center_x, (int)center_y, sd, (int)XSZ, (int)YSZ, peak, threshold, rate_map);
collapse_x = collapse_y = 1.0;
storage[0]=center_x;storage[1]=center_y;storage[2]=collapse_x;storage[3]=collapse_y;
storage[4]=peak;storage[5]=threshold;storage[6]=sd;
} else if (MAKE_GAUSSIAN==2) { // will use to read from vector
for(y = 0; y < (int)map_y; y++)
for(x = 0; x < (int)map_x; x++)
fscanf(rfp,"%f", &(rate_map[y][x]));
collapse_x = (double)XSZ / (double)map_x;
collapse_y = (double)YSZ / (double)map_y;
} else { // MAKE_GAUSSIAN==0: resuse existing map
center_x=storage[0];center_y=storage[1];collapse_x=storage[2];collapse_y=storage[3];
peak=storage[4];threshold=storage[5];sd=storage[6];
}
nmb_header_lines = found = 0;
if (fflag) {
fseek(ifp,0,SEEK_SET);
fscanf(ifp, "%d\n", &nmb_header_lines);
for (i = 1; i < nmb_header_lines ; i++) {
get_line(ifp, line, 256);
// fprintf(ofp, "%s", line);
sscanf(line,"%s%lf", keyword, &value);
if(!strcmp(keyword,"%SAMPLING_INTERVAL(samps/sec)")){
samps_per_sec = value;
found++;
}else if(!strcmp(keyword, "%SCALE_Y(RatioTracktoMapPixels)")){
scale_y = value;
found++;
}else if(!strcmp(keyword, "%FIRST_TIMESTAMP(100usUnits)")){
first_time_stamp = value;
found++;
}else if(!strcmp(keyword, "%SCALE_X(RatioTracktoMapPixels)")){
scale_x = value;
found++;
}
}
if (found != 4) {
printf("PLACE simts() Only found %d/4. File is corrupt\n",found);
hxe(); }
}
// deliberately rounding down
timestamp_units_per_samp= (int)(TIMESTAMP_UNITS_PER_SEC/samps_per_sec);
if(AngularModulation) {
SampsPerCycle = (int)(samps_per_sec * 60.0 / CyclesPerMin);
SampsPerHalfCycle = (int)(SampsPerCycle / 2.0);
}
n_samps=0; modulate = 1.0;
while (1) {
if (fflag) {
if (fread(xys,sizeof(unsigned char), 3, ifp) != 3) break;
for(i = 0; i < xys[2]; i++) fread(&mt,sizeof(int), 1, ifp); // throw times away
x = xys[0]; y = xys[1];
} else {
if (n_samps>=nloc[0]) break;
x = loc[0][n_samps]; y=loc[1][n_samps];
}
n_samps++;
// determine rate modulation
if (AngularModulation){
samps = (n_samps % SampsPerCycle);
if(samps > SampsPerHalfCycle)
samps = SampsPerCycle - samps;
modulate = (double)samps / (double)SampsPerHalfCycle;
}
// use the collapse factors in case the rate map is a different resolution than the track
s = generate_spks(modulate*rate_map[(int)(y/collapse_y)][(int)(x/ collapse_x)],\
max_rate, samps_per_sec);
// xys[2] = s;
// write x y s -- CAN KILL (x loc, y loc, s # of spikes) ****************
// fwrite(xys, sizeof(unsigned char), (size_t)3, ofp);
if (sflag) {
if (s>0) {
xx[q]=(double)(n_samps + drand48())*timestamp_units_per_samp*0.1; // 0.1 converts to ms
yy[q]=(double)s;
q++;
}
} else {
sample_time = n_samps * timestamp_units_per_samp;
if (q>0 && s>0 && jit>0.) {
sample_time+=(jit*(1.-2.*drand48()));
while(sample_time<xx[q-1]) { // don't allow regress in time
sample_time=n_samps*timestamp_units_per_samp+(jit*(1.-2.*drand48())); }
}
for (i=0; i<s && q<maxsz; i++,q++) {
if (i==0) {
xx[q]=0.1*(sample_time+drand48()*timestamp_units_per_samp);
while (q>0 && xx[q]<xx[q-1]+REFRAC) xx[q]+=drand48()*timestamp_units_per_samp;
} else {
xx[q]=xx[q-1]+REFRAC+0.1*drand48()*timestamp_units_per_samp;
}
if (DEBUG==1) printf("%d:%20.2lf:%20.2lf %d\n",q,sample_time,xx[q],s);
}
}
if (q>maxsz) break;
// write spike times -- CHANGE HERE ****************
// fwrite(time, sizeof(int), (size_t)s, ofp); // this is the time in time
}
vector_resize(vv, q);
if (sflag) vector_resize(vv1, q);
return (double)q;
}
ENDVERBATIM
:* outvec.kendall(srcvecs,tmpvecs)
VERBATIM
static double kendall (void* vv) {
ui4 i, j, k, m, n, p, q, maxsz, nx, av[VRRY], bv[VRRY], num, numb, indflag, match;
ui4 MaxBins, NumPairs, *SpikeCounts;
si1 **Diffs, dif;
Object *ob, *ob2;
double *x, *avo[VRRY], *bvo[VRRY], val[VRRY];
double IntStart, IntervalStop, AnInt, MaxT;
double tau, z, prob;
void *vva[VRRY],*vvb[VRRY];
void KendallTauPairs();
nx = vector_instance_px(vv, &x);
ob = *hoc_objgetarg(1);
ob2 = *hoc_objgetarg(2);
num = ivoc_list_count(ob);
numb = ivoc_list_count(ob2);
if (num>VRRY) hoc_execerror("kendall ****ERRA****: can only handle VRRY vectors", 0);
if (numb!=4) hoc_execerror("kendall ****ERRB****: 2nd list count should be 4:i,j,tau,p", 0);
for (i=0;i<num;i++) av[i]=list_vector_px2(ob, i, &avo[i], &vva[i]);
for (i=0;i<numb;i++) {
bv[i]=list_vector_px3(ob2, i, &bvo[i], &vvb[i],1);
if (bv[0]!=bv[i]) { printf("kendall ****ERRC2**** %d %d %d\n",i,bv[0],bv[i]);
hoc_execerror("Vectors must all be same size: ", 0); }
}
Diffs = (si1 **)calloc(num, sizeof(si1 *));
AnInt=AnalysisInterval*1e3; // convert from sec to ms
MaxT=MaxTime*1e3; // convert from sec to ms
MaxBins = (ui4) AnInt/BinSize;
NumPairs = (ui4)MaxBins*(MaxBins-1); // not just neighboring bins
SpikeCounts = (ui4 *)calloc(MaxBins, sizeof(ui4));
for(i=0;i<num;i++) {
Diffs[i] = (si1 *)calloc(NumPairs, sizeof(si1));
}
// do correlation over segments of the data between Start and Stop
// MaxTime is the full duration of the recording
IntStart = IntervalStart*1e3; // s to ms
IntervalStop = IntStart + AnInt;
while (IntervalStop <= MaxT) {
// Initialize the Spike Count bins
for (i=0;i<num;i++) {
for(j=0;j<MaxBins;j++) SpikeCounts[j]=0;
for (j=0;j<av[i];j++) {
if (avo[i][j]<IntStart) continue;
if (avo[i][j]>IntervalStop) break;
SpikeCounts[(int)((avo[i][j]-IntStart)/BinSize)]++;
}
// PairwiseDiffSigns(SpikeCounts[i], Diffs[i], MaxBins, NumPairs);
// SpikeCounts->data, Diffs->signs, MaxBins->NumEvents, NumPairs->NumPairs
for (k=0,n=0;k<MaxBins;k++) {
for (j=(k+1);j<MaxBins;j++,n++) {
if (dif=(SpikeCounts[k]-SpikeCounts[j])) {
Diffs[i][n]=(dif>0)?1:-1;
} else Diffs[i][n]=0;
}
}
}
// calculate Kendall's tau for each pair of time series
for(i=0,k=0; i<num; i++) {
for(j=i+1; j<num; j++,k++) {
KendallTauPairs(Diffs[i], Diffs[j], MaxBins, NumPairs, &tau, &z, &prob);
if (VERBOSE) printf("%g %g %d %d %0.6lf %0.3lf %0.3lf\n",\
IntStart/1e3, IntervalStop/1e3, i, j, tau, z, prob);
// i j tau prob -- save i, j, tau, prob in an NQS given by 2nd list arg
if (k>=bv[0]){printf("PLACE kendall() ERRD out of room: %d %d\n",k,bv[0]); hxe();}
bvo[0][k]=(double)i;bvo[1][k]=(double)j;bvo[2][k]=tau;bvo[3][k]=prob;
}
}
IntStart+=AnInt;
IntervalStop +=AnInt;
}
for (j=0;j<numb;j++) vector_resize(vvb[j], k);
return (double)k;
}
static double kendal2 (void* vv) {
ui4 i, n, j, k, MaxBins, NumPairs;
double *x, *y, tau, z, prob;;
si1 *Dx, *Dy;
si4 difx, dify;
ui4 LeftNum = 0, RightNum = 0;
si4 Numerator = 0;
si1 dxdy;
sf8 Variance;
ui4 NumEvents;
MaxBins = vector_instance_px(vv, &x);
if ((i=vector_arg_px(1,&y)) != MaxBins ) {printf("kenall ERRA: %d %d\n",n,MaxBins); hxe();}
NumPairs = (ui4)MaxBins*(MaxBins-1)/2; // not just neighboring bins
if(VERBOSE>1) printf("NumPairs=%d, MaxBins=%d\n",NumPairs,MaxBins);
Dx = (si1 *)calloc(NumPairs, sizeof(si1)); // sets memory to 0
if(!Dx) { printf("kendal2 ERR0: out of mem!\n"); hxe(); }
Dy = (si1 *)calloc(NumPairs, sizeof(si1));
if(!Dy) { printf("kendal2 ERR1: out of mem!\n"); hxe(); }
for (k=0,n=0;k<MaxBins;k++) for (j=(k+1);j<MaxBins;j++,n++) {
if (difx=x[k]-x[j]) Dx[n]=(difx>0)?1:-1;
if (dify=y[k]-y[j]) Dy[n]=(dify>0)?1:-1;
}
if(VERBOSE) printf("n=%d, NumPairs=%d\n",n,NumPairs);
KendallTauPairs(Dx, Dy, MaxBins, NumPairs, &tau, &z, &prob);
if(VERBOSE) printf("tau:%g z:%g p:%g\n",tau, z, prob);
free(Dx); free(Dy);
return tau;
}
// kend2() is like kendall but only looks at neighboring pairs
static double kend2 (void* vv) {
ui4 i, n, j, k, pw, MaxBins, NumPairs;
double *x, *y, tau, z, prob;;
si1 *Dx, *Dy;
si4 difx, dify;
MaxBins = vector_instance_px(vv, &x);
if ((i=vector_arg_px(1,&y)) != MaxBins ) {printf("kend2 ERRA: %d %d\n",n,MaxBins); hxe();}
if (ifarg(2)) pw=1; else pw=0; // pairwise flag
NumPairs = (ui4)(MaxBins-1); // just neighboring bins
Dx = (si1 *)calloc(NumPairs, sizeof(si1)); // sets memory to 0
Dy = (si1 *)calloc(NumPairs, sizeof(si1));
for (k=0,n=0;k<MaxBins-1;k++,n++) { j=(k+1);
if (difx=x[k]-x[j]) Dx[n]=(difx>0)?1:-1;
if (dify=y[k]-y[j]) Dy[n]=(dify>0)?1:-1;
}
KendallTauPairs(Dx, Dy, MaxBins, NumPairs, &tau, &z, &prob);
if(VERBOSE>0) printf("tau:%g z:%g p:%g\n",tau, z, prob);
free(Dx); free(Dy);
return tau;
}
// Given data arrays x[1..n] and y[1..n], this program returns Kendall's correlation tau,
// its number of standard deviations from zero as z, and its two-tailed probability level as prob.
// Small values of prob indicate a significant correlation (tau positive) or anticorrelation (tau
// negative).
//
// the necessary formulas are:
// tau = numerator / denominator
// numerator = (concordant_pairs - discordant_pairs)
// denominator = sqrt(concordant_pairs + discordant_pairs+ extra-y) * sqrt(concordant_pairs + discordant_pairs+ extra-x)
// denominator = sqrt(LeftNum) * sqrt(RightNum)
// where concordant pairs are those differences (xi - xj) that have the same sign as the differences (yi - yj)
// extra-y are those non-zero differences (yi - yj) where (xi - xj) = 0
// extra-x are those non-zero differences (xi - xj) where (yi - yj) = 0
// variance(tau) = (4N + 10) / 9N(N - 1)
// z = tau / sqrt(variance)
// Given arrays dx[1..n] and dy[1..n] which represent the signs (-1, 0, 1) of all
// pairs in the data sets x and y,
// this program returns Kendall's correlation tau. dx and dy are the result of
// calling PairwiseDiffSigns
// its number of standard deviations from zero as z, and its two-tailed probability level as prob.
// Small values of prob indicate a significant correlation (tau positive) or
// anticorrelation (tau negative).
void KendallTauPairs (si1 *dx, si1 *dy, ui4 NumEvents, ui4 NumPairs, sf8 *tau, sf8 *z, sf8 *prob)
{
sf8 erfcc(sf8 x);
ui4 i;
ui4 LeftNum = 0, RightNum = 0;
si4 Numerator=0;
si1 dxdy;
sf8 Variance;
for (i=0; i < NumPairs; i++) {
dxdy = dx[i] * dy[i];
if (dxdy) { // Neither array has a tie.
LeftNum++;
RightNum++;
dxdy > 0.0 ? Numerator++ : Numerator--;
} else { // One or both arrays have ties.
if(dy[i])
LeftNum++; //an extra y event
else if(dx[i])
RightNum++; //an extra x event
}
}
if(LeftNum && RightNum){
*tau = (sf8) Numerator / (sqrt((sf8)LeftNum) * sqrt((sf8)RightNum));
Variance = (sf8)(4.0 * NumEvents + 10.0) / (sf8)(9.0 * NumEvents * (NumEvents - 1.0));
*z = (*tau) / sqrt(Variance);
*prob = erfcc(fabs(*z)/1.4142136);
}else{
*tau = *z = 0.0;
*prob = 1.0;
}
return;
}
// Returns the complementary error function erfc(x) with fractional error
// everywhere less than 1.2 x 10^-7.
double erfcc(double x)
{
double mt,z,ans;
z=fabs(x);
mt=1.0/(1.0+0.5*z);
ans=mt*exp(-z*z-1.26551223+mt*(1.00002368+mt*(0.37409196+mt*(0.09678418+ mt*(-0.18628806+mt*(0.27886807+mt*(-1.13520398+mt*(1.48851587+ mt*(-0.82215223+mt*0.17087277)))))))));
return x >= 0.0 ? ans : 2.0-ans;
}
ENDVERBATIM
VERBATIM
static double mkgaussfield(void* vv) {
// NOT DEBUGGED
double x,y,peak,threshold,sd,d, *map;
int i,j,nx,maxsz;
nx = vector_instance_px(vv, &map);
maxsz=vector_buffer_size(vv);
if (maxsz<XSZ*YSZ){printf("mkgasussfield ERR vector too small %d<%d",maxsz,XSZ*YSZ); hxe();}
vector_resize(vv, XSZ*YSZ);
// pick the center of the field
x = y = 0.0;
while(hypot(x - (double)(XSZ/2.0), y - (double)(YSZ/2.0)) > (XSZ/2.0)){
x = (double)(XSZ * drand48());
y = (double)(YSZ * drand48());
}
peak = 60.0; // AP/s
threshold = 0.5; // minimum rate AP/s
sd = 0.2 * XSZ/2.0; // sd set to 20% of radius
for (i=0; i<YSZ; i++) { // YSZ rows
for (j=0; j<XSZ; j++) { // XSZ cols
d = hypot((double)(i-y), (double)(j-x));
map[i*XSZ+j] = peak * gaussaf(0.0, sd, d);
// if (map[i*XSZ+j]!=0) printf("BB: %d %d %f\n", i,j,map[i*XSZ+j]);
if(map[i*XSZ+j] < threshold)
map[i*XSZ+j] = 0.0;
}
}
}
ENDVERBATIM
VERBATIM
static double dumpratemap(void* vv) {
int i,j,k,nx,maxsz;
double *map;
if (rate_map==NULL) {printf("rate_map not allocated\n"); hxe();}
nx = vector_instance_px(vv, &map);
maxsz=vector_buffer_size(vv);
if (maxsz<XSZ*YSZ){printf("dumpratemap ERR vector too small %d<%d\n",maxsz,XSZ*YSZ); hxe();}
vector_resize(vv, XSZ*YSZ);
for (i=0; i<YSZ; i++) { // YSZ rows
for (j=0; j<XSZ; j++) { // XSZ cols
map[i*XSZ+j]=rate_map[i][j];
}
}
}
ENDVERBATIM
: from /u/billl/nrniv/place/fenton/read_time_series/read_time_series.c
: took out the BIG_ENDIAN byte reversal stuff
:* v.rd()
VERBATIM
static double rd (void* vv) {
unsigned char x, y, s;
unsigned int i, j, n, time, beg, nx,ny, maxnx, locflag, linflag;
unsigned int maxsz, k, nmb_header_lines, found;
int err;
double *xx, *xloc, *yloc, *xlin, value, errtol;
void *vv1, *vv2, *vv4;
char line[4196], keyword[4196];
FILE* fp, *hoc_obj_file_arg();
fp = hoc_obj_file_arg(1);
vector_instance_px(vv, &xx);
if((maxsz=vector_buffer_size(vv))<1) {maxsz=100; xx=vector_newsize(vv,maxsz);} else vector_resize(vv, maxsz);
err=-1; maxnx=locflag=linflag=0;
if (ifarg(3)) {
nx = vector_arg_px(2, &xloc); vv1=vector_arg(2);
maxnx=vector_buffer_size(vv1);
ny = vector_arg_px(3, &yloc); vv2=vector_arg(3);
k=vector_buffer_size(vv2);
if (k!=maxnx) {printf("rd_PLACE: ERRA not same size: %d %d\n",k,maxnx); hxe();}
vector_resize(vv1, maxnx); vector_resize(vv2, maxnx);
locflag=1;
} else if (ifarg(2) && hoc_is_double_arg(2)) {
errtol=*getarg(2); err=0;
} else if (ifarg(2)) { linflag=1;
vv1=vector_arg(2);
xlin=vector_newsize(vv1,maxsz);
}
fseek(fp,0,SEEK_SET);
found=0;
fscanf(fp, "%d\n", &nmb_header_lines);
for (i = 1; i < nmb_header_lines ; i++) {
get_line(fp, line, 256);
sscanf(line,"%s%lf", keyword, &value);
if(!strcmp(keyword,"%SAMPLING_INTERVAL(samps/sec)")){
samps_per_sec = value;
found++;
}else if(!strcmp(keyword, "%SCALE_Y(RatioTracktoMapPixels)")){
scale_y = value;
found++;
}else if(!strcmp(keyword, "%FIRST_TIMESTAMP(100usUnits)")){
first_time_stamp = value;
found++;
}else if(!strcmp(keyword, "%SCALE_X(RatioTracktoMapPixels)")){
scale_x = value;
found++;
}
}
if (found != 4) {
printf("PLACE rd() Only found %d/4. File is corrupt\n",found);
hxe(); }
for (n=0,j=0,k=0;fread(&x,sizeof(unsigned char), 1, fp) != 0; n++) {
fread(&y,sizeof(unsigned char), 1, fp);
fread(&s,sizeof(unsigned char), 1, fp);
for (i=0; i<s; i++) {
if(!fread(&time,sizeof(unsigned int), 1, fp))
printf("rd_PLACE: WARNZ couldnt read spiketime %d of %d at byte %d\n",i+1,s,ftell(fp));
if (k>=maxsz) {
xx=vector_newsize(vv,maxsz*=2);
if (linflag) xlin=vector_newsize(vv1,maxsz);
}
if (linflag) xlin[k]=(double)n+1.;
xx[k++]=0.1*(double)(time); // ms
if (err>=0 &&
(xx[k-1]<((n-1)*1e3/samps_per_sec-errtol)||xx[k-1]>n*1e3/samps_per_sec+errtol)) {
if (VERBOSE) {
printf("rdERRB:%d:%g (%lf:%lf) ",n,xx[k-1],((n-1)*1e3/samps_per_sec),n*1e3/samps_per_sec);
}
err++;
}
}
if (locflag) {
if (j>=maxnx) {
maxnx*=2;
xloc=vector_newsize(vv1,maxnx);
yloc=vector_newsize(vv2,maxnx);
}
xloc[j]=(double)x; yloc[j]=(double)y; j++;
}
}
vector_resize(vv, k);
if (linflag) vector_resize(vv1,k);
if (locflag) {vector_resize(vv1, j); vector_resize(vv2, j);}
if (err>0) return (double)-err; else return (double)n*1e3/samps_per_sec;
}
ENDVERBATIM
: modified from rd()
:* spks.rdtts(file,[,xloc,yloc,tloc])
: read spike times into spks vector; optionally read locations and times into 3 other vecs
VERBATIM
static double rdtts (void* vv) {
unsigned char x, y, s;
unsigned int i, j, n, time, tt, beg, nx,ny, maxnx, locflag;
unsigned int maxsz, k, nmb_header_lines, found;
int err;
double *xx, *tvec, *xloc, *yloc, value, errtol;
void *vv1, *vv2, *vv3;
char line[4196], keyword[4196];
FILE* fp, *hoc_obj_file_arg();
fp = hoc_obj_file_arg(1);
vector_instance_px(vv, &xx);
if((maxsz=vector_buffer_size(vv))<1) {maxsz=100; xx=vector_newsize(vv,maxsz);} else vector_resize(vv, maxsz);
err=-1;
if (ifarg(3)) {
nx = vector_arg_px(2, &xloc); vv1=vector_arg(2);
maxnx=vector_buffer_size(vv1);
ny = vector_arg_px(3, &yloc); vv2=vector_arg(3);
ny = vector_arg_px(4, &tvec); vv3=vector_arg(4);
xloc=vector_newsize(vv1,maxnx); yloc=vector_newsize(vv2,maxnx);
tvec=vector_newsize(vv3,maxnx); locflag=1;
} else {
maxnx=locflag=0;
if (ifarg(2)) errtol=*getarg(2); err=0;
}
fseek(fp,0,SEEK_SET);
found=0;
fscanf(fp, "%d\n", &nmb_header_lines);
for (i = 1; i < nmb_header_lines ; i++) {
get_line(fp, line, 256);
sscanf(line,"%s%lf", keyword, &value);
if(!strcmp(keyword,"%SAMPLING_INTERVAL(samps/sec)")){
samps_per_sec = value; found++;
}else if(!strcmp(keyword, "%SCALE_Y(RatioTracktoMapPixels)")){
scale_y = value; found++;
}else if(!strcmp(keyword, "%FIRST_TIMESTAMP(100usUnits)")){
first_time_stamp = value; found++;
}else if(!strcmp(keyword, "%SCALE_X(RatioTracktoMapPixels)")){
scale_x = value; found++;
}
}
if (found!=4){printf("PLACE rd() Only found %d/4. File is corrupt\n",found); hxe();}
for (n=0,j=0,k=0;fread(&tt,sizeof(unsigned int), 1, fp) != 0; n++) {
fread(&x,sizeof(unsigned char), 1, fp);
fread(&y,sizeof(unsigned char), 1, fp);
fread(&s,sizeof(unsigned char), 1, fp);
for (i=0; i<s; i++) {
if(!fread(&time,sizeof(unsigned int), 1, fp))
printf("rd_PLACE: WARNZ couldnt read spiketime %d of %d at byte %d\n",i+1,s,ftell(fp));
if (k>=maxsz) xx=vector_newsize(vv,maxsz*=2);
xx[k++]=0.1*(double)(time); // ms
if (err>=0 &&
(xx[k-1]<((n-1)*1e3/samps_per_sec-errtol)||xx[k-1]>n*1e3/samps_per_sec+errtol)) {
if (VERBOSE) {
printf("rdERRB:%d:%g (%lf:%lf) ",n,xx[k-1],((n-1)*1e3/samps_per_sec),n*1e3/samps_per_sec);
}
err++;
}
}
if (locflag) {
if (j>=maxnx) {maxnx*=2; tvec=vector_newsize(vv3,maxnx);
xloc=vector_newsize(vv1,maxnx); yloc=vector_newsize(vv2,maxnx);
}
tvec[j]=0.1*(double)tt; xloc[j]=(double)x; yloc[j]=(double)y; j++;
}
}
vector_resize(vv, k);
if (locflag) {vector_resize(vv1, j); vector_resize(vv2, j); vector_resize(vv3, j);}
if (err>0) return (double)-err; else return (double)n*1e3/samps_per_sec;
}
ENDVERBATIM
VERBATIM
//more bpf related stuff
si1 *sscanfLine(si1 *dp, si1 *s)
{
int i = 0;
*s = (si1) 0x0;
while (sscanf(dp, "%c", s + i) != EOF){
dp++;
if(*(s + i) == '\n'){
*(s + i + 1) = '\0';
return s;
}
++i;
}
*(s + i) = '\0';
return s;
}
#define DEBUGPR(x) printf("HERE %d\n", x); fflush(stdout);
si4 ReadBPFHeader (ui1 *dp, ui4 *BPFRecordSizes, si2 *Gains, si4 *AcquisitionSystem, ui4 *EEGChList, ui4 *BPFRecordTypeNumbers )
{
ui1 VerifyRecordFormat(ui1 *, ui4 *, ui1 *, si1 *, si1 *, ui1);
ui1 *FindHeaderStart(), *FindHeaderEnd(), *FindSection(), *FindSectionEnd();
void GetRecordSizes();
ui1 BPFRecordTypes[256];
ui1 *HeaderStart, *HeaderEnd, *SectionStart, *SectionEnd;
si1 KeyWordString[256], Encoding[256], **Parameters;
si4 Type, NumberOfParameters, NumEEGChannels;
si4 GetRecordTypesUsed(), FindKeyWord();
// Look for header beginning
HeaderStart = FindHeaderStart(dp);
if(HeaderStart == NULL){
fprintf(MyStdErr,"\nCan't find start of header (%s).\nWrong File Format.\n", HEADER_BEGIN_LABEL);
return(-1);
}
dp = (HeaderStart + strlen(HEADER_BEGIN_LABEL));
// Check to see if it's the BPF header or not
sscanf((si1 *)dp,"%s", KeyWordString);
if(strstr(KeyWordString, BPF_HEADER_TYPE) == NULL){
fprintf(MyStdErr,"This is not a BPF file. It is a %s file.\n", KeyWordString);
return(-1);
}
// Section: DATABASE_INFORMATION
dp = HeaderStart;
SectionStart = FindSection(dp, SECTION_DATABASE_INFORMATION);
if(SectionStart == NULL){
fprintf(MyStdErr, "\nCan't find Header section %s.\nWrong File Format.\n", SECTION_DATABASE_INFORMATION);
return(-1);
}
dp = SectionStart;
// Acquisition System
*AcquisitionSystem = ACX;
if(!FindKeyWord(SectionStart, ACQUISITION_SYSTEM_LABEL, KeyWordString, Encoding)){
printf("AcX Acquisition System\n");
}else{
Type = GetType(KeyWordString);
switch(Type){
case -1:
fprintf(MyStdErr,"KeyWord.integer format not respected in: %s\nWrong File Format.\n", KeyWordString);
return(0);
case 0: // ASCII encoded
NumberOfParameters = GetASCIIEncoding(Encoding, &Parameters);
if(NumberOfParameters == 1){
if(strcasestr(Parameters[0], "Axona")){
*AcquisitionSystem = AXONA;
if(VERBOSE) printf("AXONA Acquisition System\n");
}
}
break;
default:
fprintf(MyStdErr,"%s must be ASCII encoded\nWrong File Format.\n", ACQUISITION_SYSTEM_LABEL);
return(0);
}
}
// Section: SETUP_INFORMATION
dp = HeaderStart;
SectionStart = FindSection(dp, SECTION_SETUP_INFORMATION);
if(SectionStart == NULL){
fprintf(MyStdErr, "\nCan't find Header section %s.\nWrong File Format.\n", SECTION_SETUP_INFORMATION);
return(-1);
}
dp = SectionStart;
NumEEGChannels = GetRecordTypesUsed(BPFRecordTypeNumbers, dp, EEGChList);
if(!GetGainsList(Gains, dp)) return -1;
SectionEnd = FindSectionEnd(SectionStart);
if(SectionEnd == NULL)
fprintf(MyStdErr,"\nCan't find end of Header Section (%s).\nWrong File Format.\n", SECTION_SETUP_INFORMATION);
// Section: RECORD_FORMAT
dp = HeaderStart;
SectionStart = FindSection(dp, SECTION_RECORD_FORMAT_INFORMATION);
if(SectionStart == NULL){
fprintf(MyStdErr, "\nCan't find section %s. Wrong data format\n", SECTION_RECORD_FORMAT_INFORMATION);
return (-1);
}
// EEG RECORD ID
if(BPFRecordTypeNumbers[EEG_BPF_REC_TYPE]){
BPFRecordSizes[EEG_BPF_REC_TYPE] = BPFRecordTypeNumbers[EEG_BPF_REC_TYPE] * EEG_BPF_REC_DATA_SIZE + EEG_BPF_REC_ID_SIZE;
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, EEG_ID, EEG_RECORD_FORMAT, EEG_BPF_REC_TYPE))
return(-1);
}
// SINGLE ELECTRODE RECORD ID
if(BPFRecordTypeNumbers[SINGLE_BPF_REC_TYPE]){
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, SINGLE_ID, SINGLE_RECORD_FORMAT, SINGLE_BPF_REC_TYPE))
return(-1);
}
// STEREOTRODE RECORD ID
if(BPFRecordTypeNumbers[STEREO_BPF_REC_TYPE]){
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, STEREO_ID, STEREO_RECORD_FORMAT, STEREO_BPF_REC_TYPE))
return(-1);
}
// TETROTRODE RECORD ID
if(BPFRecordTypeNumbers[TETRODE_BPF_REC_TYPE]){
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, TETRODE_ID, TETRODE_RECORD_FORMAT, TETRODE_BPF_REC_TYPE))
return(-1);
}
// SYNC RECORD ID
if(BPFRecordTypeNumbers[SYNC_BPF_REC_TYPE]){
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, SYNC_ID, SYNC_RECORD_FORMAT, SYNC_BPF_REC_TYPE))
return(-1);
}
// ROOM POSITION RECORD ID
if(BPFRecordTypeNumbers[ROOM_POSITION_BPF_REC_TYPE]){
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, ROOM_POSITION_ID, ROOM_POSITION_RECORD_FORMAT, ROOM_POSITION_BPF_REC_TYPE))
return(-1);
}
// ARENA POSITION RECORD ID
if(BPFRecordTypeNumbers[ARENA_POSITION_BPF_REC_TYPE]){
if(!VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, ARENA_POSITION_ID, ARENA_POSITION_RECORD_FORMAT, ARENA_POSITION_BPF_REC_TYPE))
return(-1);
}
// KEY EVENT RECORD ID
if(BPFRecordTypeNumbers[KEY_EVENT_BPF_REC_TYPE]){
VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, KEY_EVENT_ID, KEY_EVENT_RECORD_FORMAT, KEY_EVENT_BPF_REC_TYPE);
// This record may not be defined
}
// INPUT EVENT RECORD ID
if(BPFRecordTypeNumbers[INPUT_EVENT_BPF_REC_TYPE]){
VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, INPUT_EVENT_ID, INPUT_EVENT_RECORD_FORMAT, INPUT_EVENT_BPF_REC_TYPE);
// This record may not be defined
}
// OUTPUT EVENT RECORD ID
if(BPFRecordTypeNumbers[OUTPUT_EVENT_BPF_REC_TYPE]){
VerifyRecordFormat(BPFRecordTypes, BPFRecordSizes, SectionStart, OUTPUT_EVENT_ID, OUTPUT_EVENT_RECORD_FORMAT, OUTPUT_EVENT_BPF_REC_TYPE);
// This record may not be defined
}
SectionEnd = FindSectionEnd(SectionStart);
if(SectionEnd == NULL)
fprintf(MyStdErr,"\nCan't find end of Header Section (%s).\nWrong File Format.\n", SECTION_RECORD_FORMAT_INFORMATION);
HeaderEnd = FindHeaderEnd(SectionEnd);
if(HeaderEnd == NULL){
fprintf(MyStdErr, "\nCan't find end of header. Wrong data format\n");
return (-1);
}
return(NumEEGChannels);
}
ui1 *FindHeaderStart(ui1 *mdata)
{
si1 String[4196], Line[4196], *HeaderStart;
ui1 *dp;
dp = mdata;
while(sscanf((si1 *)dp, "%s", String) == 1){
HeaderStart = (si1 *)strstr(String, HEADER_BEGIN_LABEL);
if(HeaderStart != NULL){
return(dp + (HeaderStart - String));
}
sscanfLine(dp, Line);
dp += strlen(Line);
}
fprintf(MyStdErr,"Can't find BPF Header begin label: %s\n", HEADER_BEGIN_LABEL);
return(NULL);
}
ui1 *FindHeaderEnd(ui1 *mdata)
{
si1 String[4196], Line[4196];
ui1 *dp;
dp = mdata;
while(sscanfLine((si1 *)dp, Line) != NULL){
dp = (ui1 *)strstr((si1 *)dp, Line);
dp += strlen(Line);
sscanf(Line, "%s", String);
if(strstr(String, HEADER_END_LABEL) != NULL){
return(dp);
}
}
fprintf(MyStdErr,"Can't find BPF Header end label: %s\n", HEADER_END_LABEL);
return(NULL);
}
ui1 *FindSection(ui1 *dp, ui1 *Section)
{
si1 Flag[4196];
while (sscanf((si1 *)dp,"%s", Flag) != EOF){
dp = (ui1 *)strstr((si1 *)dp, Flag);
dp += strlen(Flag);
if(Comment((si1 *)dp, Flag)) continue;
if(!strcmp(Flag, SECTION_BEGIN_LABEL)){
sscanf((si1 *)dp,"%s", Flag);
dp = (ui1 *)strstr((si1 *)dp, Flag);
dp += strlen(Flag);
if(!strcmp(Flag, (si1 *)Section))
return(dp);
}
}
return NULL;
}
si4 FindKeyWord(si1 *dp, si1 *KeyWord, si1 *KeyWordString, si1 *Encoding)
{
si4 SectionEnd(), n_chars;
while (sscanf(dp,"%s", KeyWordString) != EOF){
dp = (si1 *)strstr(dp, KeyWordString);
dp += strlen(KeyWordString);
if(Comment(dp, KeyWordString)) continue;
if(SectionEnd(KeyWordString)) break;
if(strstr(KeyWordString, KeyWord) != NULL){
// search until the "." for a paren;
// Do this because some BPF files don't have a space between the Keyword and the encoding info
// for example ListOfGains.0(2 2 2 2 etc
n_chars = 0;
while(*(dp-n_chars) != '.'){
if(*(dp-n_chars) == '('){
dp -= n_chars;
break;
}
n_chars++;
}
(void)sscanfLine(dp, Encoding);
return(1);
}
}
return 0;
}
ui1 *FindSectionEnd(ui1 *dp)
{
si1 Line[4196], String[4196];
while (sscanfLine(dp, Line) != EOF){
dp = (ui1 *)strstr((si1 *)dp, Line);
dp+= strlen(Line);
sscanf(Line,"%s",String);
if(SectionEnd(String)){
return(dp);
}
}
return(NULL);
}
si4 Comment(ui1 *dp, si1 *String)
{
si1 Line[4196];
if(!strcmp(String, COMMENT)){
sscanfLine(dp, Line);
dp += strlen(Line);
return(1);
}else
return(0);
}
si4 SectionEnd(si1 *String)
{
if(!strcmp(String, SECTION_END_LABEL))
return(1);
return(0);
}
si4 GetType(si1 *KeyWord)
{
si1 *TypeString;
si4 Type = -1;
TypeString = (si1 *)strchr(KeyWord, TYPE_PREFIX);
if(TypeString == NULL)
return(-1);
sscanf(TypeString+1,"%d", &Type); // +1 to step past the '.'
return(Type);
}
si4 GetBinaryEncoding(si1 *Encoding, si1 *P)
{
static si1 String[4196], *TerminateString, *Line;
// The character pair (' must start the Binary Encoding
if((Line = (si1 *)strstr(Encoding, BINARY_STARTING_MARK)) == NULL){
fprintf(MyStdErr,"a (' must start binary encodings.\n");
return (0);
}
Line+= strlen(BINARY_STARTING_MARK);;
sscanf(Line,"%c",P);
if(P == NULL){
fprintf(MyStdErr,"Can't read Binary code.\n");
return (0);
}
return(1);
}
si4 GetASCIIEncoding(si1 *Encoding, si1 ***P)
{
si4 n;
static si1 String[4196], *TerminateString, *Line, **Parameters;
Parameters = (si1 **)calloc(MAX_NUMBER_BPF_FORMAT_PARAMETERS, sizeof(si1 *));
if(Parameters == NULL){
fprintf(MyStdErr, "\nCan't allocate Parameters\n");
return (0);
}
*P = Parameters;
// a '(' must start the ASCII Encoding
if((Line = (si1 *)strstr(Encoding, ASCII_STARTING_MARK )) == NULL){
fprintf(MyStdErr,"a '(' must start ASCII encodings. Found: %s\n", Encoding);
return (0);
}
Line+= strlen(ASCII_STARTING_MARK);;
if(sscanf(Line, "%s", String) != 1){
fprintf(MyStdErr,"No ASCII string to decode\n");
return (0);
}
Parameters[0] = (si1 *)calloc(strlen(String + 1), sizeof(si1));
strcpy(Parameters[0], String);
if((TerminateString = (si1 *)strstr(Parameters[0], ASCII_TERMINATING_MARK )) != NULL){ // a ')' ends the encoding
*TerminateString = '\0';
return(1);
}
n = 1;
Line += (strlen(String) + 1);
while((*Line == ' ') || (*Line == '\t')){ // move to start of next string
Line ++;
}
while(sscanf(Line, "%s", String) != EOF){
if(!strcmp(String, ")"))
return (n);
Parameters[n] = (si1 *)calloc(strlen(String) + 1, sizeof(si1));
strcpy(Parameters[n], String);
if((TerminateString = (si1 *)strchr(Parameters[n], ')')) != NULL){ // a ')' indicates the end of the ASCII encoding
*TerminateString = '\0' ;
return(++n);
}
Line += (strlen(String) + 1);
while((*Line == ' ') || (*Line == '\t')){ // move to start of next string
Line ++;
}
n++;
if(n == MAX_NUMBER_BPF_FORMAT_PARAMETERS){
fprintf(MyStdErr,"Too many Parameters in Encoding of:\t");
return(0);
}
}
return (n);
}
ui1 VerifyRecordFormat(ui1 *RecordTypes, ui4 *RecordSizes, ui1 *SectionStart, si1 *id, si1 *RecordFormat, ui1 RecType){
ui1 *dp;
si1 KeyWordString[4196], Encoding[4196];
si4 Type;
// Get record type
dp = SectionStart;
if(!FindKeyWord((si1 *)dp, id, KeyWordString, Encoding)){
fprintf(MyStdErr,"\nCan't find keyWord (%s).\nWrong File Format.\n", id);
return 0;
}
Type = GetType(KeyWordString);
switch(Type){
case -1:
fprintf(MyStdErr,"KeyWord.integer format not respected in: %s\nWrong File Format.\n", KeyWordString);
return 0;
case 0: // ASCII encoded
fprintf(MyStdErr,"%s\nWrong File Format.\n", id);
return 0;
case RECORD_ID_ENCODING: // 1 byte Binary encoded
if(!GetBinaryEncoding(Encoding, (si1 *)RecordTypes+RecType)){
fprintf(MyStdErr,"\nCan't get RECORD ID from header.\n");
return 0;
}
if(RecordTypes[RecType] != RecType){
fprintf(MyStdErr,"\nRecord type (%c) does not match line in the header.\n", RecType);
return(0);
}
break;
default:
fprintf(MyStdErr,"Record id (%s) must be Binary encoded\nWrong File Format.\n", Type, KeyWordString);
return (0);
}
//Get record size
dp = SectionStart;
if(!FindKeyWord((si1 *)dp, RecordFormat, KeyWordString, Encoding)){
fprintf(MyStdErr,"\nCan't find keyWord (%s).\nWrong File Format.\n", RecordFormat);
return 0;
}
RecordSizes[RecType] = (ui4)GetType(KeyWordString);
/* if (RecordSizes[RecType] != (ui4)GetType(KeyWordString)){
fprintf(MyStdErr,"\nRecord size (%d) does not match line in the header (%s).\n", RecordSizes[RecType], KeyWordString);
return(0);
}
*/
return(1);
}
si4 GetRecordTypesUsed(ui4 *RecordTypeNumbers, si1 *SectionStart, ui4 *BPFEEGChannels){
ui4 GetNumberOfChannels(), GetIfPositions();
si1 *dp;
ui4 i, n;
si4 NumEEGChannels;
si4 GetEEGChannelList();
for(i = 0; i < MAX_BPF_RECORD_TYPES; i++)
RecordTypeNumbers[i] = 0;
n = 0;
dp = SectionStart;
n += GetNumberOfChannels(RecordTypeNumbers, dp, NUMBER_OF_EEG_CHANNELS, EEG_BPF_REC_TYPE);
if(n){ // There are EEG channels
NumEEGChannels = GetEEGChannelList(BPFEEGChannels, dp);
}else
BPFEEGChannels = NULL;
n += GetNumberOfChannels(RecordTypeNumbers, dp, NUMBER_OF_SINGLE_CHANNELS, SINGLE_BPF_REC_TYPE);
n += GetNumberOfChannels(RecordTypeNumbers, dp, NUMBER_OF_STEREO_CHANNELS, STEREO_BPF_REC_TYPE);
n += GetNumberOfChannels(RecordTypeNumbers, dp, NUMBER_OF_TETRODE_CHANNELS, TETRODE_BPF_REC_TYPE);
if(!n){
fprintf(MyStdErr,"\nNo indication in header that there are electrophysiological data.\nWrong File Format.\n");
return (NumEEGChannels);
}
(void) GetNumberOfChannels(RecordTypeNumbers, dp, NUMBER_OF_SYNC_CHANNELS, SYNC_BPF_REC_TYPE);
(void) GetIfPositions(RecordTypeNumbers, dp, ROOM_POSITIONS, ROOM_POSITION_BPF_REC_TYPE);
(void) GetIfPositions(RecordTypeNumbers, dp, ARENA_POSITIONS, ARENA_POSITION_BPF_REC_TYPE);
return (NumEEGChannels);
}
si4 GetEEGChannelList(ui4 *BPFEEGChannels, si1 *SectionStart){
ui1 *dp;
si1 KeyWordString[4196], Encoding[4196], **Parameters;
si4 NumberOfParameters, Type, i;
dp = (ui1 *)SectionStart;
if(!FindKeyWord((si1 *)dp, LIST_OF_EEG_CHANNELS, KeyWordString, Encoding)){
fprintf(MyStdErr,"\nCan't find keyWord (%s).\n", LIST_OF_EEG_CHANNELS);
return -1;
}
Type = GetType(KeyWordString);
switch(Type){
case -1:
fprintf(MyStdErr,"KeyWord.integer format not respected in: %s\nWrong File Format.\n", KeyWordString);
return -1;
case 0: // ASCII encoded
NumberOfParameters = GetASCIIEncoding(Encoding, &Parameters);
if(NumberOfParameters){
for(i=0; i < NumberOfParameters; i++){
BPFEEGChannels[i] = atoi(Parameters[i]) - 1; // The header count begins at 1 but the channel count begins at 0 for processing
}
}
return NumberOfParameters;
default:
fprintf(MyStdErr,"%s must be ASCII encoded\nWrong File Format.\n", KeyWordString);
return 0;
}
}
ui4 GetNumberOfChannels(BPFRecordTypeNumbers, SectionStart, KeyWord, RecType)
ui4 *BPFRecordTypeNumbers;
ui1 *SectionStart;
si1* KeyWord;
ui1 RecType;
{
ui1 *dp;
si1 KeyWordString[4196], Encoding[4196], **Parameters;
si4 Type, NumberOfParameters;
dp = SectionStart;
if(!FindKeyWord((si1 *)dp, KeyWord, KeyWordString, Encoding)){
fprintf(MyStdErr,"\nCan't find keyWord (%s).\n", KeyWord);
return(-1);
}
Type = GetType(KeyWordString);
switch(Type){
case -1:
fprintf(MyStdErr,"KeyWord.integer format not respected in: %s\nWrong File Format.\n", KeyWordString);
return(0);
case 0: // ASCII encoded
NumberOfParameters = GetASCIIEncoding(Encoding, &Parameters);
if(NumberOfParameters != 1){
fprintf(MyStdErr,"(%s).\nWrong File Format.\n", KeyWord);
return 0;
}
sscanf(Parameters[0],"%d", &(BPFRecordTypeNumbers[RecType]));
return(1);
default:
fprintf(MyStdErr,"%s must be ASCII encoded\nWrong File Format.\n", KeyWord);
return(0);
}
}
ui4 GetIfPositions(BPFRecordTypeNumbers, SectionStart, KeyWord, RecType)
ui4 *BPFRecordTypeNumbers;
ui1 *SectionStart;
si1* KeyWord;
ui1 RecType;
{
ui1 *dp;
si1 KeyWordString[4196], Encoding[4196], **Parameters;
si4 Type, NumberOfParameters;
dp = SectionStart;
if(!FindKeyWord((si1 *)dp, KeyWord, KeyWordString, Encoding)){
// fprintf(MyStdErr,"Can't find keyWord (%s). This is not critical.\n", KeyWord);
BPFRecordTypeNumbers[RecType] = 0;
return 1;
}
Type = GetType(KeyWordString);
switch(Type){
case -1:
fprintf(MyStdErr,"KeyWord.integer format not respected in: %s\nWrong File Format.\n", KeyWordString);
return(0);
case 0: // ASCII encoded
NumberOfParameters = GetASCIIEncoding(Encoding, &Parameters);
if(NumberOfParameters != 1){
fprintf(MyStdErr,"(%s).\nWrong File Format.\n", KeyWord);
return 0;
}
sscanf(Parameters[0],"%d", &(BPFRecordTypeNumbers[RecType]));
return(1);
default:
fprintf(MyStdErr,"%s must be ASCII encoded\nWrong File Format.\n", KeyWord);
return(0);
}
}
int GetGainsList(si2 *Gains, ui1 *SectionStart){
ui1 *dp;
si1 KeyWordString[4196], Encoding[4196], **Parameters;
si4 NumberOfParameters, Type, i;
dp = SectionStart;
if(!FindKeyWord((si1 *)dp, LIST_OF_GAINS, KeyWordString, Encoding)){
fprintf(MyStdErr,"\nCan't find keyWord (%s).\n", LIST_OF_GAINS);
return 1;
}
Type = GetType(KeyWordString);
switch(Type){
case -1:
fprintf(MyStdErr,"KeyWord.integer format not respected in: %s\nWrong File Format.\n", KeyWordString);
return 1;
case 0: // ASCII encoded
NumberOfParameters = GetASCIIEncoding(Encoding, &Parameters);
if(NumberOfParameters){
if(NumberOfParameters>MAX_NUMBER_OF_BPF_CHANNELS){
fprintf(MyStdErr,"GetGainsList ERRA: NumberOfParameters=%d > MAX_NUMBER_OF_BPF_CHANNELS=%d!!!\n",NumberOfParameters,MAX_NUMBER_OF_BPF_CHANNELS);
return 0;
}
// Gains = (ui2 *)calloc(MAX_NUMBER_OF_BPF_CHANNELS, sizeof(si2));
for(i=0; i < NumberOfParameters; i++){
// sscanf(Parameters[i],"%d", &(Gains[i]));
Gains[i] = atoi(Parameters[i]);
}
}
return 1;
default:
fprintf(MyStdErr,"%s must be ASCII encoded\nWrong File Format.\n", KeyWordString);
return 1;
}
return 1;
}
void bpf_instruct(){
printf("\n\nCall with a brain potential file (.bpf) format file in the pwd.\n");
printf("Outputs the EEG waveform for each channel in the BPF file to a .wfm file in the pwd.\n");
printf("The output voltage is in mV units.\n");
printf("Options:\n");
printf("-e# give the number of samples in the EEG record. Default is %d\n", NUMBER_OF_BPF_EEG_SAMPLES);
printf("-g# give the factor by which to multiply the voltages so as adjust the gain of the recorded data. The default is 1.0\n");
printf("\t\tUse this option for example if there is external gain that the recording system does not register in the BPF file header.\n");
printf("\t\tOne example is if an external amplifier was used and set to 1000x. In this case you will need to set option g to 0.001.\n");
printf("\t\tAnother example is the wireless DT recording system which has a hidden gain of 2.08\n");
}
#define PREAMP_GAIN 100.0
#define DSP_GAIN 20.0
int optind;
char *optarg;
ENDVERBATIM
VERBATIM
inline ui4 reverse_ui4(ui1* b)
{
ui1 *f;
static ui4 u;
f = (ui1 *) &u;
// b += 3;
//
// *f++ = *b--;
// *f++ = *b--;
// *f++ = *b--;
// *f = *b;
*f++ = *b++;
*f++ = *b++;
*f++ = *b++;
*f = *b;
return(u);
}
inline ui2 reverse_ui2(ui1* b)
{
ui1 *f;
static ui2 u;
f = (ui1 *) &u;
// b += 1;
//
// *f++ = *b--;
// *f = *b;
*f++ = *b++;
*f = *b;
return(u);
}
inline si2 reverse_si2(ui1* b)
{
ui1 *f;
static si2 u;
f = (ui1 *) &u;
// b += 1;
//
// *f++ = *b--;
// *f = *b;
*f++ = *b++;
*f = *b;
return(u);
}
//* probenumvec.readbpfunits("filename.bpf",outveclist,channumvec[,waveformflag])
static double readbpfunits (void* vv) {
Object* pList; void *vw;
int iMaxCols, iSz, npvec, ncvec, vlen[VRRY], ix, ix2, waveformflag;
char* fname; //input bpf file name
double *pvec, *cvec, retval, cl, pr, *vvo[VRRY];
FILE *fp;
si1 d, type, dp[MAX_BPF_REC_SIZE], *hp;
si4 ReadBPFHeader();
ui4 data_offset, RecordSizes[256];
ui1 dummy[5], key;
si1 line[256], string[256];
ui1 probe, clust, x, y;
ui2 ang, reverse_ui2();
si2 reverse_si2();
ui4 i, j, time_stamp, reverse_ui4();
si2 *ADCValue1, *ADCValue2, Gains[MAX_NUMBER_OF_BPF_CHANNELS];
si4 channel, BytesPerChannel, TetrodeBPFRecSize, ChannelsPerTetrode, SamplesPerWaveform;
sf8 ddt, volts1, volts2, energy;
si4 AcquisitionSystem;
ui4 EEGChList[MAX_NUMBER_OF_BPF_CHANNELS];
ui4 BPFRecordTypeNumbers[256];
ChannelsPerTetrode=4; SamplesPerWaveform=NUMBER_OF_BPF_TETRODE_SAMPLES;
npvec = vector_instance_px(vv,&pvec);
fname = gargstr(1);//input bpf filename
pList = *hoc_objgetarg(2); //output list of vectors
ncvec = vector_arg_px(3, &cvec);
waveformflag=(ifarg(4)?(int)*getarg(4):1);
iMaxCols = ivoc_list_count(pList);
if (iMaxCols>VRRY) {printf("readbpfunits; Can only handle VRRY==%d vecs\n",VRRY); hxe();}
retval=0.;
if (npvec!=ncvec) {
printf("pv.readbpfunits(file,outvecl,cv) ERRA: need probe#s in pv, clust#s in cv %d %d\n",\
npvec,ncvec); hxe();}
if (waveformflag) {
if (npvec*ChannelsPerTetrode!=iMaxCols) {
printf("readbpfunits ERRA0: waveform outveclist should be size %d * %d (#chans), not %d\n",\
ChannelsPerTetrode,npvec,iMaxCols); hxe();}
} else if (npvec!=iMaxCols) if (npvec*ChannelsPerTetrode!=iMaxCols) {
printf("readbpfunits ERRA0: outveclist should be size %d (#chans), not %d\n",\
npvec,iMaxCols); hxe();
}
//set pointers to list vectors
for(i=0;i<iMaxCols;i++) {
vlen[i]=0;
if (i==0) iSz=list_vector_px3(pList,i,&vvo[i],&vw); else {
if ((j=list_vector_px(pList,i,&vvo[i]))!=iSz) {
fprintf(stderr,"readbpfunits ERRF: different size vectors %d %d %d\n",i,j,iSz);hxe(); }
}
}
// set these here so that they can be determined by optional values
RecordSizes[EEG_BPF_REC_TYPE] = EEG_BPF_REC_SIZE;
RecordSizes[SINGLE_BPF_REC_TYPE] = SINGLE_BPF_REC_SIZE;
RecordSizes[STEREO_BPF_REC_TYPE] = STEREO_BPF_REC_SIZE;
RecordSizes[TETRODE_BPF_REC_TYPE] = TETRODE_BPF_REC_SIZE;
RecordSizes[SYNC_BPF_REC_TYPE] = SYNC_BPF_REC_SIZE;
RecordSizes[ROOM_POSITION_BPF_REC_TYPE] = ROOM_POSITION_BPF_REC_SIZE;
RecordSizes[ARENA_POSITION_BPF_REC_TYPE] = ARENA_POSITION_BPF_REC_SIZE;
RecordSizes[KEY_EVENT_BPF_REC_TYPE] = KEY_EVENT_BPF_REC_SIZE;
RecordSizes[INPUT_EVENT_BPF_REC_TYPE] = INPUT_EVENT_BPF_REC_SIZE;
RecordSizes[OUTPUT_EVENT_BPF_REC_TYPE] = OUTPUT_EVENT_BPF_REC_SIZE;
fp = fopen(fname, "r");
if(fp == NULL){printf("readbpfunits ERRZ: can't open %s for reading!\n", fname); goto dofree; }
while(fgets(line, 256, fp) != NULL){ sscanf(line, "%s", string);
if(!strcmp(string, HEADER_END_LABEL)) {data_offset=ftell(fp); break;}
}
hp = (si1 *)calloc(data_offset, sizeof(si1));
if(hp == NULL){fprintf(stderr,"readbpfunits ERRM: calloc failed\n"); goto dofree;}
rewind(fp);
if(fread(hp,sizeof(ui1), data_offset, fp) != data_offset){
fprintf(stderr,"readbpfunits ERRN: fread header failed\n"); goto dofree; }
fflush(stdout);
ReadBPFHeader(hp, RecordSizes, Gains, &AcquisitionSystem, EEGChList, BPFRecordTypeNumbers);
BytesPerChannel = 2 * SamplesPerWaveform;
ddt = 1.0 / (sf8)SamplesPerWaveform;
while(1) {
fseek(fp, data_offset,0);// do because in Windows fp doesn't increment properly after fread
if (fread(((void *)dp),sizeof(ui1), (size_t)1, fp) == EOF) break;
type = (si1)*dp;
data_offset += RecordSizes[type];
if(fread((void *)(dp+1),sizeof(ui1),(size_t)RecordSizes[type]-1,fp) != (RecordSizes[type]-1)){
fprintf(stderr,"Couldn't read complete record of type %c\n", type);
goto resize;
}
time_stamp = reverse_ui4(dp+1);
// printf("%c\t%d", type, time_stamp);
switch (type) {
case EEG_BPF_REC_TYPE: // ignore these
break;
case SYNC_BPF_REC_TYPE:
break;
case INPUT_EVENT_BPF_REC_TYPE:
break;
case OUTPUT_EVENT_BPF_REC_TYPE:
break;
case KEY_EVENT_BPF_REC_TYPE:
key = (si1) *(dp + BPF_KEY_EVENT_REC_OFFSET);
// printf("\t%c\n",key);
break;
case ROOM_POSITION_BPF_REC_TYPE:
case ARENA_POSITION_BPF_REC_TYPE:
// x = *(dp + BPF_POS_REC_X_OFFSET);
// y = *(dp + BPF_POS_REC_Y_OFFSET);
// ang = reverse_ui2(dp + BPF_POS_REC_ANG_OFFSET);
// printf("\t%d\t%d\t%d",x,y,ang);
break;
case SINGLE_BPF_REC_TYPE:
case STEREO_BPF_REC_TYPE:
case TETRODE_BPF_REC_TYPE:
probe = *(dp + BPF_RECORD_PROBE_OFFSET);
clust = *(dp + BPF_SPK_REC_CLUST_OFFSET);
pr=(double)probe; cl=(double)clust;
// if vec contains probe and in contains clust then save it to that vec
for (i=0.;i<npvec;i++) if (pvec[i]==pr&&cvec[i]==cl) break;
if ((ix=i)<npvec) { // add to the chosen vector
if (!waveformflag) {
if (vlen[ix]>=iSz) { // lengthen all vectors
iSz=((iSz==0)?1e4:iSz*2);
for (i=0;i<iMaxCols;i++) vvo[i]=list_vector_resize(pList,i,iSz);
}
vvo[ix][vlen[ix]]=(double)time_stamp/10.0;
vlen[ix]++;
} else {
ix*=ChannelsPerTetrode;
if (vlen[ix]+SamplesPerWaveform+10>=iSz) { // lengthen all vectors
iSz=((iSz==0)?1e4:iSz*2)+SamplesPerWaveform+10;
for (i=0;i<iMaxCols;i++) vvo[i]=list_vector_resize(pList,i,iSz);
}
for (channel=0; channel<ChannelsPerTetrode; channel++) {
ix2=ix+channel;
for(j=0; j<SamplesPerWaveform; j++) {
ADCValue2 = (si2 *)(dp + 7 + (2 * j) + (channel * BytesPerChannel));
*ADCValue2 = (si2)reverse_si2(ADCValue2);
// convert ADC value to voltage -- REMOVED
volts2 = 10.0 * (sf8)*ADCValue2 / (sf8)SHRT_MAX / (sf8)Gains[channel];
vvo[ix2][vlen[ix2]]=volts2;
vlen[ix2]++;
}
}
}
}
// printf("\t%d\t%d\t%0.4lf",probe, clust, energy);
break;
default:
printf("%c UNKNOWN RECORD TYPE", type); fflush(stdout);
break;
}
// printf("\n"); fflush(stdout);
}
resize:
for (i=0;i<iMaxCols;i++) list_vector_resize(pList,i,vlen[i]);
retval=1.0; //success
dofree: // free memory
if(VERBOSE>1) printf("readbpfunits: freeing memory\n");
if(fp) fclose(fp);
if(hp) free(hp);
return retval;
}
void SetBPFRecSizes(ui4 RecordSizes[256]) {
// set these here so that they can be determined by optional values
RecordSizes[EEG_BPF_REC_TYPE] = EEG_BPF_REC_SIZE;
RecordSizes[SINGLE_BPF_REC_TYPE] = SINGLE_BPF_REC_SIZE;
RecordSizes[STEREO_BPF_REC_TYPE] = STEREO_BPF_REC_SIZE;
RecordSizes[TETRODE_BPF_REC_TYPE] = TETRODE_BPF_REC_SIZE;
RecordSizes[SYNC_BPF_REC_TYPE] = SYNC_BPF_REC_SIZE;
RecordSizes[ROOM_POSITION_BPF_REC_TYPE] = ROOM_POSITION_BPF_REC_SIZE;
RecordSizes[ARENA_POSITION_BPF_REC_TYPE] = ARENA_POSITION_BPF_REC_SIZE;
RecordSizes[KEY_EVENT_BPF_REC_TYPE] = KEY_EVENT_BPF_REC_SIZE;
RecordSizes[INPUT_EVENT_BPF_REC_TYPE] = INPUT_EVENT_BPF_REC_SIZE;
RecordSizes[OUTPUT_EVENT_BPF_REC_TYPE] = OUTPUT_EVENT_BPF_REC_SIZE;
}
//* nq.fcdv.rdbpfu("filename.bpf",nq.vl)
// nq cols should be eg "timestamp", "probenum", "clustnum", "waveform" with vdec("waveform")
static double rdbpfu (void* vv) {
Object* pList; void *vw;
int iMaxCols, iSz, novec, cnt, vpr;
char* fname; //input bpf file name
double *ovec, retval, *vvo[4];
FILE *fp;
si1 d, type, dp[MAX_BPF_REC_SIZE], *hp;
si4 ReadBPFHeader();
ui4 data_offset, RecordSizes[256], filebytes;
ui1 dummy[5], key;
si1 line[256], string[256];
ui1 probe, clust, x, y;
ui2 ang, reverse_ui2();
si2 reverse_si2();
ui4 i, j, time_stamp, reverse_ui4();
si2 *ADCValue1, *ADCValue2, Gains[MAX_NUMBER_OF_BPF_CHANNELS], *wavep;
si4 chan, BytesPerChannel, TetrodeBPFRecSize, ChannelsPerTetrode, SamplesPerWaveform;
sf8 ddt, volts1, volts2, energy;
si4 AcquisitionSystem;
ui4 EEGChList[MAX_NUMBER_OF_BPF_CHANNELS];
ui4 BPFRecordTypeNumbers[256];
ChannelsPerTetrode=4; SamplesPerWaveform=NUMBER_OF_BPF_TETRODE_SAMPLES;
novec = vector_instance_px(vv,&ovec);
fname = gargstr(1);//input bpf filename
pList = *hoc_objgetarg(2); //output list of vectors
iMaxCols = ivoc_list_count(pList);
if (iMaxCols!=4) {printf("rdbpfu: NQS veclist should be size 4\n"); hxe();}
retval=0.;
//set pointers to list vectors
for(i=0;i<iMaxCols;i++) {
if (i==0) iSz=list_vector_px3(pList,i,&vvo[i],&vw); else {
if ((j=list_vector_px3(pList,i,&vvo[i],&vw))!=iSz) {
fprintf(stderr,"rdbpfu ERRF: different size vectors %d %d %d\n",i,j,iSz);hxe(); }
}
}
SetBPFRecSizes(RecordSizes);// set these here so that they can be determined by optional values
if((fp=fopen(fname, "r"))==NULL){printf("rdbpfu ERRZ: can't open %s for reading!\n", fname); goto dofree; }
while(fgets(line, 256, fp) != NULL){ sscanf(line, "%s", string);
if(!strcmp(string, HEADER_END_LABEL)) {data_offset=ftell(fp); break;}
}
hp = (si1 *)calloc(data_offset, sizeof(si1));
if(hp == NULL){fprintf(stderr,"rdbpfu ERRM: calloc failed\n"); goto dofree;}
fseek(fp,0,SEEK_END); // go to end of file
filebytes = ftell(fp); if(VERBOSE>1) printf("filebytes=%d\n",filebytes); // get the # of bytes
rewind(fp); // go back to beginning of file
if(fread(hp,sizeof(ui1), data_offset, fp) != data_offset){ // read header into buffer
fprintf(stderr,"rdbpfu ERRN: fread header failed\n"); goto dofree; }
fflush(stdout);
ReadBPFHeader(hp, RecordSizes, Gains, &AcquisitionSystem, EEGChList, BPFRecordTypeNumbers);//parse header
BytesPerChannel = 2 * SamplesPerWaveform;
ddt = 1.0/(sf8)SamplesPerWaveform;
cnt=vpr=0;
while (1) { //read structures from file
fseek(fp, data_offset,0);// do because in Windows fp doesn't increment properly after fread (??)
if (fread(((void *)dp),sizeof(ui1), (size_t)1, fp) == EOF) break;
type=(si1)*dp;
data_offset += RecordSizes[type]; // increment pointer
if(data_offset >= filebytes) { // check for EOF here before reading next record, which could be past end of file
if(VERBOSE>1) printf("data_offset=%d >= filebytes=%d\n",data_offset,filebytes);
break;
}
if(fread((void *)(dp+1),sizeof(ui1),(size_t)RecordSizes[type]-1,fp)!=(RecordSizes[type]-1)){
fprintf(stderr,"rdbpfu: Couldn't read complete record of type %c, @ offset=%d\n", type,data_offset);
goto resize;
}
switch (type) {
case INPUT_EVENT_BPF_REC_TYPE: // ignore these
case OUTPUT_EVENT_BPF_REC_TYPE:
case EEG_BPF_REC_TYPE:
case SYNC_BPF_REC_TYPE:
case ROOM_POSITION_BPF_REC_TYPE:
case ARENA_POSITION_BPF_REC_TYPE:
case KEY_EVENT_BPF_REC_TYPE:
break;
case SINGLE_BPF_REC_TYPE:
case STEREO_BPF_REC_TYPE:
case TETRODE_BPF_REC_TYPE:
time_stamp = reverse_ui4(dp+1);
probe = *(dp + BPF_RECORD_PROBE_OFFSET);
clust = *(dp + BPF_SPK_REC_CLUST_OFFSET);
// if vec contains probe and in contains clust then save it to that vec
if (cnt>=iSz) { // lengthen all vectors
iSz=((iSz==0)?1e4:iSz*2);
for (i=0;i<iMaxCols;i++) vvo[i]=list_vector_resize(pList,i,iSz);
}
vvo[0][cnt]=(double)time_stamp/10.0;
vvo[1][cnt]=(double)probe;
vvo[2][cnt]=(double)clust;
vvo[3][cnt]=(double)vpr;
cnt++;
if (vpr+ChannelsPerTetrode*SamplesPerWaveform>=novec) {
novec=(novec==0?10000:10*novec);
ovec=vector_newsize(vv,novec);
}
wavep = (si2*) (dp + 7);//temporary si2 pointer to waveform
for (chan=0;chan<ChannelsPerTetrode;chan++) for(j=0;j<SamplesPerWaveform;j++,wavep++) {
ovec[vpr++] = (double) 10.0 * wavep[0] / SHRT_MAX;
}
// printf("\t%d\t%d\t%0.4lf",probe, clust, energy);
break;
default:
printf("%c UNKNOWN RECORD TYPE", type); fflush(stdout);
break;
}
}
resize:
if(VERBOSE>1) printf("vpr = %d, novec = %d\n",vpr,novec);
for (i=0;i<iMaxCols;i++) list_vector_resize(pList,i,cnt);
vector_resize(vv,vpr);
retval=1.0; //success
dofree: // free memory
if(VERBOSE>1) printf("rdbpfu: freeing memory\n");
if(fp) fclose(fp);
if(hp) free(hp);
return retval;
}
// vector.readbpfeeg("filename.bpf",list_of_output_vectors [,records_per_sample,gain_adjust])
static double readbpfeeg(void* vv) {
Object* pList = 0;
int iMaxCols = 0;
double** vvo = 0;
double retval = 0.0;
FILE *fp = 0;
si1 file[4196], *c;
si1 type, dp[MAX_BPF_REC_SIZE], *hp = 0;
si4 d, channel, EEGBPFRecSize, EEGChannels = NUMBER_OF_BPF_EEG_CHANNELS, SamplesPerRecord = NUMBER_OF_BPF_EEG_SAMPLES;
si4 BytesPerChannel;
ui4 data_offset, RecordSizes[256];
si1 line[4096], string[4096];
ui1 probe, clust, x, y;
ui2 ang;
ui4 i, time_stamp;
si2 *ADCValue, Gains[MAX_NUMBER_OF_BPF_CHANNELS];
sf4 volts, AdjustGain = 1.0;
si4 AcquisitionSystem;
ui4 EEGChList[MAX_NUMBER_OF_BPF_CHANNELS];
ui4 BPFRecordTypeNumbers[256];
char* fname; //input bpf file name
double* pThisV=0; //stores eeg channel #s
int iThisVSz=0 , iSz=0, iOffsetStart=0,iOffsetCur=0;
iThisVSz = vector_instance_px(vv,&pThisV);
// set these here so that they can be determined by optional values
RecordSizes[EEG_BPF_REC_TYPE] = EEG_BPF_REC_SIZE;
RecordSizes[SINGLE_BPF_REC_TYPE] = SINGLE_BPF_REC_SIZE;
RecordSizes[STEREO_BPF_REC_TYPE] = STEREO_BPF_REC_SIZE;
RecordSizes[TETRODE_BPF_REC_TYPE] = TETRODE_BPF_REC_SIZE;
RecordSizes[SYNC_BPF_REC_TYPE] = SYNC_BPF_REC_SIZE;
RecordSizes[ROOM_POSITION_BPF_REC_TYPE] = ROOM_POSITION_BPF_REC_SIZE;
RecordSizes[ARENA_POSITION_BPF_REC_TYPE] = ARENA_POSITION_BPF_REC_SIZE;
RecordSizes[INPUT_EVENT_BPF_REC_TYPE] = INPUT_EVENT_BPF_REC_SIZE;
RecordSizes[OUTPUT_EVENT_BPF_REC_TYPE] = OUTPUT_EVENT_BPF_REC_SIZE;
fname = gargstr(1);//input bpf filename
pList = *hoc_objgetarg(2); //output list of vectors
iMaxCols = ivoc_list_count(pList); //output list count (should be >= EEGChannels+1)
if(ifarg(3)) SamplesPerRecord = *getarg(3); //optional argument for SamplesPerRecord
if(ifarg(4)) AdjustGain = (sf4) *getarg(4); //optional argument to adjust gain
fp = fopen(fname, "r");
if(fp == NULL){
printf("readbpfeeg ERRZ: can't open %s for reading!\n", fname);
goto dofree;
}
while(get_line(fp, line, 4096) != NULL){
sscanf(line, "%s", string);
if(!strcmp(string, HEADER_END_LABEL)){
data_offset = ftell(fp);
break;
}
}
hp = (si1 *)calloc(data_offset, sizeof(si1));
if(hp == NULL){
fprintf(stderr,"readbpfeeg ERRM: calloc failed\n");
goto dofree;
}
rewind(fp);
if(fread(hp,sizeof(ui1), data_offset, fp) != data_offset){
fprintf(stderr,"readbpfeeg ERRN: fread header failed\n");
goto dofree;
}
fflush(stdout);
//something wrong here
EEGChannels = ReadBPFHeader(hp, RecordSizes, Gains, &AcquisitionSystem, EEGChList, BPFRecordTypeNumbers);
if(EEGChannels < 1){
printf("readbpfeeg ERRB: No EEG Channels.\n");
goto dofree;
} else if(VERBOSE) printf("readbpfeeg: EEGChannels=%d\n",EEGChannels);
EEGBPFRecSize = EEG_BPF_REC_INFO_SIZE + (SamplesPerRecord * 2 * EEGChannels);
if(iMaxCols<EEGChannels+1){
fprintf(stderr,"readbpfeeg ERRD: need EEGChannels+1=%d for eegdata+time, List.count=%d!!\n",EEGChannels+1,iMaxCols);
goto dofree;
} else if(VERBOSE) printf("readbpfeeg: iMaxCols=%d\n",iMaxCols);
if(iThisVSz<EEGChannels){
printf("readbpfeeg ERRG: This vector sz %d < EEGChannels %d!\n",iThisVSz,EEGChannels);
goto dofree;
} else if(VERBOSE) printf("readbpfeeg: iThisVSz=%d\n",iThisVSz);
vvo = (double**)malloc(iMaxCols*sizeof(double*)); //alloc mem for pointers to list vectors
if(!vvo){
fprintf(stderr,"readbpfeeg ERRE: couldn't alloc mem for vvo!\n");
goto dofree;
}
//set pointers to list vectors
for(i=0;i<EEGChannels;i++){
pThisV[i]=1+EEGChList[i]; //store channel #s
if((iSz=list_vector_px(pList,i,&vvo[i]))<SamplesPerRecord){
fprintf(stderr,"readbpfeeg ERRF: list vectors size %d < BPF SamplesPerRecord %d!\n",iSz,SamplesPerRecord);
goto dofree;
} else if(VERBOSE>1) printf("readbpfeeg: iSz %d = %d\n",i,iSz);
}
if((iSz=list_vector_px(pList,i,&vvo[i]))<SamplesPerRecord)
{ fprintf(stderr,"readbpfeeg ERRT: time Vector invalid!\n");
goto dofree;
}
BytesPerChannel = 2 * SamplesPerRecord;
if(VERBOSE>1) printf("readbpfeeg: BytesPerChannel=%d SamplesPerRecord=%d\n",BytesPerChannel,SamplesPerRecord);
i = 0;
clock_t startt,endt;
startt=clock();
if(AcquisitionSystem == ACX)
{ while(1)
{ if(fread(((void *)dp),sizeof(ui1), (size_t)1, fp) == EOF)
break;
type = (si1)*dp;
if(fread((void *)(dp+1),sizeof(ui1),(size_t)RecordSizes[type]-1,fp) != (RecordSizes[type]-1))
break;
if(type!=EEG_BPF_REC_TYPE) continue;
for(channel=0; channel < EEGChannels; channel++)
{ iOffsetCur=iOffsetStart;
for(i=0; i < SamplesPerRecord; i++, iOffsetCur++)
{ ADCValue = (si2 *)(dp + 5 + (2 * i) + (channel * BytesPerChannel));
*ADCValue = (si2)reverse_si2(ADCValue);
// convert ADC value to voltage in mV
vvo[channel][iOffsetCur]=AdjustGain*10.*(sf4)*ADCValue/(sf4)SHRT_MAX/(sf4)Gains[EEGChList[channel]]*1e3;
}
}
iOffsetStart += SamplesPerRecord; //increment offset
}
}
else if(AcquisitionSystem == AXONA)
{ while(1)
{ if(fread(((void *)dp),sizeof(ui1), (size_t)1, fp) == EOF)
break;
type = (si1)*dp;
if(fread((void *)(dp+1),sizeof(ui1),(size_t)RecordSizes[type]-1,fp) != (RecordSizes[type]-1))
break;
if(type!=EEG_BPF_REC_TYPE) continue;
for(channel=0; channel < EEGChannels; channel++)
{ iOffsetCur=iOffsetStart;
for(i=0; i < SamplesPerRecord; i++,iOffsetCur++)
{ ADCValue = (si2 *)(dp + 5 + (2 * i) + (channel * BytesPerChannel));
*ADCValue = (si2)reverse_si2(ADCValue);
// convert ADC value to voltage in mV
vvo[channel][iOffsetCur]=1000.0*(AdjustGain*(sf4)*ADCValue/(sf4)SHRT_MAX);
}
}
iOffsetStart += SamplesPerRecord; //increment offset
}
}
int iT = EEGChannels, iLim = iOffsetStart; //create time index in minutes
double tcur = 0.0, tinc = (1/2e3)/60.0; //assumes sampling rate == 2Khz, which is default
for(i=0;i<iLim;i++)
{ vvo[iT][i] = tcur;
tcur += tinc;
}
endt=clock();
printf("main loop time = %gms\n", 1000.0*((double)(endt-startt))/CLOCKS_PER_SEC);
if(VERBOSE) printf("offset=%d\n",iOffsetStart);
retval=1.0;//success
//free memory
dofree:
return retval; // ???????????????????? why return before free ?????
if(VERBOSE>1) printf("readbpfeeg: freeing memory\n");
if(fp) fclose(fp);
if(vvo) free(vvo);
if(hp) free(hp);
return retval;
}
ENDVERBATIM
: bpfeegchlab(filepath,numchannels,str1,str2,...)
: return # of channel labels found (includes empty labels which are usually for SYNC channels)
FUNCTION bpfeegchlab () {
VERBATIM
double dRet;
dRet=-1.0;
FILE* fp = 0;
char* fname = 0, **plabs = 0, line[4096],string[4096],*pch;
fname = gargstr(1); //input bpf filename
fp = fopen(fname,"r");
if(!fp){
printf("bpfeegchlab ERRA: couldn't open %s for reading!\n",fname);
goto CHLCLEANUP;
}
int iChans = (int)*getarg(2), idx , jdx, iLine = 0, iMaxLine = 100 , iLabLen = strlen(CHANNEL_LABEL), iCH = 0;
if(iChans < 1){
printf("bpfeegchlab ERRB: invalid number of channels: %d!\n",iChans);
goto CHLCLEANUP;
}
plabs = (char**)calloc(iChans,sizeof(char*));
if(!plabs){
printf("bpfeegchlab ERRC: out of memory!\n");
goto CHLCLEANUP;
}
for(idx=0;idx<iChans;idx++){
plabs[idx] = gargstr(idx+3);
if(!plabs[idx]){
printf("bpfeegchlab ERRD: couldn't get input string %d!\n",idx+3);
goto CHLCLEANUP;
}
plabs[idx][0]=0; //set to empty string
}
while(iLine++<iMaxLine && iCH<iChans && get_line(fp,line,4096)!=NULL){
if(!strncmp(line,CHANNEL_LABEL,iLabLen)){
if(VERBOSE) printf("bpfeegchlab: %s\n",line);
pch=strchr(line,'\"');
if(!pch) continue;
pch++;
if(!*pch) continue;
idx=0;
while(*pch && *pch!='\"') plabs[iCH][idx++]=*pch++;
plabs[iCH][idx]=0; //terminate string with 0
if(VERBOSE) printf("bpfeegchlab: label %d = %s\n",iCH,plabs[iCH]);
iCH++;
} else if(VERBOSE>1) printf("bpfeegchlab: L%d=%s\n",iLine,line);
}
dRet=(double)iCH; //return # of channel labels found
CHLCLEANUP:
if(fp) fclose(fp);
if(plabs) free(plabs);
return dRet;
ENDVERBATIM
}
: get # of SYNC channels in a BPF file
FUNCTION bpfsyncchan () {
VERBATIM
FILE *fp = 0;
si1 file[4196], *c;
si1 type, dp[MAX_BPF_REC_SIZE], *hp = 0;
si4 d, channel, EEGBPFRecSize, EEGChannels = NUMBER_OF_BPF_EEG_CHANNELS;
ui4 data_offset, RecordSizes[256];
si1 line[4096], string[4096];
ui4 i;
si2 Gains[MAX_NUMBER_OF_BPF_CHANNELS];
sf4 AdjustGain = 1.0;
si4 AcquisitionSystem;
ui4 EEGChList[MAX_NUMBER_OF_BPF_CHANNELS];
ui4 BPFRecordTypeNumbers[256];
char* fname; //input bpf file name
// set these here so that they can be determined by optional values
RecordSizes[EEG_BPF_REC_TYPE] = EEG_BPF_REC_SIZE;
RecordSizes[SINGLE_BPF_REC_TYPE] = SINGLE_BPF_REC_SIZE;
RecordSizes[STEREO_BPF_REC_TYPE] = STEREO_BPF_REC_SIZE;
RecordSizes[TETRODE_BPF_REC_TYPE] = TETRODE_BPF_REC_SIZE;
RecordSizes[SYNC_BPF_REC_TYPE] = SYNC_BPF_REC_SIZE;
RecordSizes[ROOM_POSITION_BPF_REC_TYPE] = ROOM_POSITION_BPF_REC_SIZE;
RecordSizes[ARENA_POSITION_BPF_REC_TYPE] = ARENA_POSITION_BPF_REC_SIZE;
RecordSizes[INPUT_EVENT_BPF_REC_TYPE] = INPUT_EVENT_BPF_REC_SIZE;
RecordSizes[OUTPUT_EVENT_BPF_REC_TYPE] = OUTPUT_EVENT_BPF_REC_SIZE;
fname = gargstr(1);//input bpf filename
fp = fopen(fname, "r");
if(fp == NULL){
printf("bpfsyncchan ERRA: Can't open %s for reading\n", fname);
goto dofree;
}
while(get_line(fp, line, 4096) != NULL){
sscanf(line, "%s", string);
if(!strcmp(string, HEADER_END_LABEL)){
data_offset = ftell(fp);
break;
}
}
hp = (si1 *)calloc(data_offset, sizeof(si1));
if(hp == NULL){
fprintf(stderr,"bpfsyncchan ERRB: calloc failed\n");
goto dofree;
}
rewind(fp);
if(fread(hp,sizeof(ui1), data_offset, fp) != data_offset){
fprintf(stderr,"bpfsyncchan ERRC: fread header failed\n");
goto dofree;
}
fflush(stdout);
EEGChannels = ReadBPFHeader(hp, RecordSizes, Gains, &AcquisitionSystem, EEGChList, BPFRecordTypeNumbers); //something wrong here
if(VERBOSE) printf("bpfsyncchan : SYNCChannels=%d\n",BPFRecordTypeNumbers[SYNC_BPF_REC_TYPE]);
dofree:
if(VERBOSE>1) printf("bpfsyncchan: freeing memory\n");
if(fp) fclose(fp);
if(hp) free(hp);
return (double) (BPFRecordTypeNumbers[SYNC_BPF_REC_TYPE] >= 0 ? BPFRecordTypeNumbers[SYNC_BPF_REC_TYPE] : 0);
ENDVERBATIM
}
: get # of eeg channels in a BPF file
FUNCTION bpfeegchan(){
VERBATIM
FILE *fp = 0;
si1 file[4196], *c;
si1 type, dp[MAX_BPF_REC_SIZE], *hp = 0;
si4 d, channel, EEGBPFRecSize, EEGChannels = NUMBER_OF_BPF_EEG_CHANNELS;
ui4 data_offset, RecordSizes[256];
si1 line[4096], string[4096];
ui4 i;
si2 Gains[MAX_NUMBER_OF_BPF_CHANNELS];
sf4 AdjustGain = 1.0;
si4 AcquisitionSystem;
ui4 EEGChList[MAX_NUMBER_OF_BPF_CHANNELS];
ui4 BPFRecordTypeNumbers[256];
char* fname; //input bpf file name
// set these here so that they can be determined by optional values
RecordSizes[EEG_BPF_REC_TYPE] = EEG_BPF_REC_SIZE;
RecordSizes[SINGLE_BPF_REC_TYPE] = SINGLE_BPF_REC_SIZE;
RecordSizes[STEREO_BPF_REC_TYPE] = STEREO_BPF_REC_SIZE;
RecordSizes[TETRODE_BPF_REC_TYPE] = TETRODE_BPF_REC_SIZE;
RecordSizes[SYNC_BPF_REC_TYPE] = SYNC_BPF_REC_SIZE;
RecordSizes[ROOM_POSITION_BPF_REC_TYPE] = ROOM_POSITION_BPF_REC_SIZE;
RecordSizes[ARENA_POSITION_BPF_REC_TYPE] = ARENA_POSITION_BPF_REC_SIZE;
RecordSizes[INPUT_EVENT_BPF_REC_TYPE] = INPUT_EVENT_BPF_REC_SIZE;
RecordSizes[OUTPUT_EVENT_BPF_REC_TYPE] = OUTPUT_EVENT_BPF_REC_SIZE;
fname = gargstr(1);//input bpf filename
fp = fopen(fname, "r");
if(fp == NULL){
printf("bpfeegchan ERRA: Can't open %s for reading\n", fname);
goto dofree;
}
while(get_line(fp, line, 4096) != NULL){
sscanf(line, "%s", string);
if(!strcmp(string, HEADER_END_LABEL)){
data_offset = ftell(fp);
break;
}
}
hp = (si1 *)calloc(data_offset, sizeof(si1));
if(hp == NULL){
fprintf(stderr,"bpfeegchan ERRB: calloc failed\n");
goto dofree;
}
rewind(fp);
if(fread(hp,sizeof(ui1), data_offset, fp) != data_offset){
fprintf(stderr,"bpfeegchan ERRC: fread header failed\n");
goto dofree;
}
fflush(stdout);
EEGChannels = ReadBPFHeader(hp, RecordSizes, Gains, &AcquisitionSystem, EEGChList, BPFRecordTypeNumbers); //something wrong here
if(EEGChannels < 1){
printf("bpfeegchan ERRD: No EEG Channels.\n");
goto dofree;
} if(VERBOSE) printf("bpfeegchan : EEGChannels=%d\n",EEGChannels);
dofree:
if(VERBOSE>1) printf("bpfeegchan: freeing memory\n");
if(fp) fclose(fp);
if(hp) free(hp);
return (double)EEGChannels;
ENDVERBATIM
}
: get total # of samples for a single EEG channel for the full duration of recording of a BPF file
FUNCTION bpfeegsamp(){
VERBATIM
FILE *fp = 0;
si1 file[4196], *c;
si1 type, dp[MAX_BPF_REC_SIZE], *hp = 0;
si4 d, channel, EEGBPFRecSize, EEGChannels = NUMBER_OF_BPF_EEG_CHANNELS;
//samples per channel for single unit of time - if 2000Hz recording, then SamplesPerRecord = 2000
si4 SamplesPerRecord = NUMBER_OF_BPF_EEG_SAMPLES;
si4 BytesPerChannel;
ui4 data_offset, RecordSizes[256];
si1 line[4096], string[4096];
ui1 probe, clust, x, y;
ui2 ang;
ui4 i;
si2 *ADCValue, Gains[MAX_NUMBER_OF_BPF_CHANNELS];
sf4 volts, **EEGdata = 0, AdjustGain = 1.0;
si4 AcquisitionSystem;
ui4 EEGChList[MAX_NUMBER_OF_BPF_CHANNELS];
ui4 BPFRecordTypeNumbers[256];
char* fname; //input bpf file name
int iTotalSamplesPerChannel=0; //samples per eeg channel for full BPF recording
// set these here so that they can be determined by optional values
RecordSizes[EEG_BPF_REC_TYPE] = EEG_BPF_REC_SIZE;
RecordSizes[SINGLE_BPF_REC_TYPE] = SINGLE_BPF_REC_SIZE;
RecordSizes[STEREO_BPF_REC_TYPE] = STEREO_BPF_REC_SIZE;
RecordSizes[TETRODE_BPF_REC_TYPE] = TETRODE_BPF_REC_SIZE;
RecordSizes[SYNC_BPF_REC_TYPE] = SYNC_BPF_REC_SIZE;
RecordSizes[ROOM_POSITION_BPF_REC_TYPE] = ROOM_POSITION_BPF_REC_SIZE;
RecordSizes[ARENA_POSITION_BPF_REC_TYPE] = ARENA_POSITION_BPF_REC_SIZE;
RecordSizes[INPUT_EVENT_BPF_REC_TYPE] = INPUT_EVENT_BPF_REC_SIZE;
RecordSizes[OUTPUT_EVENT_BPF_REC_TYPE] = OUTPUT_EVENT_BPF_REC_SIZE;
fname = gargstr(1);//input bpf filename
if(ifarg(2)) SamplesPerRecord = *getarg(2); //optional argument for SamplesPerRecord
fp = fopen(fname, "r");
if(fp == NULL)
{ printf("bpfeegsamp ERRA: Can't open %s for reading\n", fname);
goto dofree;
}
while(get_line(fp, line, 4096) != NULL)
{ sscanf(line, "%s", string);
if(!strcmp(string, HEADER_END_LABEL))
{ data_offset = ftell(fp);
break;
}
}
hp = (si1 *)calloc(data_offset, sizeof(si1));
if(hp == NULL)
{ fprintf(stderr,"bpfeegsamp ERRB: calloc failed\n");
goto dofree;
}
rewind(fp);
if(fread(hp,sizeof(ui1), data_offset, fp) != data_offset)
{ fprintf(stderr,"bpfeegsamp ERRC: fread header failed\n");
goto dofree;
}
fflush(stdout);
EEGChannels = ReadBPFHeader(hp, RecordSizes, Gains, &AcquisitionSystem, EEGChList, BPFRecordTypeNumbers); //something wrong here
if(EEGChannels < 1)
{ printf("bpfeegsamp ERRD: No EEG Channels.\n");
goto dofree;
} if(VERBOSE) printf("bpfeegsamp: EEGChannels=%d\n",EEGChannels);
EEGBPFRecSize = EEG_BPF_REC_INFO_SIZE + (SamplesPerRecord * 2 * EEGChannels);
BytesPerChannel = 2 * SamplesPerRecord;
i = 0;
while(1)
{
if(fread(((void *)dp),sizeof(ui1), (size_t)1, fp) == EOF)
break;
type = (si1)*dp;
data_offset += RecordSizes[type];
if(fread((void *)(dp+1),sizeof(ui1),(size_t)RecordSizes[type]-1,fp) != (RecordSizes[type]-1))
break;
if(type == EEG_BPF_REC_TYPE)
iTotalSamplesPerChannel+=SamplesPerRecord;
}
if(VERBOSE) printf("bpfeegsamp: TotalSamplesPerChannel=%d\n",iTotalSamplesPerChannel);
//free memory
dofree:
if(VERBOSE>1) printf("bpfeegsamp: freeing memory\n");
if(fp) fclose(fp);
if(hp) free(hp);
return (double)iTotalSamplesPerChannel;
ENDVERBATIM
}
VERBATIM
// * readbdat - read the extra/intra-cellular .dat files provided by buzsaki
// vector.readbdat("filename.dat",list_of_output_vectors,numchannels)
static double readbdat (void* vv) {
FILE* fp;
ListVec* pList;
double** pLV, retval;
int i,j,iChans,*pLen,iMinSz;
short* pbuf;
char* fname;
retval=0.0; fp=0; pList=0; pbuf=0;
fname=gargstr(1);
fp=fopen(fname,"rb");
if(!fp){
printf("readbdat ERRA: couldn't open %s for reading!\n",fname);
goto CLEANBD;
}
pList=AllocListVec(*hoc_objgetarg(2));
if(!pList){
printf("readbdat ERRB: couldn't initialize list vec arg 2!\n");
goto CLEANBD;
}
pLV = pList->pv; pLen = pList->plen;
iChans = (int)*getarg(3);
if(iChans < 1) {
printf("readbdat ERRC: num channels must be > 0!\n",iChans);
goto CLEANBD;
}
if(pList->isz < iChans) {
printf("readbdat ERRD: ListVec arg 2 size %d < # of channels %d!\n",pList->isz,iChans);
goto CLEANBD;
}
j=0;
pbuf = (short*)malloc(sizeof(short)*iChans);
iMinSz=pLen[0];
for(i=0;i<iChans;i++) if(pLen[i]<iMinSz) iMinSz=pLen[i];
while(1) {
if(fread(pbuf,sizeof(short), iChans, fp)!=iChans) break;
if(j < iMinSz) for(i=0;i<iChans;i++) pLV[i][j]=pbuf[i];
j++;
}
if(j > iMinSz) printf("readbdat WARNA: could only read %d samples, file has %d samples\n",iMinSz,j);
retval=j;
CLEANBD:
if(pList) FreeListVec(&pList); if(fp) fclose(fp); if(pbuf) free(pbuf);
return retval;
}
ENDVERBATIM
: get the probe number from a TS file-name -- convention is filename ends in
: P%dC%d , so the # after P is the probe number. probe,cluster are digits
FUNCTION GetProbeNum () {
VERBATIM
char* fname,buf[10]; int i,j,sz,ploc;
fname=gargstr(1);
sz=strlen(fname);
for(i=sz-1;i>=0;i--) if(fname[i]=='P') break;
if(i-1<=0 || fname[i-1] != '.') return -1;
ploc=i; buf[0]=j=0;
for(i=ploc+1;i<sz && fname[i]>='0' && fname[i]<='9' && j<10;i++) buf[j++]=fname[i];
if(i>=sz || fname[i]!='C') return -1;
if(j<10) buf[j]=0; else buf[9]=0;
return atoi(buf);
ENDVERBATIM
}
: get the cluster number from a TS file-name -- convention is filename ends in
: P%dC%d , so the # after C is the cluster number. probe,cluster are digits
FUNCTION GetClustNum () {
VERBATIM
char* fname,buf[10]; int i,j,sz,cloc;
fname=gargstr(1);
sz=strlen(fname);
for(i=sz-2;i>=0;i--) if(fname[i]=='C') break;
cloc=i; buf[0]=j=0;
if(i+1>=sz) return -1;
for(i=cloc+1;i<sz && fname[i]>='0' && fname[i]<='9' && j<10;i++) buf[j++]=fname[i];
if(j<10) buf[j]=0; else buf[9]=0;
return atoi(buf);
ENDVERBATIM
}
:* PROCEDURE install()
PROCEDURE install () {
if (INSTALLED==1) {
printf("$Id: place.mod,v 1.141 2011/07/06 15:24:40 samn Exp $\n")
} else {
INSTALLED=1
VERBATIM
install_vector_method("simts_PLACE", simts); // need to give _PLACE explicitly
install_vector_method("rd_PLACE", rd);
install_vector_method("rdtts_PLACE", rdtts);
install_vector_method("kendall",kendal2);
install_vector_method("kend2",kend2);
install_vector_method("kendall_PLACE", kendall);
install_vector_method("mkgaussfield_PLACE", mkgaussfield);
install_vector_method("dumpratemap_PLACE", dumpratemap);
install_vector_method("glob", glob); // not particularly placeish
install_vector_method("cumul", cumul); // not particularly placeish
install_vector_method("readbpfunits",readbpfunits);
install_vector_method("rdbpfu",rdbpfu);
install_vector_method("readbpfeeg",readbpfeeg);
install_vector_method("readbdat",readbdat);
ENDVERBATIM
}
}
