Add channel rereferencing across language bindings

cpp_package/examples/signal_processing/src/downsampling.cpp

@@ -45,6 +45,14 @@ int main (int argc, char *argv[])
         int filtered_size = 0;
         std::vector<int> eeg_channels = BoardShim::get_eeg_channels (board_id);
 
+        // Re-reference every EEG channel using the sample-wise mean of the first two EEG channels.
+        // The operation changes data in-place and may be applied before other signal processing.
+        if (eeg_channels.size () >= 2)
+        {
+            std::vector<int> reference_channels = {eeg_channels[0], eeg_channels[1]};
+            DataFilter::reference (data, eeg_channels, reference_channels);
+        }
+
         for (int i = 0; i < eeg_channels.size (); i++)
         {
             std::cout << "Data from :" << eeg_channels[i] << " before downsampling " << std::endl;

cpp_package/src/data_filter.cpp

@@ -7,6 +7,65 @@
 #include "data_handler.h"
 
 
+void DataFilter::reference (BrainFlowArray<double, 2> &data,
+    const std::vector<int> &channels_to_reference, const std::vector<int> &reference_channels)
+{
+    const int rows = data.get_size (0);
+    const int cols = data.get_size (1);
+    if ((rows <= 0) || (cols <= 0) || channels_to_reference.empty () ||
+        reference_channels.empty ())
+    {
+        throw BrainFlowException (
+            "data and channel lists must be non-empty",
+            (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+    }
+
+    for (const int channel : channels_to_reference)
+    {
+        if ((channel < 0) || (channel >= rows))
+        {
+            throw BrainFlowException (
+                "channel index is out of range",
+                (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+        }
+    }
+    for (const int channel : reference_channels)
+    {
+        if ((channel < 0) || (channel >= rows))
+        {
+            throw BrainFlowException (
+                "reference channel index is out of range",
+                (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+        }
+    }
+
+    // Compute the complete reference signal before mutating any channel. This keeps the result
+    // correct when a reference channel is also present in channels_to_reference.
+    std::vector<double> reference_signal (cols, 0.0);
+    for (const int channel : reference_channels)
+    {
+        const double *channel_data = data.get_address (channel);
+        for (int sample = 0; sample < cols; sample++)
+        {
+            reference_signal[sample] += channel_data[sample];
+        }
+    }
+    const double scale = 1.0 / static_cast<double> (reference_channels.size ());
+    for (double &sample : reference_signal)
+    {
+        sample *= scale;
+    }
+
+    for (const int channel : channels_to_reference)
+    {
+        double *channel_data = data.get_address (channel);
+        for (int sample = 0; sample < cols; sample++)
+        {
+            channel_data[sample] -= reference_signal[sample];
+        }
+    }
+}
+
 double DataFilter::get_oxygen_level (double *ppg_ir, double *ppg_red, int data_len,
     int sampling_rate, double coef1, double coef2, double coef3)
 {

cpp_package/src/inc/data_filter.h

@@ -28,6 +28,16 @@ class DataFilter
     /// write user defined string to BrainFlow logger
     static void log_message (int log_level, const char *format, ...);
 
+    /**
+     * re-reference selected channels in-place
+     * @param data input 2d array, rows are channels and columns are samples
+     * @param channels_to_reference channel rows from which to subtract the reference signal
+     * @param reference_channels channel rows whose sample-wise mean defines the reference signal
+     */
+    static void reference (BrainFlowArray<double, 2> &data,
+        const std::vector<int> &channels_to_reference,
+        const std::vector<int> &reference_channels);
+
     /// perform low pass filter in-place
     static void perform_lowpass (double *data, int data_len, int sampling_rate, double cutoff,
         int order, int filter_type, double ripple);

csharp_package/brainflow/brainflow/data_filter.cs

@@ -76,6 +76,61 @@ public static void set_log_file (string log_file)
                 throw new BrainFlowError (res);
             }
         }
+
+        /// <summary>
+        /// Re-reference selected channels in-place using the sample-wise mean of reference channels.
+        /// </summary>
+        /// <param name="data">Rows are channels and columns are samples.</param>
+        /// <param name="channels_to_reference">Rows from which to subtract the reference signal.</param>
+        /// <param name="reference_channels">Rows whose mean defines the reference signal.</param>
+        public static void reference (double[,] data, int[] channels_to_reference, int[] reference_channels)
+        {
+            if ((data == null) || (channels_to_reference == null) || (reference_channels == null) ||
+                (data.GetLength (0) == 0) || (data.GetLength (1) == 0) ||
+                (channels_to_reference.Length == 0) || (reference_channels.Length == 0))
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+
+            int rows = data.GetLength (0);
+            int cols = data.GetLength (1);
+            foreach (int channel in channels_to_reference)
+            {
+                if ((channel < 0) || (channel >= rows))
+                {
+                    throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+                }
+            }
+            foreach (int channel in reference_channels)
+            {
+                if ((channel < 0) || (channel >= rows))
+                {
+                    throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+                }
+            }
+
+            // Snapshot the reference before changing any row so overlapping channel lists are safe.
+            double[] reference_signal = new double[cols];
+            foreach (int channel in reference_channels)
+            {
+                for (int sample = 0; sample < cols; sample++)
+                {
+                    reference_signal[sample] += data[channel, sample];
+                }
+            }
+            for (int sample = 0; sample < cols; sample++)
+            {
+                reference_signal[sample] /= reference_channels.Length;
+            }
+
+            foreach (int channel in channels_to_reference)
+            {
+                for (int sample = 0; sample < cols; sample++)
+                {
+                    data[channel, sample] -= reference_signal[sample];
+                }
+            }
+        }
         // accord GetRow returns a copy instead pointer, so we can not easily update data in place like in other bindings
 
         /// <summary>

csharp_package/brainflow/examples/downsampling/downsampling.cs

@@ -22,6 +22,15 @@ static void Main (string[] args)
             board_shim.stop_stream ();
             double[,] unprocessed_data = board_shim.get_board_data ();
             int[] eeg_channels = BoardShim.get_eeg_channels (board_id);
+
+            // Re-reference every EEG channel using the sample-wise mean of the first two EEG channels.
+            // The operation changes unprocessed_data in-place.
+            if (eeg_channels.Length >= 2)
+            {
+                int[] reference_channels = {eeg_channels[0], eeg_channels[1]};
+                DataFilter.reference (unprocessed_data, eeg_channels, reference_channels);
+            }
+
             board_shim.release_session ();
 
             for (int i = 0; i < eeg_channels.Length; i++)

docs/SwiftAPIParity.rst

@@ -29,8 +29,8 @@ Implemented:
 
 - Filters, noise removal, and detrending: lowpass, highpass, bandpass, bandstop, environmental
   noise removal, rolling filter, and detrend.
-- Transforms and features: downsampling, wavelet transform/inverse/denoising, CSP, windowing,
-  FFT/IFFT, PSD/Welch, band powers, and ICA.
+- Transforms and features: channel re-referencing, downsampling, wavelet
+  transform/inverse/denoising, CSP, windowing, FFT/IFFT, PSD/Welch, band powers, and ICA.
 - Helpers: standard deviation, railed percentage, oxygen level, heart rate, peak detection,
   nearest power of two, file I/O, reshape helpers, logging, and version.
 

java_package/brainflow/src/main/java/brainflow/DataFilter.java

@@ -192,6 +192,74 @@ public static void set_log_file (String log_file) throws BrainFlowError
         }
     }
 
+    /**
+     * Re-reference selected channels in-place using the sample-wise mean of reference channels.
+     *
+     * @param data rows are channels and columns are samples
+     * @param channels_to_reference rows from which to subtract the reference signal
+     * @param reference_channels rows whose mean defines the reference signal
+     */
+    public static void reference (double[][] data, int[] channels_to_reference, int[] reference_channels)
+            throws BrainFlowError
+    {
+        if ((data == null) || (channels_to_reference == null) || (reference_channels == null)
+                || (data.length == 0) || (data[0] == null) || (data[0].length == 0)
+                || (channels_to_reference.length == 0) || (reference_channels.length == 0))
+        {
+            throw new BrainFlowError ("Invalid reference arguments",
+                    BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+
+        final int rows = data.length;
+        final int cols = data[0].length;
+        for (double[] row : data)
+        {
+            if ((row == null) || (row.length != cols))
+            {
+                throw new BrainFlowError ("Data must be a rectangular matrix",
+                        BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+            }
+        }
+        for (int channel : channels_to_reference)
+        {
+            if ((channel < 0) || (channel >= rows))
+            {
+                throw new BrainFlowError ("Channel index is out of range",
+                        BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+            }
+        }
+        for (int channel : reference_channels)
+        {
+            if ((channel < 0) || (channel >= rows))
+            {
+                throw new BrainFlowError ("Reference channel index is out of range",
+                        BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+            }
+        }
+
+        // Snapshot the reference before changing any row so overlapping channel lists are safe.
+        double[] reference_signal = new double[cols];
+        for (int channel : reference_channels)
+        {
+            for (int sample = 0; sample < cols; sample++)
+            {
+                reference_signal[sample] += data[channel][sample];
+            }
+        }
+        for (int sample = 0; sample < cols; sample++)
+        {
+            reference_signal[sample] /= reference_channels.length;
+        }
+
+        for (int channel : channels_to_reference)
+        {
+            for (int sample = 0; sample < cols; sample++)
+            {
+                data[channel][sample] -= reference_signal[sample];
+            }
+        }
+    }
+
     /**
      * calc stddev
      */

java_package/brainflow/src/main/java/brainflow/examples/Downsampling.java

@@ -30,6 +30,14 @@ public static void main (String[] args) throws Exception
         board_shim.release_session ();
 
         int[] eeg_channels = BoardShim.get_eeg_channels (board_id);
+        // Re-reference every EEG channel using the sample-wise mean of the first two EEG channels.
+        // The operation changes data in-place.
+        if (eeg_channels.length >= 2)
+        {
+            int[] reference_channels = {eeg_channels[0], eeg_channels[1]};
+            DataFilter.reference (data, eeg_channels, reference_channels);
+        }
+
         for (int i = 0; i < eeg_channels.length; i++)
         {
             System.out.println ("Original data:");

julia_package/brainflow/src/data_filter.jl

@@ -141,6 +141,43 @@ end
 WaveletType = Union{WaveletTypes, Integer}
 
 
+"""
+    reference(data, channels_to_reference, reference_channels)
+
+Re-reference selected channel rows in-place using the sample-wise mean of
+`reference_channels`. Channel indices follow Julia's one-based convention.
+"""
+function reference(data::AbstractMatrix{<:AbstractFloat}, channels_to_reference,
+                   reference_channels)
+    rows, cols = size(data)
+    if rows == 0 || cols == 0 || isempty(channels_to_reference) || isempty(reference_channels)
+        throw(BrainFlowError("Data and channel lists must be non-empty", Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    if !all(channel -> channel isa Integer && 1 <= channel <= rows, channels_to_reference)
+        throw(BrainFlowError("Channel index is out of range", Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    if !all(channel -> channel isa Integer && 1 <= channel <= rows, reference_channels)
+        throw(BrainFlowError("Reference channel index is out of range", Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+
+    # Snapshot the reference before changing any row so overlapping channel lists are safe.
+    reference_signal = zeros(eltype(data), cols)
+    for channel in reference_channels
+        for sample in 1:cols
+            reference_signal[sample] += data[channel, sample]
+        end
+    end
+    reference_signal ./= length(reference_channels)
+
+    for channel in channels_to_reference
+        for sample in 1:cols
+            data[channel, sample] -= reference_signal[sample]
+        end
+    end
+    return nothing
+end
+
+
 @brainflow_rethrow function perform_lowpass(data, sampling_rate::Integer, cutoff::Float64, order::Integer,
     filter_type::FiltType, ripple::Float64)
     ccall((:perform_lowpass, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),

julia_package/brainflow/test/downsampling.jl

@@ -16,9 +16,15 @@ BrainFlow.release_session(board_shim)
 eeg_channels = BrainFlow.get_eeg_channels(BrainFlow.SYNTHETIC_BOARD)
 sampling_rate = BrainFlow.get_sampling_rate(BrainFlow.SYNTHETIC_BOARD)
 
+# Re-reference every EEG channel using the sample-wise mean of the first two EEG channels.
+# The operation changes data in-place.
+if length(eeg_channels) >= 2
+    BrainFlow.reference(data, eeg_channels, eeg_channels[1:2])
+end
+
 data_first_channel = data[eeg_channels[1], :]
 println("Original Data First Channel")
 println(data_first_channel)
 downsampled_data = BrainFlow.perform_downsampling(data_first_channel, 3, BrainFlow.EACH)
 println("After Downsampling")
-println(downsampled_data)
+println(downsampled_data)

julia_package/brainflow/test/julia_tests.jl

@@ -43,6 +43,17 @@ end
     @test Int32(BrainFlow.BLACKMAN_HARRIS) == 3
 end
 
+@testset "channel reference" begin
+    data = [1.0 2.0; 3.0 4.0; 5.0 6.0]
+    BrainFlow.reference(data, [1, 3], [1, 2])
+    @test data == [-1.0 -1.0; 3.0 4.0; 3.0 3.0]
+
+    @test_throws BrainFlow.BrainFlowError BrainFlow.reference(copy(data), Int[], [1])
+    @test_throws BrainFlow.BrainFlowError BrainFlow.reference(copy(data), [1], Int[])
+    @test_throws BrainFlow.BrainFlowError BrainFlow.reference(copy(data), [4], [1])
+    @test_throws BrainFlow.BrainFlowError BrainFlow.reference(copy(data), [1], [4])
+end
+
 @testset "presets" begin
     presets = BrainFlow.get_board_presets(BrainFlow.CYTON_BOARD)
     @test presets[1] == Int32(BrainFlow.DEFAULT_PRESET)
@@ -52,4 +63,4 @@ end
     params = BrainFlowModelParams(BrainFlow.RESTFULNESS, BrainFlow.DEFAULT_CLASSIFIER)
     @test params.metric == BrainFlow.RESTFULNESS
     @test params.classifier == BrainFlow.DEFAULT_CLASSIFIER
-end
+end