/* Last edited on 2024-01-04 14:20:48 by stolfi */ /* ---------------------------------------------------------------------- */ /* neuromat_filter_bandpass_gauss.c */ void neuromat_filter_bandpass_gauss_compute_parms ( double flo0, double flo1, double fhi1, double fhi0, double *sigma_lo_P, double *sigma_hi_P, double *mag_P, bool_t verbose ) { demand((!isnan(flo0)) && (!isnan(flo1)) && (!isnan(fhi1)) && (!isnan(fhi0)), "invalid NAN arguments"); demand(flo1 < fhi1, "cannot have {flo1 >= fhi1}"); demand(flo0 >= 0, "invalid {flo0}"); double tiny = 1.0e-8; double sigma_lo = NAN, sigma_hi = NAN, mag = NAN; if ((flo0 > 0) && (fhi0 < +INF)) { demand(flo0 < flo1, "invalid {flo0, flo1}"); demand(fhi1 < fhi0, "invalid {fhi1, fhi0}"); assert((flo1 < +INF) && (fhi1 > 0)); /* Since {flo1 < fhi1}. */ double fmid = (flo1 + fhi1)/2; /* Freq of approximate maximum. */ mag = 1.0; sigma_lo = (flo0 + flo1)/2; sigma_hi = (fhi0 + fhi1)/2; fprintf(stderr, " %s: sigma_lo = %16.12f sigma_hi = %16.12f\n", __FUNCTION__, sigma_lo, sigma_hi); for (int32_t pass = 0; pass < 3; pass++) { double Whi_flo0 = neuromat_filter_lowpass_gauss(flo0, sigma_hi); sigma_lo = neuromat_filter_lowpass_gauss_compute_sigma(flo0, Whi_flo0 - tiny/mag); if (verbose) { fprintf(stderr, " %s: Whi(flo0) = %16.12f sigma_lo = %16.12f\n", __FUNCTION__, Whi_flo0, sigma_lo); } double Wlo_fhi0 = neuromat_filter_lowpass_gauss(fhi0, sigma_lo); sigma_hi = neuromat_filter_lowpass_gauss_compute_sigma(fhi0, tiny/mag + Wlo_fhi0); if (verbose) { fprintf(stderr, " %s: Wlo(fhi0) = %16.12f sigma_hi = %16.12f\n", __FUNCTION__, Wlo_fhi0, sigma_hi); } double gmid = neuromat_filter_bandpass_gauss(fmid, sigma_lo, sigma_hi, 1.0); mag = 1.0/gmid; if (verbose) { fprintf(stderr, " %s: gmid = %16.12f mag = %16.2f\n", __FUNCTION__, gmid, mag); } } } else if (fhi0 < +INF) { demand(flo1 == 0, "{flo1} must be zero when {flo0} is zero"); sigma_lo = 0; demand(fhi1 < fhi0, "invalid {fhi1, fhi0}"); assert(fhi1 > 0); /* Since {flo1 < fhi1}. */ sigma_hi = neuromat_filter_lowpass_gauss_compute_sigma(fhi0, tiny); mag = 1.0; } else if (flo0 > 0) { demand(fhi1 == +INF, "{fhi1} must be {+INF} when {fhi0} is {+INF}"); sigma_hi = +INF; demand(flo0 < flo1, "invalid {fhi1, fhi0}"); assert(flo1 < +INF); /* Since {flo1 < fhi1}. */ sigma_lo = neuromat_filter_lowpass_gauss_compute_sigma(flo0, 1-tiny); mag = 1.0; } else { demand(flo1 == 0, "{flo1} must be zero when {flo0} is zero"); demand(fhi1 == +INF, "{fhi1} must be {+INF} when {fhi0} is {+INF}"); sigma_lo = 0; sigma_hi = +INF; mag = 1.0; } if (verbose) { fprintf(stderr, " %s: flo0 = %16.12f flo1 = %16.12f fhi1 = %16.12f fhi0 = %16.12f", __FUNCTION__, flo0, flo1, fhi1, fhi0); fprintf(stderr, " --> sigma_lo = %16.12f sigma_hi = %16.12f mag = %16.12f\n", sigma_lo, sigma_hi, mag); } (*sigma_lo_P) = sigma_lo; (*sigma_hi_P) = sigma_hi; (*mag_P) = mag; } demand(fsmp > 0, "invalid {fsmp}"); /* Reduce the parameter {zf} to {[-(zsmp/2) _ +(zsmp/2)]}: */ zf = zf - zsmp*floor(zf/zsmp + 0.5); double zsmp = fsmp/sigma; /* Range of summation for folded-over tails: */ int32_t M = (int32_t)ceil(9.0/zsmp); double gtf = 0, gt0 = 0; for (int32_t k = M; k > 0; k--) { double zkfp = zf + k*zsmp; double gkfp = exp(-0.5*zkfp*zkfp); double zkfm = zf - k*zsmp; double gkfm = exp(-0.5*zkfm*zkfm); double zk0 = k*zsmp; double gk0 = exp(-0.5*zk0*zk0); gtf += (gkfp+gkfm); gt0 += 2*gk0; } return (gf + gtf)/(1 + gt0); /* ---------------------------------------------------------------------- */ /* neuromat_filter_bandpass_log_gauss.c */ /* Range of summation for folded-over tails: */ int32_t M = (int32_t)ceil(9.0/zsmp); double gtf = 0, gt0 = 0; for (int32_t k = M; k > 0; k--) { double zkfp = zf + k*zsmp; double gkfp = exp(-0.5*zkfp*zkfp); double zkfm = zf - k*zsmp; double gkfm = exp(-0.5*zkfm*zkfm); double zk0 = k*zsmp; double gk0 = exp(-0.5*zk0*zk0); gtf += (gkfp+gkfm); gt0 += 2*gk0; } /* ---------------------------------------------------------------------- */ /* neuromat_median_filter.c */ auto int32_t *compx(const void *a, const void *b, void *data); /* Assumes {a,b} point to integers. Compares {x[*a]} with {x[*b]}. */ double neuromat_median_filter_compute_sample(int32_t nt, double x[], int32_t nwt, double wt[], int32_t kx[]) { /* Collect the indices {kx[0..nwt-1]} of the window samples, sorted by {x} value: */ for (int32_t j = -hwdi; j <= +hwdi; j++) { kx[j + hwdi] = it + j; } /* Sort {kx[0..nwt-1]} by increasing {x[kx[k]]}: */ /* Compute the weighted median: */ double sumw = 0; int32_t km = 0; while ((km < nwt) && (sumw < 0.5)) { sumw += wt[hwt + kx[km]-it]; km++; } return x[km]; } /* ---------------------------------------------------------------------- */ auto void skip_line(void); /* Consumes the remaining characters up to and including the EOL. Bombs out if EOF before EOL. */ /* Internal implementations: */ void skip_line(void) { int32_t r; do { r = fgetc(rd); } while ((r != EOF) && (r != '\n')); if (r == EOF) { fprintf(stderr, "** EOF found while skipping line %d\n", nl+1); exit(1); } } void neuromat_eeg_pca_remove_patterns ( int nt, int ne, double **vin, int np, double *P, double **coef, double **vot ); /* The procedure assumes that {coef[it][ip]} is a coefficient for pattern {ip} in frame {it}, for {it} in {0..nt-1} and {ip} in {0..np-1} The procedure assumes that {vot} is a C-style matrix with {nt} rows and at least {ne} columns. The program subtracts from each row {vi[0..ne-1] = vin[it][0..ne-1]} the fitted components, namely {ci*P} where {ci[0..np-1] = coef[it][0..np-1]}. It stores the residual into {vo[0..ne-1] = vot[it][0..ne-1]}. */ void neuromat_eeg_pca_remove_patterns ( int nt, int ne, double **vin, int np, double *P, double **coef, double **vot ) { demand(ne >= 0, "invalid {ne}"); demand(np >= 0, "invalid {np}"); demand(np <= ne, "more components than electrodes"); int it; for (it = 0; it < nt; it++) { double *vi = vin[it]; /* Original electrode values {vi[0..ne-1]}. */ double *cf = coef[it]; /* Coefficients of each pattern {cf[0..np-1]}. */ double *vo = vot[it]; /* Residual is {vo[0..ne-1]}. */ double vr[ne]; /* Reconstructed frame from fitted patterns. */ } } /* ---------------------------------------------------------------------- */ int neuromat_image_timeline_position(double tlo, double t, double thi, sign_t round, int xsz); /* Maps affinely {t} from {tlo__thi} to {0_xsz} and rounds in the direction {round}. */ float_image_t *neuromat_image_make_slider_bar(int NC, int NX, int NY, int *bxminP, int *bxmaxP) { demand(NY % 2 == 0, "slider bar image height must be multiple of 2"); demand(NY >= 4, "slider bar image height too small"); demand(NX > NY, "slider bar image too narrow"); /* Size and position of trigger axis bar (pixels): */ int bysize = 2; int bymin = (NY - bysize)/2; int bymax = NY - 1 - bymin; /* Paint the axis bar: */ int bxmin = bymin; int bxmax = NX - 1 - bxmin; float_image_t *img = float_image_new(NC, NX, NY); float_image_fill(img, 0.0); float vbar = 0.50; float_image_fill_rectangle(img, bxmin, bxmax, bymin, bymax, vbar); (*bxminP) = bxmin; (*bxmaxP) = bxmax; return img; } void neuromat_eeg_image_paint_marker_tics ( float_image_t *img, int bxmin, int bxmax, int mymin, int mymax, int nt, int nc, double **val, int ic_mark ); /* Assumes that that {val[0..nt-1][0..nc-1]} is a measurement set with {nc} channels sampled at {nt} equally spaced sampling times. Paints vertical tic marks into {img} at the positions that correesponds to the indices {it} when the marker channel {val[it][ic_mark]} is nonzero. The tics will span pixel rows {mymin..mymax} inclusive, will lie within columns {bxmin..bxmax} inclusive. Specifically, each tic mark is centered at the fractional abscissa {neuromat_image_slider_bar_position(bxmin, bxmax, it+0.5, nt)}. Tic marks are at least one pixel wide, and closely spaced tic marks (less than 1 pixel apart) may be coalesced. */ void neuromat_eeg_image_paint_marker_tics ( float_image_t *img, int bxmin, int bxmax, int mymin, int mymax, int nt, int nc, double **val, int ic_mark ) { auto void paint_marker(double xa, double xb); /* Paint a marker from fractional column indices {xa} to {xb}. */ /* Plot frames where marker is nonzero, merging any markers that are less than 1 pixel apart: */ double xup = -INF; /* Fractional abscissa of last start-of-mark. */ double xdn = -INF; /* Fractional abscissa of last end-of-mark. */ int it; for (it = 0; it < nt; it++) { double vmark = val[it][ic_mark]; demand(! isnan(vmark), "marker channel is NAN"); if (vmark != 0) { double xlo = neuromat_image_slider_bar_position(bxmin, bxmax, (double)it, nt); double xhi = neuromat_image_slider_bar_position(bxmin, bxmax, (double)it + 1, nt); if (xhi - xlo < 1.0) { /* Fatten marker to 1 pixel: */ double xmd = (xlo + xhi)/2; xlo = xmd - 0.5; xhi = xmd + 0.5; } if (xlo >= xdn + 1.0) { /* Paint accumulated previous marker, if any: */ if (xdn > xup) { paint_marker(xup, xdn); } /* Start new marker: */ xup = xlo; } /* Current marker ends here so far: */ xdn = xhi; } } /* Paint last marker, if any: */ if (xdn > xup) { paint_marker(xup, xdn); } void paint_marker(double xa, double xb) { int NC = (int)img->sz[0]; double hwd = 0; int c; for (c = 0; c < NC; c++) { float vmark = 0.80f; float_image_paint_rectangle(img, c, xa, xb, (double)mymin, (double)mymax+1, hwd, vmark, NAN, 6); } } } double neuromat_image_slider_bar_position(int bxmin, int bxmax, double t, int nt); /* Returns the fractional column index along the slider bar that corresponds to a particular time {t} in the range {[0 _ nt]}. Assumes that the slider bar spans pixel columns {bxmin..bxmax} inclusive. Namely, if {t=0} returns {bxmin}, and if {t=nt} returns {bxmax+1}. Other values of {t} are interpolated linearly. */ double neuromat_image_slider_bar_position(int bxmin, int bxmax, double t, int nt) { /* Linear interpolation along slider bar: */ demand(bxmin <= bxmax, "bad image column range"); double bxsize = bxmax + 1 - bxmin; return bxmin + (t*bxsize)/nt; } /* ---------------------------------------------------------------------- */ /* Find the two largest weights {wmax,wsec} and the index of {wmax}: */ double wmax = -INF; double wsec = -INF; int imax = -1; for (ie = 0; ie < ne; ie++) { double wi = wraw[ie]; assert (wi >= 0); if (wi > wmax) { wsec = wmax; wmax = wi; imax = ie; } else if (wi > wsec) { wsec = wi; } } assert(wsec > 0); /* Normalize the weights so that the second-largest is 1.0: */ for (ie = 0; ie < ne; ie++) { wraw[ie] /= wsec; } wmax = wmax/wsec; wsec = 1.0; int nb = (wmax < 1.0e-7 ? 3 : 2); /* Number of indeterminate coefficients. */ /* ---------------------------------------------------------------------- */ typedef enum { neuromat_eeg_image_basis_type_POLY, } neuromat_eeg_image_basis_type_t; switch(btype) { case neuromat_eeg_image_basis_type_SHEPARD: { neuromat_eeg_image_basis_fill_shepard(msk, ne, bas, ctr, rad, pos3D); } break; case neuromat_eeg_image_basis_type_POLY: { neuromat_eeg_image_basis_fill_poly(msk, ne, bas, ctr, rad, pos3D); } break; case neuromat_eeg_image_basis_type_RADIAL: { double rho = 2.0; neuromat_eeg_image_basis_fill_radial(msk, ne, bas, ctr, rad, rho, pos3D); } break; case neuromat_eeg_image_basis_type_VORONOI: { neuromat_eeg_image_basis_fill_voronoi(msk, ne, bas, ctr, rad, pos3D); } break; default: assert(FALSE); } void neuromat_eeg_image_basis_fill_poly(float_image_t *msk, int ne, float_image_t *bas[], r3_t pos3D[]); /* The basis consists of bivariate polynomials. It usually has sample values less than 0 and greater than 1. */ void neuromat_eeg_image_basis_fill_poly(int ne, float_image_t *msk, float_image_t *bas[], r2_t *ctr, r2_t *rad, r3_t pos3D) { bool_t verbose = FALSE; if (ne == 0) { return; } /* We start with a raw basis polynomials in {X,Y}, where {X,Y} range over {[-1 _ +1]} in the rectangle {ctr ± rad}. The basis consists of all monomials of degree up to {g} except the constant polynomial, shifted so that their integral over the rectangle of radius {rad} is zero. For degree {g=3} that gives exactly 20 basis elements. */ /* Find degree of polynomial {g}: */ int g = 0; /* Degree of polynomial. */ int nm = 0; /* Number of basis monomials. */ while (nm < ne) { g++; nm = (g + 1)*(g + 2)/2 - 1; } if (verbose) { fprintf(stderr, "using polynomial basis of degree %d\n", g); } if (verbose) { fprintf(stderr, "there are %d monomials not counting the unit one\n", nm); } assert(nm == ne); /* Otherwise we will need to optimize... */ /* Enumerate monomial basis: */ double eX[nm], eY[nm]; /* Exponents of X,Y in monomials. */ double C[nm]; /* Constant term in each monomial for zero integral. */ int kg; int km = 0; for (kg = 1; kg <= g; kg++) { int ix; for (ix = 0; ix <= kg; ix++) { int iy = kg - ix; /* Integrals of {X^ix} and {Y^iy} over {[-1 _ +1]}: */ double IX = ((ix % 2) == 0 ? 2.0/(ix + 1.0) : 0.0); double IY = ((iy % 2) == 0 ? 2.0/(iy + 1.0) : 0.0); double CXY = -IX*IY; /* Save the polynomial data: */ if (verbose) { fprintf(stderr, " M_%-3d = X^%d*Y^%d%+16.12f\n", km, ix, iy, CXY); } eX[km] = ix; eY[km] = iy; C[km] = CXY; km++; } } assert(km == nm); auto double monomial(int j, double X, double Y); /* evaluates monomial number {j} at point {X,Y}, assuming the schematic head is the unit disk. */ /* Now fill the images with the Lagrangian basis derived from the monomials: */ neuromat_eeg_image_basis_fill_lagrangian(ne, bas, ctr, rad, pos, monomial); return; /* Implementation of local functions: */ double monomial(int j, double X, double Y) { return pow(X, eX[j])*pow(Y, eY[j]) + C[j]; } } void neuromat_eeg_image_basis_fill_radial(float_image_t *msk, int ne, float_image_t *bas[], r2_t *ctr, r2_t *rad, r3_t pos3D, double rho) { demand(msk->sz[0] == 1, "mask must be monochromatic"); int NX = (int)msk->sz[1]; int NY = (int)msk->sz[2]; bool_t verbose = FALSE; if (ne == 0) { return; } /* We start with a basis of radial functions on the units sphere. Each element is centered at one electrode position and has a gauss times sinc shape with first root at the nearest neighbor. */ double cx = ctr->c[0]; double cy = ctr->c[1]; double rx = rad->c[0]; double ry = rad->c[1]; /* Compute normalized electrode coords {Xe[0..ne-1],Ye[0..ne-1]}: */ double Xe[ne], Ye[ne]; int ie; for (ie = 0; ie < ne; ie++) { /* Map electrode {ie} to unit disk schematic head: */ r2_t *posi = &(pos[ie]); Xe[ie] = (posi->c[0] - cx)/rx; Ye[ie] = (posi->c[1] - cy)/ry; } /* Define the scale radius {Re[0..ne-1]} from electrodes to nearest neighbors: */ double Re[ne]; for (ie = 0; ie < ne; ie++) { /* Set {Re[ie]} to distance to nearest neighbor: */ Re[ie] = +INF; int je; for (je = 0; je < ne; je++) { if (ie != je) { double dij = hypot(Xe[ie] - Xe[je], Ye[ie] - Ye[je]); if (dij < Re[ie]) { Re[ie] = dij; } } } if (verbose) { fprintf(stderr, "electrode C%03d at ( %+7.3f %+7.3f ) radius %8.4f\n", ie+1, Xe[ie], Ye[ie], Re[ie]); } } /* Now fill the images with the Lagrangian basis derived from the radial basis: */ auto double mexican(double d, double sigma, double tau); /* The mother function of the radial basis, with width parameter {sigma} and node spacing {tau}, evaluated at distance {d} from the center. */ auto double radial(int j, double X, double Y); /* The raw radial element centered at electrode {j}. */ neuromat_eeg_image_basis_fill_lagrangian(ne, bas, ctr, rad, pos, radial); return; /* Implementation of local functions: */ double radial(int j, double X, double Y) { double d = hypot(X - Xe[j], Y - Ye[j]); return mexican(d, rho*Re[j], Re[j]); } double mexican(double d, double sigma, double tau) { if (d == 0) { return 1.0; } double tbell = d/sigma; if (tbell >= 8.5) { return 0.0; } double bell = exp(-tbell*tbell/2); double tsinc = d/tau; double sinc = sin(M_PI*tsinc)/(M_PI*tsinc); return bell*sinc; } } /* ---------------------------------------------------------------------- */ double neuromat_eeg_geom_3D_dist_sqr(r2_t *p, r2_t *q); /* A distance between points {p} and {q} on the sphere, squared. Specifically, the result is {2*(1 - cos(t))} where {t} is the angle between the two unit vectors {P} and {Q} of {R^3} corresponding to {p} and {q} by {neuromat_eeg_geom_3D_from_2D}. If {p} and {q} are close, the result is approximately the square of the Euclidean distance from {P} to {Q}. In particular, the result is 0 if {P==Q} (that is, {p==q}), and 4 (maximum) if {P} and {Q} are diametrally opposite. */ double neuromat_eeg_geom_3D_dist_sqr(r2_t *p, r2_t *q) { r3_t P = neuromat_eeg_geom_3D_from_2D(p); r3_t Q = neuromat_eeg_geom_3D_from_2D(q); double ct = r3_dot(&P, &Q); /* Co-sine of angle. */ return 2*fmax(1 - ct, 0); } /* ---------------------------------------------------------------------- */ void neuromat_filter_apply(int nt, int ne, double **val, double fsmp, neuromat_filter_t *gain, bool_t padding) { demand(fsmp > 0.0, "invalid sampling frequency {fsmp}"); /* Compute size {ns} of signal includeing mirror padding at ends: */ int nr; /* Num samples in end-mirrored segments. */ int ns; /* Total samples, including mirrored ends and zero-padding. */ if (padding) { int nr_max = 1024; /* Max signal replication at ends. */ nr = (int)ceil(fsmp); /* Signal replication at ends: one second's worth. */ if (nr > nr_max) { nr = nr_max; } if (nr > nt) { nr = nt; } ns = nt + 2*nr; /* Round {ns} up to a multiple of 1024: */ ns = ((ns + 1023)/1024)*1024; } else { nr = 0; ns = nt; } /* Allocate the FFT work areas and plans: */ double *in = (double*) fftw_malloc(sizeof(double)*(ns+4)); double *out = (double*) fftw_malloc(sizeof(double)*(ns+4)); fftw_plan pd = fftw_plan_r2r_1d(ns, in, out, FFTW_DHT, FFTW_ESTIMATE | FFTW_DESTROY_INPUT); fftw_plan pi = fftw_plan_r2r_1d(ns, out, in, FFTW_DHT, FFTW_ESTIMATE | FFTW_DESTROY_INPUT); /* Apodizing half-window for replicated ends: */ double *W = notnull(malloc(nr*sizeof(double)), "no mem"); int r; for (r = 0; r < nr; r++) { W[r] = 0.5*(1 + cos(M_PI*(r+0.5)/nr)); } /* Fourier transfer function: */ double *G = neuromat_filter_tabulate_gain(ns, fsmp, gain); int it, ie; for (ie = 0; ie < ne; ie++) { /* Copy the signal {ie} into the FFT buffer, with mirrored apodized ends and zero padding: */ for (it = 0; it < nt; it++) { double vt = val[it][ie]; if (it < nr) { in[nr - 1 - it] = vt*W[it]; } in[nr + it] = vt; if (it >= nt - nr) { in[nr + 2*nt - it] = vt*W[nt - 1 - it]; } } for (it = nt + 2*nr; it < ns; it++) { in[it] = 0.0; } /* Transform to frequency domain: */ fftw_execute(pd); /* Apply frequency filter: */ int kf0; for (kf0 = 0; kf0 <= ns-kf0; kf0++) { int kf1 = (ns - kf0) % ns; /* The other Hartley element with same absolute freq. */ double h0 = out[kf0]; double G0 = G[kf0]; if (kf1 == kf0) { /* Coef {kf0} is the only one with that frequency: */ out[kf0] = h0*G0; } else { /* Coefs {kf0,kf1} have the same frequency and get mixed: */ double h1 = out[kf1]; double G1 = G[kf1]; /* Apply filter: */ double Gp = G0 + G1; double Gm = G0 - G1; out[kf0] = 0.5*(h0*Gp + h1*Gm); out[kf1] = 0.5*(h1*Gp - h0*Gm); } } /* Return to time domain: */ fftw_execute(pi); /* Extract signal and return it to {val}: */ for (it = 0; it < nt; it++) { val[it][ie] = in[nr + it]/ns; } } fftw_destroy_plan(pd); fftw_destroy_plan(pi); fftw_free(out); fftw_free(in); if (W != NULL) { free(W); } free(G); return; } /* ---------------------------------------------------------------------- */ double neuromat_filter_lowpass_soft(double f, double fa, double fb); /* An improvised lowpass filter. If {fa==fb==0}, returns */ double neuromat_filter_lowpass_soft(double f, double fa, double fb) { demand((0 <= fa) && (fa <= fb), "critical frequencies out of order"); if (f >= fb) { return 0.0; } else if (f <= fa) { return 1.0; } else { double zw = fb - fa; double z = f - fa; assert(zw > 0); double r = z/zw; double s = 1-r*r; if (s < 0.025) { return 0.0; } else { return exp((s-1)/s); } } } /* ---------------------------------------------------------------------- */ int neuromat_eeg_frame_buffer_get { /* If {raw=FALSE}, assumes that {rd} is in plain EEG-20 text format, as per {neuromat_eeg_frame_read}, positioned after any header lines. If {raw=TRUE}, assumes {rd} is in the ".raw" (binary) EEG-128 format, as per {neuromat_eeg_raw_frame_read}. In that case, {version} specifies the sample format, either 16-bit signed integers ({version=2}) or 32-bit IEEE single-precision floats ({version=4}). */ int nl = (nlP == NULL ? 0 : (*nlP)); /* Number of file lines already read (incl. comments and headers). */ int nf = (nfP == NULL ? 0 : (*nfP)); /* Number of data frames already read. */ if (raw) { double unit = 1.0; nfr = neuromat_eeg_raw_frame_read(rd, version, unit, buf->nc, NULL, frm); nf++; } else { nfr = neuromat_eeg_frame_read(rd, buf->nc, frm, &nl, &nf); } /* Update line and frame counters: */ (*nlP) = nl; (*nfP) = nf; } /* ---------------------------------------------------------------------- */ void neuromat_locate_runs ( int nt, /* Number of data frames. */ int nc, /* Number of data channels. */ int ic_trig, /* Index of trigger channel. */ double **val, /* The EEG dataset ({nt} frames with {nc} channels each). */ int nr, /* Expected number of runs in file. */ int np, /* Number of trigger pulses per run. */ int it_pulse_ini[], /* (OUT) Index of first frame in each pulse. */ int it_pulse_fin[], /* (OUT) Index of last frame in each pulse. */ int it_run_ini[], /* (OUT) Index of first frame in each run. */ int it_run_fin[] /* (OUT) Index of last frame in each run. */ ) { demand((ic_trig) && (ic_trig < nc), "invalid trigger channel index"); demand(np >0, "each run must have at least one trigger pulse"); int ir = 0; /* Index of next run. */ int kp = 0; /* Number of trigger-up events seen since last run. */ double vtr_prev = 0; /* Trigger value in previous frame. */ int it; for (it = 0; it <= nt; it++) { double vtr = (it < nt ? val[it][ic_trig] : 0.0); demand(vtr >= 0, "negative trigger value"); if ((vtr_prev == 0) && (vtr > 0)) { /* Trigger up-event: */ demand(it > 0, "incomplete trigger pulse at start of file"); assert(kp < np); it_pulse_ini[ir*np + kp] = it; if (kp == 0) { /* Start of run: */ demand(ir < nr, "too many trigger events"); it_run_ini[ir] = it; } kp++; } else if ((vtr_prev > 0) && (vtr == 0)) { /* Trigger down-event: */ demand(it < nt, "incomplete trigger pulse at end of file"); assert(kp > 0); assert(kp <= np); assert(it > 0); it_pulse_fin[ir*np + (kp-1)] = it-1; if (kp == np) { /* End of run: */ assert(ir < nr); it_run_fin[ir] = it-1; kp = 0; ir++; } } vtr_prev = vtr; } demand(kp == 0, "incomplete run at end of file"); demand(ir == nr, "insufficient number of runs in file"); } /* ---------------------------------------------------------------------- */ void neuromat_eeg_func_basis_compute_shepard_weights(int32_t ne, double bval[], r3_t *p3D, r3_t pos3D[], double erad[]); /* Evaluates {ne} Shepard radial elements, for centers {pos3D[0..ne-1]} and scaling radii {erad[0..ne-1]}, at the point {*p}. Assumes that all those points are on the unit sphere. Returns the weights in {bval[0..ne-1]}, not normalized. */ void neuromat_eeg_func_basis_compute_local_affine_approxs(int32_t ne, double bval[], double wtel[], r3_t *p3D, r3_t pos3D[], double U[]); /* For each {k} in {0..ne-1}, adjusts by least squares an affine real function {F[k]} of {\RR^3} to the data samples {(pos3D[i], (i==k))} with weights {wtel[0..ne-1]}, for {i} in {0..ne-1}. Then sets {bval[k]} to {Fk](p3D)}. The matrix {U} should have {4*ne} elements and is used internally by the call. */ void neuromat_eeg_func_basis_eval_radial(int32_t ne, double bval[], r3_t *p3D, r3_t pos3D[], double erad[], neuromat_eeg_func_basis_mother_t mother) { /* Set {rval[0..ne-1]} to the raw radial basis values of each element at {p3D}: */ for (int32_t ie = 0; ie < ne; ie++) { bval[ie] = mother(ie, p3D); } } void neuromat_eeg_func_basis_compute_local_affine_approxs(int32_t ne, double bval[], double wtel[], r3_t *p3D, r3_t pos3D[], double U[]) { int32_t nb = 4; /* Number of elements in the LSQ fitting basis. */ auto double phi(int32_t j, r3_t *q); /* Computes a basis function {\phi[j]} at point {q} of {\RR^3}, for {j} in {0..nb-1}. Namely, {\phi[0](q) = 1}, {\phi[1+j](q) = q.c[j] - p3D.c[j]}. */ /* Let {\phi[j]} for {j} in {0..nb-1} be a basis for the affine real functions on {\RR^3}, specifically {\phi[0](p) = 1} and {\phi[1+j](p) = p.c[j]} for {j} in {0..2}. We compute for each {k} in {0..ne-1} a real affine function {s[k]} on {\RR^3} that best approximates (in the LSQ sense) the canonical data for electrode {k}, namely 1 at {pos3D[k]} and 0 at each {pos3D[i]} with {i!=k}, where the weights for these data points are {wtel[0..ne-1]}. The approximation, for each {k}, will be { s[k] = \sum_{j=0}^{3} U[j][k]*\phi[k] } */ double phi(int32_t j, r3_t *q) { assert((j >= 0) && (j < nb)); if (j == 0) { return 1.0; } else { double dj = q->c[j-1] - p3D->c[j-1]; return dj; } } auto void generate_case(int32_t i, int32_t nv, double vi[], int32_t nf, double fi[], double *wiP); /* Generates the LSQ basis values {vi[0..nv-1]} at electrode {i}, the corresponding function values {fi[0..nf-1]}, and the weight {*wiP}. Namely, sets {vi[j] = \phi[j](pos3D[i])}, {fi[k] = (k==i)}, and {*wiP = wtel[i]}. Requires {nv==nb} and {nf==ne}. */ /* Compute the linear function: */ bool_t verbose_lsq = FALSE; lsq_fit(nb, ne, ne, generate_case, U, verbose_lsq); /* Now evaluate the linear functions at {p}: */ /* Should be just {U[0,ie]} since {\phi(p) = (1,0,0,0)}, but let's play safe: */ double phip[nb]; for (int32_t j = 0; j < nb; j++) { phip[j] = phi(j, p3D); } rmxn_map_row(nb, ne, phip, U, bval); return; void generate_case(int32_t i, int32_t nv, double vi[], int32_t nf, double fi[], double *wiP) { assert(nv == nb); assert(nf == ne); assert((i >= 0) && (i < ne)); for (int32_t j = 0; j < nv; j++) { vi[j] = phi(j, &(pos3D[i])); } for (int32_t k = 0; k < nf; k++) { fi[k] = (k == i ? 1 : 0); } (*wiP) = wtel[i]; } } void neuromat_eeg_func_basis_eval_shepard_0(int32_t ne, double bval[], r3_t *p3D, r3_t pos3D[], double erad[]) { /* Set {bval[0..ne-1]} to the Shepard weights of each element at {p3D}: */ neuromat_eeg_func_basis_compute_shepard_weights(ne, bval, p3D, pos3D, erad); } void neuromat_eeg_func_basis_eval_shepard_1(int32_t ne, double bval[], r3_t *p3D, r3_t pos3D[], double erad[], double U[]) { /* Set {bval[0..ne-1]} to the Shepard weights of each element at {p3D}: */ neuromat_eeg_func_basis_compute_shepard_weights(ne, bval, p3D, pos3D, erad); /* Now compute the actual values by LSQ approximation: */ neuromat_eeg_func_basis_compute_local_affine_approxs(ne, bval, bval, p3D, pos3D, U); }