Add Python DataFilter.rereference_channels for EEG/EXG

.github/workflows/run_unix.yml

@@ -390,6 +390,7 @@ jobs:
         sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/decomposition_regression.py
         sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/signal_metrics_regression.py
         sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/binding_validation_regression.py
+        sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/rereference_channels.py
     - name: Signal Processing Regression Cpp
       run: $GITHUB_WORKSPACE/cpp_package/examples/signal_processing/build/signal_processing_regression
       env:

.github/workflows/run_windows.yml

@@ -250,6 +250,7 @@ jobs:
         python python_package\examples\tests\decomposition_regression.py
         python python_package\examples\tests\signal_metrics_regression.py
         python python_package\examples\tests\binding_validation_regression.py
+        python python_package\examples\tests\rereference_channels.py
       shell: cmd
     - name: Signal Processing Regression Csharp
       run: dotnet run --project csharp_package\brainflow\examples\signal_processing_regression

cpp_package/examples/signal_processing/src/signal_processing_regression.cpp

@@ -71,5 +71,41 @@ int main ()
     lengths[0]++;
     reject (
         [&] { DataFilter::perform_inverse_wavelet_transform (coefficients, lengths, 256, 3, 3); });
+
+    BrainFlowArray<double, 2> reref_data (4, 3);
+    double values[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 3; j++)
+            reref_data.at (i, j) = values[i * 3 + j];
+    DataFilter::rereference_channels (reref_data, {0, 1}, {0, 1});
+    for (int j = 0; j < 3; j++)
+        if (std::abs (reref_data.at (0, j) + 1.5) > 1e-12 ||
+            std::abs (reref_data.at (1, j) - 1.5) > 1e-12 ||
+            std::abs (reref_data.at (2, j) - values[6 + j]) > 1e-12)
+            throw std::runtime_error ("Two-channel average reference failed");
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 3; j++)
+            reref_data.at (i, j) = values[i * 3 + j];
+    DataFilter::rereference_channels (reref_data, {0, 1, 2, 3}, {1});
+    for (int j = 0; j < 3; j++)
+        if (std::abs (reref_data.at (1, j)) > 1e-12 ||
+            std::abs (reref_data.at (0, j) - (values[j] - values[3 + j])) > 1e-12)
+            throw std::runtime_error ("Single-channel reference failed");
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 3; j++)
+            reref_data.at (i, j) = values[i * 3 + j];
+    DataFilter::rereference_channels (reref_data, {0, 1, 2, 3}, {0, 1, 2, 3});
+    for (int j = 0; j < 3; j++)
+    {
+        double mean = (values[j] + values[3 + j] + values[6 + j] + values[9 + j]) / 4.0;
+        for (int i = 0; i < 4; i++)
+            if (std::abs (reref_data.at (i, j) - (values[i * 3 + j] - mean)) > 1e-12)
+                throw std::runtime_error ("Common average reference failed");
+    }
+    reject ([&] { DataFilter::rereference_channels (reref_data, {0}, {4}); });
+    reject ([&] { DataFilter::rereference_channels (reref_data, {0, 0}, {1}); });
+
     std::cout << "Signal-processing binding validation passed." << std::endl;
 }

cpp_package/src/data_filter.cpp

@@ -670,6 +670,28 @@ double *DataFilter::get_activity_index (const double *accel_x, const double *acc
     return output;
 }
 
