diff --git a/processors/spikefeatures/CMakeLists.txt b/processors/spikefeatures/CMakeLists.txt new file mode 100644 index 0000000..30f8902 --- /dev/null +++ b/processors/spikefeatures/CMakeLists.txt @@ -0,0 +1,10 @@ + +ADD_LIBRARY(spikefeatures "spikefeatures.cpp" ) + +find_package(OpenMP) +if(OpenMP_CXX_FOUND) + TARGET_LINK_LIBRARIES(spikefeatures PUBLIC OpenMP::OpenMP_CXX) + +else() + TARGET_LINK_LIBRARIES(spikefeatures) +endif() diff --git a/processors/spikefeatures/doc.yaml b/processors/spikefeatures/doc.yaml new file mode 100644 index 0000000..a1b3eca --- /dev/null +++ b/processors/spikefeatures/doc.yaml @@ -0,0 +1,48 @@ +Description: Detect spikes in a TimeSeries data stream and emits a column data type packet containing all selected features for the decoding. + +Long description: possible features available are timestamp, amplitude, slope, channel index, channel depth. These features are disabled by default and need to be enabled in the graph options. + +Input ports: + - name: data + type: TimeSerieType + slots: 1-385 + description: + +Output port: + - name: data + type: ColumnData + slots: 1-385 + description: data packet containing spikes with the selected features. The packet is expandable - aka nsamples is resized for the real number of spikes. + +Options: + - &threshold + name: threshold + type: double + default: 70.0 + description: Spike detection threshold in data units. + - name: invert signal + type: bool + default: True + description: Invert signal before spike detection. + - &peaklifetime + name: peak lifetime + type: unsigned int + default: 1 samples + description: Peak life time in samples + - name: features + type : vector of string + default: empty, required option + description: list of features to use in output + - name: channelmap + type: map of string / double + default: empty, mandatory when the depth feature is activated + description: mapping between channel number and channel depth + +States: + Static: + - <<: *threshold + shared: true + external access: write + - <<: *peaklifetime + shared: true + external access: write diff --git a/processors/spikefeatures/spikefeatures.cpp b/processors/spikefeatures/spikefeatures.cpp new file mode 100644 index 0000000..35b6ce5 --- /dev/null +++ b/processors/spikefeatures/spikefeatures.cpp @@ -0,0 +1,216 @@ +// --------------------------------------------------------------------- +// This file is part of falcon-core. +// +// Copyright (C) 2021-now Neuro-Electronics Research Flanders +// +// Falcon-server is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Falcon-server is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with falcon-core. If not, see . +// --------------------------------------------------------------------- + + +#include "spikefeatures.hpp" + +SpikeFeatures::SpikeFeatures() : IProcessor(){ + + add_option("features", features_, "Selection of features to compute.", true); + + add_option("channeldepths", channel_pos_, "Relation between channel number and depth."); + + add_option(THRESHOLD, initial_threshold_, + "Spike detection threshold in data units."); + add_option(PEAK_LIFETIME, initial_peak_lifetime_, + "Peak life time in samples"); + + add_option("invert_signal", invert_signal_, + "invert a signal to detect negative spikes"); + +} + +void SpikeFeatures::Configure(const GlobalContext &context){ + features_labels_ = {}; + default_features_ = YAML::Load( "{ " + "time: false," + "amplitude: false," + "slope: false," + "channel index: false," + "depth: false" + "}"); + + + for(auto label: features_()){ + + if (!default_features_[label]){ + throw ProcessingConfigureError("The " + label + " is not implemented.", name()); + } + + default_features_[label] = true; + features_labels_.push_back(label); + } + + if(default_features_["depth"].as() and channel_pos_().size() == 0 ){ + throw ProcessingConfigureError("If the depth feature is selected, the channelmap option " + "needs to give the mapping between channel label and depth.", name()); + + } +} + +void SpikeFeatures::CreatePorts(){ + data_in_port_ = create_input_port>( + TimeSeriesType::Capabilities(ChannelRange(1, MAX_NCHANNELS)), + PortInPolicy(SlotRange(1, MAX_NCHANNELS))); + + + decoding_out_port_ = create_output_port>( + ColumnsType::Parameters(features_labels_, 0, true), + PortOutPolicy(SlotRange(1, MAX_NCHANNELS))); + + threshold_ = create_static_state(THRESHOLD, initial_threshold_(), true, + Permission::WRITE); + + peak_lifetime_ = create_static_state(PEAK_LIFETIME, initial_peak_lifetime_(), + true, Permission::WRITE); +} + +void SpikeFeatures::CompleteStreamInfo() { + // check if we have the same number of input and output slots + + if (data_in_port_->number_of_slots() != decoding_out_port_->number_of_slots()) { + auto err_msg = "Number of output slots (" + + std::to_string(decoding_out_port_->number_of_slots()) + + ") on port '" + decoding_out_port_->name() + + "' does not match number of input slots (" + + std::to_string(data_in_port_->number_of_slots()) + + ") on port '" + data_in_port_->name() + "'."; + throw ProcessingStreamInfoError(err_msg, name()); + } + + for (slot_ = 0; slot_ < data_in_port_->number_of_slots(); ++slot_) { + decoding_out_port_->streaminfo(slot_).set_stream_parameters(data_in_port_->streaminfo(slot_)); + decoding_out_port_->streaminfo(slot_).set_parameters(ColumnsType::Parameters(features_labels_, 0, true)); + + if(default_features_["depth"].as()){ + for(auto chan: data_in_port_->prototype(slot_).labels()){ + try { + channel_pos_().at(chan); + } catch (const std::out_of_range& oor) { + throw ProcessingConfigureError("The depth feature is selected but the channelmap given in input " + "does not contain the corresponding depth to the channel " + chan, name()); + } + } + } + } +} + +void SpikeFeatures::Prepare(GlobalContext &context) { + // Create an independent spike detector by data stream input + spike_detectors_.clear(); + + for (slot_ = 0; slot_ < data_in_port_->number_of_slots(); ++slot_) { + + if(invert_signal_()){ + spike_detectors_.push_back( + std::make_unique(data_in_port_->prototype(slot_).ncolumns(), + initial_threshold_(), + initial_peak_lifetime_(), + dsp::algorithms::SpikeDetectionSign::DOWN)); + }else{ + spike_detectors_.push_back( + std::make_unique(data_in_port_->prototype(slot_).ncolumns(), + initial_threshold_(), + initial_peak_lifetime_(), + dsp::algorithms::SpikeDetectionSign::UP)); + } + + + } +} + + +void SpikeFeatures::Process(ProcessingContext &context) { +#pragma omp parallel +{ + TimeSeriesType::Data *data_in = nullptr; + ColumnsType::Data *data_out = nullptr; + std::string channel_label; + std::vector amp; + double max_index; + std::string feature_str; + size_t sample = 0; + bool should_break=false; + while (!context.terminated() and !should_break) { + #pragma omp for nowait + for (slot_ = 0; slot_ < data_in_port_->number_of_slots(); ++slot_) { + + if (!data_in_port_->slot(slot_)->RetrieveData(data_in) ){ + should_break=true; + continue; + } + // claim output data buckets + data_out = decoding_out_port_->slot(slot_)->ClaimData(true); + spike_detectors_[slot_]->set_threshold(threshold_->get()); + spike_detectors_[slot_]->set_peak_life_time(peak_lifetime_->get()); + + unsigned int spike_number = 0; + for (sample = 0; sample < data_in->nsamples();sample++) { + if (spike_detectors_[slot_]->is_spike::Data::sample_iterator>( + data_in->sample_timestamp(sample), + data_in->begin_sample(sample))) { + + amp = spike_detectors_[slot_]->amplitudes_detected_spike(); + max_index = std::max_element(amp.begin(), amp.end()) - amp.begin(); + channel_label = data_in->labels()[max_index]; + feature_str = "("; + + if(default_features_["time"].as()){ + data_out->set_data_sample(spike_number, "time", sample+data_in->hardware_timestamp()); + feature_str += "time: " + std::to_string(sample+data_in->hardware_timestamp()); + } + + if(default_features_["channel index"].as()){ + + data_out->set_data_sample(spike_number, "channel index", std::stod(channel_label)); + feature_str += "channel index: " + channel_label; + } + + if(default_features_["amplitude"].as()){ + data_out->set_data_sample(spike_number, "amplitude", amp[max_index]); + feature_str += "amplitude: " + std::to_string(data_out->data_sample(spike_number, "amplitude")); + } + + if(default_features_["slope"].as()){ + data_out->set_data_sample(spike_number, "slope", spike_detectors_[slot_]->slopes_detected_spike()[max_index]); + feature_str += "slope: " + std::to_string(data_out->data_sample(spike_number, "slope")); + } + + if(default_features_["depth"].as()){ + channel_label = data_in->labels()[max_index]; + data_out->set_data_sample(spike_number, "depth", channel_pos_().at(channel_label)); + feature_str += "depth: " + std::to_string(data_out->data_sample(spike_number, "depth")); + } + + LOG(DEBUG) << name() << " Spike detected : " + feature_str + ") ts= " << sample << " channel =" << channel_label; + spike_number++; + } + + + } + data_out->CloneTimestamps(*data_in); + decoding_out_port_->slot(slot_)->PublishData(); + data_in_port_->slot(slot_)->ReleaseData(); + } + } + } +} + + +REGISTERPROCESSOR(SpikeFeatures) diff --git a/processors/spikefeatures/spikefeatures.hpp b/processors/spikefeatures/spikefeatures.hpp new file mode 100644 index 0000000..280e556 --- /dev/null +++ b/processors/spikefeatures/spikefeatures.hpp @@ -0,0 +1,66 @@ +// --------------------------------------------------------------------- +// This file is part of falcon-core. +// +// Copyright (C) 2021-now Neuro-Electronics Research Flanders +// +// Falcon-server is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Falcon-server is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with falcon-core. If not, see . +// --------------------------------------------------------------------- + +#pragma once +#include "dsp/algorithms.hpp" +#include "iprocessor.hpp" +#include "columnsdata/columnsdata.hpp" +#include "timeseriesdata/timeseriesdata.hpp" + +class SpikeFeatures : public IProcessor { + // CONSTRUCTOR and OVERLOADED METHODS + public: + SpikeFeatures(); + void Configure(const GlobalContext &context) override; + void CreatePorts() override; + void CompleteStreamInfo() override; + void Prepare(GlobalContext &context) override; + void Process(ProcessingContext &context) override; + + // DATA PORTS + protected: + PortIn> *data_in_port_; + PortOut> *decoding_out_port_; + +protected: + std::vector> spike_detectors_; + + + // STATES +protected: + StaticState *threshold_; + StaticState *peak_lifetime_; + + // CONSTANTS +public: + const uint32_t MAX_NCHANNELS=384; + const uint32_t MAX_NSPIKES=100; + const std::string PEAK_LIFETIME = "peak lifetime"; + const std::string THRESHOLD = "threshold"; + std::vector features_labels_; + YAML::Node default_features_; + int slot_ = 0; + +private: + options::Double initial_threshold_{70.}; + options::Measurement initial_peak_lifetime_{1, "sample"}; + options::Value, false> features_{}; + options::Value, false> channel_pos_{}; + options::Bool invert_signal_{true}; +};