Hardware efficient automatic thresholding for NEO-based neural spike detection (C version)

neo.cpp

/*
 * C implementation of
 * Yang, Yuning, and Andrew J. Mason. "Hardware efficient automatic
 * thresholding for NEO-based neural spike detection." IEEE Transactions
 * on Biomedical Engineering 64.4 (2016): 826-833.
 *
 * Writted by Lennart P. L. Landsmeer
 * NeurocomputingLab, Erasmus MC
 * */
#include <stdint.h>
#include <stdio.h>

#define STDCAL_WINDOWSIZE 127
#define STDCAL_Z (1 << 5) // for division instead of right shift
#define STDCAL_GAIN (1 << 6) // for division instead of right shift

// float version of the neo filter for comparison
struct neo_filter_f { float y = 0; float yprev = 0; float yprevprev = 0; float sneo = 0; float sneoprev = 0; }; int neo_stream_f(struct neo_filter_f * state, float x, float thresh) { const float alpha_filter1 = 1/4.f; const float alpha_filter2 = 3/32.f; state->y = (1 - alpha_filter1) * state->y + alpha_filter1 * x; float neo = state->yprev * state->yprev - state->yprevprev * state->y; float sneoprev = state->sneo; state->sneo = (1 - alpha_filter2) * state->sneo  + alpha_filter2 * neo; state->yprevprev = state->yprev; state->yprev = state->y; return state->sneo >= thresh && sneoprev < thresh; }

struct rms_filter_s16 {
    // Section III.B
    int16_t y = 0;  // exponential filter of neural signal
    int16_t yprev = 0; // needed for zerocross detection
    uint32_t zerocross = 0; // no of zerocrossing in window (2^15)
    uint32_t sample_count = 0; // sample count needed know when we hit a window edge
    float output = 0; // actual output (a float, in the order of ~0.25)
};
void rms_stream_s16(rms_filter_s16 * state, int16_t x) {
    // x: raw neural signal
    float Nc = 9.5;
    state->y = (3 * state->y + 1 * x) / 4; // exp(alpha=1/4)
    state->sample_count += 1;
    if ((state->yprev < 0 && state->y >= 0) ||
        (state->yprev > 0 && state->y <= 0)) {
        state->zerocross += 1;
    }
    state->yprev = state->y;
    if (state->sample_count % (1 << 15) == 0) {
        // (100 z)^2/2^40*Nc
        state->output = (float) state->zerocross * (float) state->zerocross * 9.094947017729282e-09f * Nc;
        state->zerocross = 0;
    }
}

struct stdcal_filter_s16 {
    // Section III.C
    // missing: convergence detection
    int sample_count = 0; // sample counter (run filter each windowsize samples)
    int16_t y = 0; // exponential filter of neural signal
    int16_t window[STDCAL_WINDOWSIZE]; // we fill up a buffer (window)
    int16_t stdn = 0; // final result (just right shifted filterout)
    int16_t delta2 = 0; // previous delta value, shifted right (eq 19)
    int16_t filterout = 0; // the digital integrator (eq 19)
};
void stdcal_stream_s16(stdcal_filter_s16 * state, int16_t x) {
    // x: raw neural signal
    state->y = (3 * state->y + 1 * x) / 4; // exp(alpha=1/4)
    state->window[state->sample_count++ % STDCAL_WINDOWSIZE] = state->y;
    // apply filter 1/4 here
    if (state->sample_count % STDCAL_WINDOWSIZE != 0 || state->sample_count <= STDCAL_WINDOWSIZE) return;
    int delta = - STDCAL_WINDOWSIZE * 0.159;
    for (int i = 0; i < STDCAL_WINDOWSIZE; i++) {
        delta += state->window[i] >= state->stdn ? 1 : 0;
    }
    state->filterout = state->filterout + delta - state->delta2;
    state->stdn = state->filterout / STDCAL_GAIN;
    state->delta2 = delta / STDCAL_GAIN;
}

struct neo_filter_s16 {
    // Section II (2)
    int16_t y = 0; // first filter, alpha = 1/4
    int16_t yprev = 0; // time lag one of first filter for neo
    int16_t yprevprev = 0; // time lag two of first filter for neo
    int16_t sneo = 0; // second filter on neo, alpha = 3/32
    int16_t sneoprev = 0; // time lag one on second filter for hit detection
};
int neo_stream_s16(struct neo_filter_s16 * state, int16_t x, int16_t thresh) {
    // x: raw neural signal
    // thresh: calculated threshold value
    // returns: 1 if spike detected, else 0
    state->y = (3 * state->y + 1 * x) / 4; // exp(alpha=1/4)
    int16_t neo = state->yprev * state->yprev - state->yprevprev * state->y;
    int16_t sneoprev = state->sneo;
    state->sneo = (29 * state->sneo  + 3 * neo) / 32; // exp(alpha=3/32)
    state->yprevprev = state->yprev;
    state->yprev = state->y;
    return state->sneo >= thresh && sneoprev < thresh;
}

int main() {
    // setup i/o
    unsigned int x;
    int16_t x16;
    FILE * f = fopen("data.txt", "r");
    FILE * o = fopen("cout.txt", "w");
    FILE * c = fopen("cstd.txt", "w");

    // obtain threshold (execute every 5 mins)
    rms_filter_s16 rms;
    stdcal_filter_s16 stdcal;
    while (1 == fscanf(f, "%x", &x)) {
        x16 = x > 511 ? x - 1024 : x;
        stdcal_stream_s16(&stdcal, x16);
        rms_stream_s16(&rms, x16);
        fprintf(c, "%d %f\n", stdcal.stdn, rms.output);
    }
    int16_t thresh = stdcal.stdn * stdcal.stdn * rms.output;

    rewind(f);
    // actual spike detection
    neo_filter_f neof;
    neo_filter_s16 neo;
    while (1 == fscanf(f, "%x", &x)) {
        x16 = x > 511 ? x - 1024 : x;
        int hit  = neo_stream_s16(&neo, x16, thresh);
        int hitf  = neo_stream_f(&neof, x16, thresh);
        fprintf(o, "%d\t%d\t|\t%f\t%d\n", neo.sneo, hit, neof.sneo, hitf);
    }

    // cleanup
    fclose(f);
    fclose(o);
    fclose(c);
}