Add activity index calculation from accelerometer data (fixes #642)

.github/workflows/run_unix.yml

@@ -378,6 +378,8 @@ jobs:
       run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/band_power.py
     - name: BandPowerAll Python
       run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/band_power_all.py
+    - name: Activity Index Python
+      run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/activity_index.py
     - name: Denoising Cpp
       run: $GITHUB_WORKSPACE/cpp_package/examples/signal_processing/build/denoising
       env:

cpp_package/src/data_filter.cpp

@@ -585,6 +585,45 @@ double DataFilter::get_railed_percentage (double *data, int data_len, int gain)
     return output;
 }
 
+double *DataFilter::get_activity_index (const double *accel_x, const double *accel_y,
+    const double *accel_z, int data_len, int sampling_rate, int period,
+    double noise_var_x, double noise_var_y, double noise_var_z, int *output_len)
+{
+    if ((data_len <= 0) || (sampling_rate <= 0))
+    {
+        throw BrainFlowException ("invalid arguments for get_activity_index",
+            (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+    }
+    if (period <= 0)
+    {
+        period = data_len - (data_len % sampling_rate);
+    }
+    if ((period < sampling_rate) || (period > data_len) || (period % sampling_rate != 0))
+    {
+        throw BrainFlowException ("invalid period for get_activity_index",
+            (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+    }
+    int num_epochs = data_len / period;
+    if (num_epochs <= 0)
+    {
+        throw BrainFlowException ("data length shorter than epoch period",
+            (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR);
+    }
+    double *output = new double[num_epochs];
+    int res = ::get_activity_index (accel_x, accel_y, accel_z, data_len, sampling_rate, period,
+        noise_var_x, noise_var_y, noise_var_z, output);
+    if (res != (int)BrainFlowExitCodes::STATUS_OK)
+    {
+        delete[] output;
+        throw BrainFlowException ("unable to calculate activity index", res);
+    }
+    if (output_len != NULL)
+    {
+        *output_len = num_epochs;
+    }
+    return output;
+}
+
 std::string DataFilter::get_version ()
 {
     char version[64];

cpp_package/src/inc/data_filter.h

@@ -236,6 +236,24 @@ class DataFilter
         BrainFlowArray<double, 2>, BrainFlowArray<double, 2>>
     perform_ica (const BrainFlowArray<double, 2> &data, int num_components);
 
+    /**
+     * calculate activity index from 3-axis accelerometer data
+     * @param accel_x input 1d array
+     * @param accel_y input 1d array
+     * @param accel_z input 1d array
+     * @param data_len size of array
+     * @param sampling_rate sampling rate in Hz
+     * @param period epoch length in samples (defaults to full integer seconds if <= 0)
+     * @param noise_var_x baseline rest noise variance for X axis (default 0.0)
+     * @param noise_var_y baseline rest noise variance for Y axis (default 0.0)
+     * @param noise_var_z baseline rest noise variance for Z axis (default 0.0)
+     * @param output_len pointer to int to store number of epochs calculated
+     * @return pointer to array of activity indices
+     */
+    static double *get_activity_index (const double *accel_x, const double *accel_y,
+        const double *accel_z, int data_len, int sampling_rate, int period = 0,
+        double noise_var_x = 0.0, double noise_var_y = 0.0, double noise_var_z = 0.0,
+        int *output_len = NULL);
 
     /// get brainflow version
     static std::string get_version ();

csharp_package/brainflow/brainflow/data_filter.cs

@@ -1,4 +1,4 @@
-using brainflow.math;
+using brainflow.math;
 
 using System;
 using System.Numerics;
@@ -585,6 +585,53 @@ public static void write_file (double[,] data, string file_name, string file_mod
             return result;
         }
 
+        /// <summary>
+        /// calculate activity index from 3-axis accelerometer data using Bai et al. (2016) formulation
+        /// </summary>
+        public static double[] get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int sampling_rate, int period = 0, double noise_var_x = 0.0, double noise_var_y = 0.0, double noise_var_z = 0.0)
+        {
+            if (accel_x == null || accel_y == null || accel_z == null)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            if ((accel_x.Length != accel_y.Length) || (accel_x.Length != accel_z.Length))
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            if (accel_x.Length == 0)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            if (sampling_rate <= 0 || accel_x.Length < sampling_rate)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            if (period <= 0)
+            {
+                period = accel_x.Length - (accel_x.Length % sampling_rate);
+            }
+            if (period < sampling_rate || accel_x.Length < period || (period % sampling_rate != 0))
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            if (noise_var_x < 0.0 || noise_var_y < 0.0 || noise_var_z < 0.0)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            int num_epochs = accel_x.Length / period;
+            if (num_epochs == 0)
+            {
+                throw new BrainFlowError ((int)BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR);
+            }
+            double[] output = new double[num_epochs];
+            int res = DataHandlerLibrary.get_activity_index (accel_x, accel_y, accel_z, accel_x.Length, sampling_rate, period, noise_var_x, noise_var_y, noise_var_z, output);
+            if (res != (int)BrainFlowExitCodes.STATUS_OK)
+            {
+                throw new BrainFlowError (res);
+            }
+            return output;
+        }
+
         /// <summary>
         /// calculate nearest power of two
         /// </summary>

csharp_package/brainflow/brainflow/data_handler_library.cs

@@ -1,4 +1,4 @@
-using System.Runtime.InteropServices;
+using System.Runtime.InteropServices;
 
 namespace brainflow
 {
@@ -187,6 +187,8 @@ public static extern int perform_wavelet_denoising (double[] data, int data_len,
         public static extern int get_heart_rate (double[] ppg_ir, double[] ppg_red, int data_size, int sampling_rate, int fft_size, double[] output);
         [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 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)]
         public static unsafe extern int perform_lowpass (double* data, int len, int sampling_rate, double cutoff, int order, int filter_type, double ripple);
@@ -297,6 +299,8 @@ public static extern int perform_wavelet_denoising (double[] data, int data_len,
         public static extern int get_heart_rate (double[] ppg_ir, double[] ppg_red, int data_size, int sampling_rate, int fft_size, double[] output);
         [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 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)]
         public static unsafe extern int perform_lowpass (double* data, int len, int sampling_rate, double cutoff, int order, int filter_type, double ripple);
@@ -389,6 +393,19 @@ public static int perform_ica (double[] data, int rows, int cols, int num_compon
             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 ())
+            {
+                case LibraryEnvironment.x64:
+                    return DataHandlerLibrary64.get_activity_index (accel_x, accel_y, accel_z, data_len, sampling_rate, period, noise_var_x, noise_var_y, noise_var_z, output);
+                case LibraryEnvironment.x86:
+                    return DataHandlerLibrary32.get_activity_index (accel_x, accel_y, accel_z, data_len, sampling_rate, period, noise_var_x, noise_var_y, noise_var_z, output);
+            }
+
+            return (int)BrainFlowExitCodes.GENERAL_ERROR;
+        }
+
         public static int perform_lowpass (double[] data, int len, int sampling_rate, double cutoff, int order, int filter_type, double ripple)
         {
             switch (PlatformHelper.get_library_environment ())

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

@@ -104,6 +104,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 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);
+
         int get_version_data_handler (byte[] version, int[] len, int max_len);
 
         int log_message_data_handler (int log_level, String message);
@@ -1073,6 +1076,67 @@ public static double[][] read_file (String file_name) throws BrainFlowError
         return reshape_data_to_2d (num_rows[0], num_cols[0], data_arr);
     }
 
+    /**
+     * calculate activity index from 3-axis accelerometer data using Bai et al. (2016) formulation
+     */
+    public static double[] get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int sampling_rate, int period,
+            double noise_var_x, double noise_var_y, double noise_var_z) throws BrainFlowError
+    {
+        if (accel_x == null || accel_y == null || accel_z == null)
+        {
+            throw new BrainFlowError ("Null pointer passed", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        if ((accel_x.length != accel_y.length) || (accel_x.length != accel_z.length))
+        {
+            throw new BrainFlowError ("Array lengths do not match", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        if (accel_x.length == 0)
+        {
+            throw new BrainFlowError ("Input arrays must not be empty", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        if (sampling_rate <= 0 || accel_x.length < sampling_rate)
+        {
+            throw new BrainFlowError ("Invalid sampling rate or data length shorter than 1 second", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        if (period <= 0)
+        {
+            period = accel_x.length - (accel_x.length % sampling_rate);
+        }
+        if (period < sampling_rate || accel_x.length < period || (period % sampling_rate != 0))
+        {
+            throw new BrainFlowError ("Invalid period or data length is shorter than period", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        if (noise_var_x < 0 || noise_var_y < 0 || noise_var_z < 0)
+        {
+            throw new BrainFlowError ("Noise variances must be non-negative", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        int num_epochs = accel_x.length / period;
+        if (num_epochs == 0)
+        {
+            throw new BrainFlowError ("Data length is shorter than period", BrainFlowExitCode.INVALID_ARGUMENTS_ERROR.get_code ());
+        }
+        double[] output = new double[num_epochs];
+        int ec = instance.get_activity_index (accel_x, accel_y, accel_z, accel_x.length, sampling_rate, period,
+                noise_var_x, noise_var_y, noise_var_z, output);
+        if (ec != BrainFlowExitCode.STATUS_OK.get_code ())
+        {
+            throw new BrainFlowError ("Failed to calculate activity index", ec);
+        }
+        return output;
+    }
+
+    public static double[] get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int sampling_rate, int period)
+            throws BrainFlowError
+    {
+        return get_activity_index (accel_x, accel_y, accel_z, sampling_rate, period, 0.0, 0.0, 0.0);
+    }
+
+    public static double[] get_activity_index (double[] accel_x, double[] accel_y, double[] accel_z, int sampling_rate)
+            throws BrainFlowError
+    {
+        return get_activity_index (accel_x, accel_y, accel_z, sampling_rate, 0, 0.0, 0.0, 0.0);
+    }
+
     public static double[] reshape_data_to_1d (int num_rows, int num_cols, double[][] buf)
     {
         double[] output_buf = new double[num_rows * num_cols];

julia_package/brainflow/src/data_filter.jl

@@ -461,3 +461,33 @@ end
             psd[1], psd[2], length(psd[1]), Float64(freq_start), Float64(freq_end), band_power)
     return band_power[1]
 end
+
+@brainflow_rethrow function get_activity_index(accel_x, accel_y, accel_z, sampling_rate::Integer, period::Integer=0, noise_var_x::Real=0.0, noise_var_y::Real=0.0, noise_var_z::Real=0.0)
+    if (length(accel_x) != length(accel_y)) || (length(accel_x) != length(accel_z))
+        throw(BrainFlowError(string("Arrays lengths must match ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    data_len = length(accel_x)
+    if data_len == 0
+        throw(BrainFlowError(string("Input arrays must not be empty ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    if sampling_rate <= 0 || data_len < sampling_rate
+        throw(BrainFlowError(string("Invalid sampling rate or data length shorter than 1 second ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    if period <= 0
+        period = data_len - mod(data_len, sampling_rate)
+    end
+    if period < sampling_rate || data_len < period || mod(period, sampling_rate) != 0
+        throw(BrainFlowError(string("Invalid period or data length is shorter than period ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    if noise_var_x < 0.0 || noise_var_y < 0.0 || noise_var_z < 0.0
+        throw(BrainFlowError(string("Noise variances must be non-negative ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    num_epochs = div(data_len, period)
+    if num_epochs == 0
+        throw(BrainFlowError(string("Data length is shorter than period ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
+    end
+    output = Vector{Float64}(undef, num_epochs)
+    ccall((:get_activity_index, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Cint, Cint, Cint, Float64, Float64, Float64, Ptr{Float64}),
+            accel_x, accel_y, accel_z, Int32(data_len), Int32(sampling_rate), Int32(period), Float64(noise_var_x), Float64(noise_var_y), Float64(noise_var_z), output)
+    return output
+end

matlab_package/brainflow/DataFilter.m

@@ -460,6 +460,51 @@ function write_file(data, file_name, file_mode)
             data =  transpose(reshape(data_array.Value(1, 1:data_count.Value), [num_cols.Value, num_rows.value]));
         end
 
+        function output = get_activity_index(accel_x, accel_y, accel_z, sampling_rate, period, noise_var_x, noise_var_y, noise_var_z)
+            % calculate activity index using Bai et al. (2016) formulation
+            if nargin < 4
+                error('accel_x, accel_y, accel_z, and sampling_rate are required');
+            end
+            if isempty(accel_x) || isempty(accel_y) || isempty(accel_z)
+                error('Input arrays must not be empty');
+            end
+            if size(accel_x, 2) ~= size(accel_y, 2) || size(accel_x, 2) ~= size(accel_z, 2)
+                error('Length of accel_x, accel_y, and accel_z must match');
+            end
+            if floor(sampling_rate) ~= sampling_rate || sampling_rate <= 0
+                error('Sampling rate must be a positive integer');
+            end
+            data_len = size(accel_x, 2);
+            if data_len < sampling_rate
+                error('Data length must be at least one second');
+            end
+            if nargin < 5 || period <= 0
+                period = data_len - mod(data_len, sampling_rate);
+            end
+            if floor(period) ~= period
+                error('Period must be an integer');
+            end
+            if period < sampling_rate || period > data_len || mod(period, sampling_rate) ~= 0
+                error('Period must be an integer multiple of sampling rate and <= data_len');
+            end
+            if nargin < 6; noise_var_x = 0.0; end
+            if nargin < 7; noise_var_y = 0.0; end
+            if nargin < 8; noise_var_z = 0.0; end
+            if noise_var_x < 0 || noise_var_y < 0 || noise_var_z < 0
+                error('Noise variances must be non-negative');
+            end
+            task_name = 'get_activity_index';
+            temp_input_x = libpointer('doublePtr', accel_x);
+            temp_input_y = libpointer('doublePtr', accel_y);
+            temp_input_z = libpointer('doublePtr', accel_z);
+            lib_name = DataFilter.load_lib();
+            num_epochs = floor(data_len / period);
+            temp_output = libpointer('doublePtr', zeros(1, num_epochs));
+            exit_code = calllib(lib_name, task_name, temp_input_x, temp_input_y, temp_input_z, data_len, sampling_rate, period, noise_var_x, noise_var_y, noise_var_z, temp_output);
+            DataFilter.check_ec(exit_code, task_name);
+            output = temp_output.Value;
+        end
+
     end
 
 end