+void DataFilter::rereference_channels (BrainFlowArray<double, 2> &data,
+    const std::vector<int> &channels, const std::vector<int> &reference_channels)
+{
+    validate_channels (data, channels);
+    validate_channels (data, reference_channels);
+    if (channels.size () > (size_t)std::numeric_limits<int>::max () ||
+        reference_channels.size () > (size_t)std::numeric_limits<int>::max ())
+    {
+        throw BrainFlowException (
+            "invalid channel dimensions", (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+    }
+    std::vector<int> channels_copy (channels.begin (), channels.end ());
+    std::vector<int> reference_copy (reference_channels.begin (), reference_channels.end ());
+    int res = ::rereference_channels (data.get_raw_ptr (), data.get_size (0), data.get_size (1),
+        channels_copy.data (), (int)channels_copy.size (), reference_copy.data (),
+        (int)reference_copy.size ());
+    if (res != (int)BrainFlowExitCodes::STATUS_OK)
+    {
+        throw BrainFlowException ("failed to rereference channels", res);
+    }
+}
+
 std::string DataFilter::get_version ()
 {
     char version[64];

cpp_package/src/inc/data_filter.h

@@ -273,6 +273,16 @@ class DataFilter
         double noise_var_x = 0.0, double noise_var_y = 0.0, double noise_var_z = 0.0,
         int *output_len = NULL);
 
+    /**
+     * re-reference selected EEG/EXG channels in-place
+     * @param data 2d array with shape (channels, samples), modified in-place
+     * @param channels zero-based channel indices to update
+     * @param reference_channels zero-based channel indices used to compute the mean reference
+     * Overlap between channels and reference_channels is allowed (common average reference).
+     */
+    static void rereference_channels (BrainFlowArray<double, 2> &data,
+        const std::vector<int> &channels, const std::vector<int> &reference_channels);
+
     /// get brainflow version
     static std::string get_version ();
 };

csharp_package/brainflow/brainflow/data_filter.cs

@@ -794,6 +794,34 @@ public static Tuple<double[,], double[,], double[,], double[,]> perform_ica (dou
             return return_data;
         }
 
+        public static void rereference_channels (double[,] data, int[] channels, int[] reference_channels)
+        {
+            if (data == null || channels == null || reference_channels == null)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            int rows = data.GetLength (0);
+            int cols = data.GetLength (1);
+            if (rows <= 0 || cols <= 0 || (long)rows * cols > int.MaxValue)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            double[] data_1d = data.Flatten ();
+            int res = DataHandlerLibrary.rereference_channels (data_1d, rows, cols, channels, channels.Length,
+                    reference_channels, reference_channels.Length);
+            if (res != (int)BrainFlowExitCodes.STATUS_OK)
+            {
+                throw new BrainFlowError (res);
+            }
+            for (int i = 0; i < rows; i++)
+            {
+                for (int j = 0; j < cols; j++)
+                {
+                    data[i, j] = data_1d[i * cols + j];
+                }
+            }
+        }
+
         /// <summary>
         /// Calculate normalized mean powers for bands 2-4, 4-8, 8-13, 13-30, 30-45 Hz.
         /// Uses get_custom_band_powers preprocessing and minimum retained data length.

csharp_package/brainflow/brainflow/data_handler_library.cs

@@ -190,6 +190,8 @@ public static extern int perform_wavelet_denoising (double[] data, int data_len,
         [DllImport ("DataHandler", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
         public static extern int perform_ica (double[] data, int rows, int cols, int num_components, double[] w, double[] k, double[] a, double[] s);
         [DllImport ("DataHandler", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
+        public static extern int rereference_channels (double[] data, int rows, int cols, int[] channels, int channels_len, int[] reference_channels, int reference_channels_len);
+        [DllImport ("DataHandler", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
         public static extern int get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int data_len, int sampling_rate, int period, double noise_var_x, double noise_var_y, double noise_var_z, double[] output);
         // unsafe methods working with pointers
         [DllImport ("DataHandler", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
@@ -304,6 +306,8 @@ public static extern int perform_wavelet_denoising (double[] data, int data_len,
         [DllImport ("DataHandler32", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
         public static extern int perform_ica (double[] data, int rows, int cols, int num_components, double[] w, double[] k, double[] a, double[] s);
         [DllImport ("DataHandler32", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
+        public static extern int rereference_channels (double[] data, int rows, int cols, int[] channels, int channels_len, int[] reference_channels, int reference_channels_len);
+        [DllImport ("DataHandler32", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
         public static extern int get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int data_len, int sampling_rate, int period, double noise_var_x, double noise_var_y, double noise_var_z, double[] output);
         // unsafe methods working with pointers
         [DllImport ("DataHandler32", SetLastError = true, CallingConvention = CallingConvention.Cdecl)]
@@ -397,6 +401,19 @@ public static int perform_ica (double[] data, int rows, int cols, int num_compon
             return (int)BrainFlowExitCodes.GENERAL_ERROR;
         }
 
+        public static int rereference_channels (double[] data, int rows, int cols, int[] channels, int channels_len, int[] reference_channels, int reference_channels_len)
+        {
+            switch (PlatformHelper.get_library_environment ())
+            {
+                case LibraryEnvironment.x64:
+                    return DataHandlerLibrary64.rereference_channels (data, rows, cols, channels, channels_len, reference_channels, reference_channels_len);
+                case LibraryEnvironment.x86:
+                    return DataHandlerLibrary32.rereference_channels (data, rows, cols, channels, channels_len, reference_channels, reference_channels_len);
+            }
+
+            return (int)BrainFlowExitCodes.GENERAL_ERROR;
+        }
+
         public static int get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int data_len, int sampling_rate, int period, double noise_var_x, double noise_var_y, double noise_var_z, double[] output)
         {
             switch (PlatformHelper.get_library_environment ())

docs/UserAPI.rst

@@ -136,6 +136,11 @@ Data files accept rectangular tab- or comma-separated finite numeric values, inc
 a final line without a newline. Invalid tokens and insufficient destination capacity
 return errors. Writing uses enough decimal precision to round-trip doubles.
 
+``DataFilter.rereference_channels`` subtracts the element-wise mean of
+``reference_channels`` from the selected ``channels`` of a ``(channels, samples)``
+array in place. Channel indices are zero-based. Overlap between the two index lists
+is allowed and supports common average referencing.
+
 Python API Reference
 ----------------------
 

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

@@ -137,6 +137,9 @@ int get_heart_rate (double[] ppg_ir, double[] ppg_red, int len, int sampling_rat
         int perform_ica (double[] data, int rows, int cols, int num_components, double[] w, double[] k, double[] a,
                 double[] s);
 
+        int rereference_channels (double[] data, int rows, int cols, int[] channels, int channels_len,
+                int[] reference_channels, int reference_channels_len);
+
         int get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int data_len, int sampling_rate, int period,
                 double noise_var_x, double noise_var_y, double noise_var_z, double[] output);
 
@@ -988,6 +991,41 @@ public static List<double[][]> perform_ica (double[][] data, int num_components,
         return res;
     }
 
+    public static void rereference_channels (double[][] data, int[] channels, int[] reference_channels)
+            throws BrainFlowError
+    {
+        if (data == null || data.length == 0 || channels == null || reference_channels == null)
+        {
+            throw new BrainFlowError ("invalid args for rereference_channels",
+                    BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        int cols = data[0].length;
+        for (int i = 0; i < data.length; i++)
+        {
+            if (data[i] == null || data[i].length != cols)
+            {
+                throw new BrainFlowError ("invalid args for rereference_channels",
+                        BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+            }
+        }
+        int rows = data.length;
+        double[] data_1d = new double[checked_product (rows, cols)];
+        for (int i = 0; i < rows; i++)
+        {
+            System.arraycopy (data[i], 0, data_1d, i * cols, cols);
+        }
+        int ec = instance.rereference_channels (data_1d, rows, cols, channels, channels.length,
+                reference_channels, reference_channels.length);
+        if (ec != BrainFlowExitCode.STATUS_OK.get_code ())
+        {
+            throw new BrainFlowError ("Failed to rereference_channels", ec);
+        }
+        for (int i = 0; i < rows; i++)
+        {
+            System.arraycopy (data_1d, i * cols, data[i], 0, cols);
+        }
+    }
+
     /**
      * get PSD
      * 

julia_package/brainflow/src/data_filter.jl

@@ -540,6 +540,26 @@ function perform_ica(data, num_components::Integer)
     return perform_ica_select_channels(data, num_components, channels)
 end
 
+@brainflow_rethrow function rereference_channels(data, channels, reference_channels)
+    if !(data isa StridedMatrix{Float64})
+        throw(BrainFlowError("In-place processing requires a StridedMatrix{Float64}", Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    shape = size(data)
+    if !(0 < length(data) <= typemax(Cint)) || prod(shape) > typemax(Cint)
+        throw(BrainFlowError("Expected a nonempty real matrix", Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    channels = _dsp_channels(data, channels)
+    reference_channels = _dsp_channels(data, reference_channels)
+    channels_c = Int32.(channels .- 1)
+    reference_c = Int32.(reference_channels .- 1)
+    data_1d = copy(reshape(transpose(data), (1, prod(shape))))
+    ccall((:rereference_channels, DATA_HANDLER_INTERFACE), Cint,
+            (Ptr{Float64}, Cint, Cint, Ptr{Int32}, Cint, Ptr{Int32}, Cint),
+            data_1d, Int32(shape[1]), Int32(shape[2]), channels_c, length(channels_c), reference_c, length(reference_c))
+    copyto!(data, transpose(reshape(data_1d, shape[2], shape[1])))
+    return
+end
+
 @brainflow_rethrow function get_psd(data, sampling_rate::Integer, window::WinType)
     data = _dsp_signal(data)
 

matlab_package/brainflow/DataFilter.m

@@ -597,6 +597,22 @@ function write_file(data, file_name, file_mode)
             DataFilter.check_ec(exit_code, task_name);
             output = temp_output.Value;
         end
+
+        function data = rereference_channels(data, channels, reference_channels)
+            % re-reference selected channels in-place
+            channels = DataFilter.checked_channels(data, channels);
+            reference_channels = DataFilter.checked_channels(data, reference_channels);
+            task_name = 'rereference_channels';
+            transposed = transpose(data);
+            flatten = transpose(transposed(:));
+            temp_input = libpointer('doublePtr', flatten);
+            temp_channels = libpointer('int32Ptr', int32(channels - 1));
+            temp_reference = libpointer('int32Ptr', int32(reference_channels - 1));
+            lib_name = DataFilter.load_lib();
+            exit_code = calllib(lib_name, task_name, temp_input, size(data, 1), size(data, 2), temp_channels, numel(channels), temp_reference, numel(reference_channels));
+            DataFilter.check_ec(exit_code, task_name);
+            data = reshape(temp_input.Value, [size(data, 2), size(data, 1)])';
+        end
 
     end
 

nodejs_package/brainflow/data_filter.ts

@@ -66,6 +66,7 @@ class DataHandlerDLL extends DataHandlerFunctions
         this.detectPeaksZScore = this.lib.func(CLike.detect_peaks_z_score);
         this.performIca = this.lib.func(CLike.perform_ica);
         this.getActivityIndex = this.lib.func(CLike.get_activity_index);
+        this.rereferenceChannels = this.lib.func(CLike.rereference_channels);
         this.getCsp = this.lib.func(CLike.get_csp);
         this.detrend = this.lib.func(CLike.detrend);
         this.calcStddev = this.lib.func(CLike.calc_stddev);
@@ -783,4 +784,29 @@ export class DataFilter
         }
         return output;
     }
+
+    public static rereferenceChannels(
+        data: number[][], channels: number[], referenceChannels: number[]): void
+    {
+        const columns = DataFilter.checkMatrix(data, true);
+        DataFilter.checkChannels(data, channels, columns);
+        DataFilter.checkChannels(data, referenceChannels, columns);
+        const flat = data.flat();
+        const res = DataHandlerDLL.getInstance().rereferenceChannels(
+            flat, data.length, columns, channels, channels.length, referenceChannels,
+            referenceChannels.length);
+        if (res !== BrainFlowExitCodes.STATUS_OK)
+        {
+            throw new BrainFlowError (res, 'Could not rereference channels');
+        }
+        var i: number;
+        var j: number;
+        for (i = 0; i < data.length; i++)
+        {
+            for (j = 0; j < columns; j++)
+            {
+                data[i][j] = flat[i * columns + j];
+            }
+        }
+    }
 }

nodejs_package/brainflow/functions.types.ts

@@ -372,6 +372,8 @@ export enum DataHandlerCLikeFunctions {
         'int perform_ica (double *data, int rows, int cols, int num_components, _Inout_ double *w_mat, _Inout_ double *k_mat, _Inout_ double *a_mat, _Inout_ double *s_mat)',
     get_activity_index =
         'int get_activity_index (double *accel_x, double *accel_y, double *accel_z, int data_len, int sampling_rate, int period, double noise_var_x, double noise_var_y, double noise_var_z, _Inout_ double *activity_index)',
+    rereference_channels =
+        'int rereference_channels (_Inout_ double *data, int rows, int cols, int *channels, int channels_len, int *reference_channels, int reference_channels_len)',
     get_csp =
         'int get_csp (const double *data, const double *labels, int n_epochs, int n_channels, int n_times, _Inout_ double *output_w, _Inout_ double *output_d)',
     get_railed_percentage =
@@ -450,6 +452,8 @@ export class DataHandlerFunctions
     getActivityIndex!: (accelX: number[], accelY: number[], accelZ: number[], dataLen: number,
         samplingRate: number, period: number, noiseVarX: number, noiseVarY: number, noiseVarZ: number,
         activityIndex: number[]) => BrainFlowExitCodes;
+    rereferenceChannels!: (data: number[], rows: number, cols: number, channels: number[],
+        channelsLen: number, referenceChannels: number[], referenceChannelsLen: number) => BrainFlowExitCodes;
     getCsp!: (data: number[], labels: number[], nEpochs: number, nChannels: number, nTimes: number,
         outputW: number[], outputD: number[]) => BrainFlowExitCodes;
     detrend!: (rawData: number[], dataLen: number, detrendOperation: number) => BrainFlowExitCodes;