Created
March 25, 2023 12:24
-
-
Save llandsmeer/355bed4f47d332ce1dd096aa232b60c0 to your computer and use it in GitHub Desktop.
Hardware efficient automatic thresholding for NEO-based neural spike detection (C version)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
* 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); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment