From d27a5be819e69469fb40b332335b1c35f8854b62 Mon Sep 17 00:00:00 2001
From: shayo <shay.ohayon@gmail.com>
Date: Wed, 18 Feb 2015 10:24:00 -0500
Subject: [PATCH] Added Standard Deviation Display to Spike Sorter
Useful to estimate the noise levels on channels and to make a somewhat
objective decision on threshold.
Computation is very fast (running statistics). Right click clears up the
stats. Display is on the top left corner.
Currently only works for single channel electrodes (i.e., always
displays the 0th channel statistics).
---
Source/Processors/SpikeSorter/SpikeSorter.cpp | 32 ++++++++--
Source/Processors/SpikeSorter/SpikeSorter.h | 61 ++++++++++++++++++-
.../SpikeSorter/SpikeSorterCanvas.cpp | 17 +++---
3 files changed, 95 insertions(+), 15 deletions(-)
diff --git a/Source/Processors/SpikeSorter/SpikeSorter.cpp b/Source/Processors/SpikeSorter/SpikeSorter.cpp
index 813bdc28c..3fab6eb2a 100644
--- a/Source/Processors/SpikeSorter/SpikeSorter.cpp
+++ b/Source/Processors/SpikeSorter/SpikeSorter.cpp
@@ -221,6 +221,8 @@ Electrode::~Electrode()
delete voltageScale;
delete channels;
delete spikeSort;
+ delete runningStats;
+
}
Electrode::Electrode(int ID, UniqueIDgenerator *uniqueIDgenerator_, PCAcomputingThread *pth, String _name, int _numChannels, int *_channels, float default_threshold, int pre, int post, float samplingRate , int sourceNodeId)
@@ -238,6 +240,7 @@ Electrode::Electrode(int ID, UniqueIDgenerator *uniqueIDgenerator_, PCAcomputing
isActive = new bool[numChannels];
channels = new int[numChannels];
voltageScale = new double[numChannels];
+ runningStats = new RunningStat[numChannels];
depthOffsetMM = 0.0;
advancerID = -1;
@@ -860,6 +863,25 @@ void SpikeSorter::handleEvent(int eventType, MidiMessage& event, int sampleNum)
// }
// }
+
+float SpikeSorter::getSelectedElectrodeNoise()
+{
+ if (electrodes.size() == 0)
+ return 0.0;
+
+ // TODO, change "0" to active channel to support tetrodes.
+ return electrodes[currentElectrode]->runningStats[0].StandardDeviation();
+}
+
+
+void SpikeSorter::clearRunningStatForSelectedElectrode()
+{
+ if (electrodes.size() == 0)
+ return;
+ // TODO, change "0" to active channel to support tetrodes.
+ electrodes[currentElectrode]->runningStats[0].Clear();
+}
+
void SpikeSorter::process(AudioSampleBuffer& buffer,
MidiBuffer& events)
{
@@ -893,23 +915,25 @@ void SpikeSorter::process(AudioSampleBuffer& buffer,
{
sampleIndex++;
-
+
// cycle through channels
for (int chan = 0; chan < electrode->numChannels; chan++)
{
// std::cout << " channel " << chan << std::endl;
-
+
if (*(electrode->isActive+chan))
{
//float v = getNextSample(currentChannel);
int currentChannel = electrode->channels[chan];
+ float currentValue = getNextSample(currentChannel);
+ electrode->runningStats[chan].Push(currentValue);
bool bSpikeDetectedPositive = electrode->thresholds[chan] > 0 &&
- (getNextSample(currentChannel) > electrode->thresholds[chan]); // rising edge
+ (currentValue > electrode->thresholds[chan]); // rising edge
bool bSpikeDetectedNegative = electrode->thresholds[chan] < 0 &&
- (getNextSample(currentChannel) < electrode->thresholds[chan]); // falling edge
+ (currentValue < electrode->thresholds[chan]); // falling edge
if (bSpikeDetectedPositive || bSpikeDetectedNegative)
{
diff --git a/Source/Processors/SpikeSorter/SpikeSorter.h b/Source/Processors/SpikeSorter/SpikeSorter.h
index bb520b868..8650b05d3 100644
--- a/Source/Processors/SpikeSorter/SpikeSorter.h
+++ b/Source/Processors/SpikeSorter/SpikeSorter.h
@@ -81,6 +81,63 @@ private:
int globalUniqueID;
};
+/* snatched from http://www.johndcook.com/blog/standard_deviation/ */
+class RunningStat
+{
+public:
+ RunningStat() : m_n(0) {}
+
+ void Clear()
+ {
+ m_n = 0;
+ }
+
+ void Push(double x)
+ {
+ m_n++;
+
+ // See Knuth TAOCP vol 2, 3rd edition, page 232
+ if (m_n == 1)
+ {
+ m_oldM = m_newM = x;
+ m_oldS = 0.0;
+ }
+ else
+ {
+ m_newM = m_oldM + (x - m_oldM) / m_n;
+ m_newS = m_oldS + (x - m_oldM)*(x - m_newM);
+
+ // set up for next iteration
+ m_oldM = m_newM;
+ m_oldS = m_newS;
+ }
+ }
+
+ int NumDataValues() const
+ {
+ return m_n;
+ }
+
+ double Mean() const
+ {
+ return (m_n > 0) ? m_newM : 0.0;
+ }
+
+ double Variance() const
+ {
+ return ((m_n > 1) ? m_newS / (m_n - 1) : 0.0);
+ }
+
+ double StandardDeviation() const
+ {
+ return sqrt(Variance());
+ }
+
+private:
+ int m_n;
+ double m_oldM, m_newM, m_oldS, m_newS;
+};
+
class Electrode
{
public:
@@ -106,6 +163,7 @@ class Electrode
double *voltageScale;
//float PCArange[4];
+ RunningStat *runningStats;
SpikeHistogramPlot* spikePlot;
SpikeSortBoxes* spikeSort;
PCAcomputingThread *computingThread;
@@ -181,8 +239,9 @@ public:
/** Creates the SpikeSorterEditor. */
AudioProcessorEditor* createEditor();
+ float getSelectedElectrodeNoise();
+ void clearRunningStatForSelectedElectrode();
-
//void addNetworkEventToQueue(StringTS S);
void postEventsInQueue(MidiBuffer& events);
diff --git a/Source/Processors/SpikeSorter/SpikeSorterCanvas.cpp b/Source/Processors/SpikeSorter/SpikeSorterCanvas.cpp
index abfc13528..4379d7ab9 100644
--- a/Source/Processors/SpikeSorter/SpikeSorterCanvas.cpp
+++ b/Source/Processors/SpikeSorter/SpikeSorterCanvas.cpp
@@ -1203,7 +1203,7 @@ bool WaveformAxes::checkThreshold(const SpikeObject& s)
void WaveformAxes::clear()
{
-
+ processor->clearRunningStatForSelectedElectrode();
spikeBuffer.clear();
spikeIndex = 0;
int numSamples=40;
@@ -1668,16 +1668,13 @@ void WaveformAxes::paint(Graphics& g)
if (drawGrid)
drawWaveformGrid(g);
+
+ double noise = processor->getSelectedElectrodeNoise();
+ String d = "STD: " + String(noise, 2) + "uV";
+ g.setFont(Font("Small Text", 13, Font::plain));
+ g.setColour(Colours::white);
- if (channel == 0)
- {
- //double depth = processor->getSelectedElectrodeDepth();
- //String d = "Depth: "+String(depth,4) +" mm";
- //g.setFont(Font("Small Text", 13, Font::plain));
- // g.setColour(Colours::white);
-
- //g.drawText(d,10,10,150,20,Justification::left,false);
- }
+ g.drawText(d, 10, 10, 150, 20, Justification::left, false);
// draw the grid lines for the waveforms
// draw the threshold line and labels
--
GitLab