diff --git a/.github/workflows/run_unix.yml b/.github/workflows/run_unix.yml index f0d10f9d8..688da4702 100644 --- a/.github/workflows/run_unix.yml +++ b/.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: diff --git a/after_processing.png b/after_processing.png new file mode 100644 index 000000000..766373894 Binary files /dev/null and b/after_processing.png differ diff --git a/before_processing.png b/before_processing.png new file mode 100644 index 000000000..b7572dbb9 Binary files /dev/null and b/before_processing.png differ diff --git a/cpp_package/src/data_filter.cpp b/cpp_package/src/data_filter.cpp index 561fcc124..ecc571bb6 100644 --- a/cpp_package/src/data_filter.cpp +++ b/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]; diff --git a/cpp_package/src/inc/data_filter.h b/cpp_package/src/inc/data_filter.h index de66122de..ceaa2eafa 100644 --- a/cpp_package/src/inc/data_filter.h +++ b/cpp_package/src/inc/data_filter.h @@ -236,6 +236,24 @@ class DataFilter BrainFlowArray, BrainFlowArray> perform_ica (const BrainFlowArray &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 (); diff --git a/csharp_package/brainflow/brainflow/data_filter.cs b/csharp_package/brainflow/brainflow/data_filter.cs index 98fce81c2..54ad41d90 100644 --- a/csharp_package/brainflow/brainflow/data_filter.cs +++ b/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; } + /// + /// 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 = 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; + } + /// /// calculate nearest power of two /// diff --git a/csharp_package/brainflow/brainflow/data_handler_library.cs b/csharp_package/brainflow/brainflow/data_handler_library.cs index 067e7ad73..2107f5386 100644 --- a/csharp_package/brainflow/brainflow/data_handler_library.cs +++ b/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 ()) diff --git a/java_package/brainflow/src/main/java/brainflow/DataFilter.java b/java_package/brainflow/src/main/java/brainflow/DataFilter.java index 13bb57db3..1fb538730 100644 --- a/java_package/brainflow/src/main/java/brainflow/DataFilter.java +++ b/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]; diff --git a/julia_package/brainflow/src/data_filter.jl b/julia_package/brainflow/src/data_filter.jl index 70b7efa00..8e2baf974 100644 --- a/julia_package/brainflow/src/data_filter.jl +++ b/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 diff --git a/matlab_package/brainflow/DataFilter.m b/matlab_package/brainflow/DataFilter.m index aefa98f39..cf1c94f8f 100644 --- a/matlab_package/brainflow/DataFilter.m +++ b/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 \ No newline at end of file diff --git a/nodejs_package/brainflow/data_filter.ts b/nodejs_package/brainflow/data_filter.ts index b032298f4..05660cece 100644 --- a/nodejs_package/brainflow/data_filter.ts +++ b/nodejs_package/brainflow/data_filter.ts @@ -65,6 +65,7 @@ class DataHandlerDLL extends DataHandlerFunctions this.lib.func(CLike.restore_data_from_wavelet_detailed_coeffs); 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.getCsp = this.lib.func(CLike.get_csp); this.detrend = this.lib.func(CLike.detrend); this.calcStddev = this.lib.func(CLike.calc_stddev); @@ -607,4 +608,60 @@ export class DataFilter } return output[0]; } + + public static getActivityIndex( + accelX: number[], accelY: number[], accelZ: number[], + samplingRate: number, period: number = 0, + noiseVarX: number = 0, noiseVarY: number = 0, noiseVarZ: number = 0): number[] + { + if (accelX.length !== accelY.length || accelX.length !== accelZ.length) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "arrays lengths must match"); + } + if (accelX.length === 0) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "input arrays must not be empty"); + } + if (!Number.isInteger(samplingRate) || samplingRate <= 0 || accelX.length < samplingRate) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "invalid sampling rate or data shorter than 1 second"); + } + if (period <= 0) + { + period = accelX.length - (accelX.length % samplingRate); + } + if (!Number.isInteger(period)) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "period must be an integer"); + } + if (period < samplingRate || accelX.length < period || (period % samplingRate !== 0)) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "invalid period or data length shorter than period"); + } + if (noiseVarX < 0 || noiseVarY < 0 || noiseVarZ < 0) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "noise variances must be non-negative"); + } + const numEpochs = Math.trunc(accelX.length / period); + if (numEpochs === 0) + { + throw new BrainFlowError ( + BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR, "data length is shorter than period"); + } + const output = [...new Array (numEpochs).fill(0)]; + const res = DataHandlerDLL.getInstance().getActivityIndex( + accelX, accelY, accelZ, accelX.length, samplingRate, period, + noiseVarX, noiseVarY, noiseVarZ, output); + if (res !== BrainFlowExitCodes.STATUS_OK) + { + throw new BrainFlowError (res, 'Could not calc activity index'); + } + return output; + } } diff --git a/nodejs_package/brainflow/functions.types.ts b/nodejs_package/brainflow/functions.types.ts index 177b54aa9..9740931c2 100644 --- a/nodejs_package/brainflow/functions.types.ts +++ b/nodejs_package/brainflow/functions.types.ts @@ -370,6 +370,8 @@ export enum DataHandlerCLikeFunctions { 'int detect_peaks_z_score (double *data, int data_len, int lag, double threshold, double influence, _Inout_ double *output)', perform_ica = '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)', 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 = @@ -445,6 +447,8 @@ export class DataHandlerFunctions influence: number, output: number[]) => BrainFlowExitCodes; performIca!: (data: number[], rows: number, cols: number, numComponents: number, wMat: number[], kMat: number[], aMat: number[], sMat: number[]) => BrainFlowExitCodes; + getActivityIndex!: (accelX: number[], accelY: number[], accelZ: number[], dataLen: number, + period: number, activityIndex: 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; diff --git a/python_package/brainflow/data_filter.py b/python_package/brainflow/data_filter.py index 33608363f..10d11fb5b 100644 --- a/python_package/brainflow/data_filter.py +++ b/python_package/brainflow/data_filter.py @@ -536,6 +536,21 @@ def __init__(self): ndpointer(ctypes.c_double) ] + self.get_activity_index = self.lib.get_activity_index + self.get_activity_index.restype = ctypes.c_int + self.get_activity_index.argtypes = [ + ndpointer(ctypes.c_double), + ndpointer(ctypes.c_double), + ndpointer(ctypes.c_double), + ctypes.c_int, + ctypes.c_int, + ctypes.c_int, + ctypes.c_double, + ctypes.c_double, + ctypes.c_double, + ndpointer(ctypes.c_double) + ] + self.get_version_data_handler = self.lib.get_version_data_handler self.get_version_data_handler.restype = ctypes.c_int self.get_version_data_handler.argtypes = [ @@ -1293,6 +1308,67 @@ def perform_ifft(cls, data): return output + @classmethod + def get_activity_index( + cls, + accel_x, + accel_y, + accel_z, + sampling_rate: int, + period: int = 0, + noise_var_x: float = 0.0, + noise_var_y: float = 0.0, + noise_var_z: float = 0.0, + ): + """get activity index from 3-axis accelerometer data using the Bai et al. (2016) formulation + + :param accel_x: acceleration X data + :type accel_x: NDArray[Shape["*"], Float64] + :param accel_y: acceleration Y data + :type accel_y: NDArray[Shape["*"], Float64] + :param accel_z: acceleration Z data + :type accel_z: NDArray[Shape["*"], Float64] + :param sampling_rate: sampling rate of accelerometer in Hz + :type sampling_rate: int + :param period: epoch length in samples (defaults to full integer seconds of data if 0) + :type period: int + :param noise_var_x: baseline rest noise variance for X axis + :type noise_var_x: float + :param noise_var_y: baseline rest noise variance for Y axis + :type noise_var_y: float + :param noise_var_z: baseline rest noise variance for Z axis + :type noise_var_z: float + :return: activity index values + :rtype: NDArray[Shape["*"], Float64] + """ + check_memory_layout_row_major(accel_x, 1) + check_memory_layout_row_major(accel_y, 1) + check_memory_layout_row_major(accel_z, 1) + if not (accel_x.shape[0] == accel_y.shape[0] == accel_z.shape[0]): + raise BrainFlowError('invalid shapes', BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value) + data_len = accel_x.shape[0] + if data_len == 0: + raise BrainFlowError('input arrays must not be empty', BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value) + if sampling_rate <= 0 or data_len < sampling_rate: + raise BrainFlowError('invalid sampling rate or data shorter than 1 second', BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value) + if period <= 0: + period = data_len - (data_len % sampling_rate) + if period < sampling_rate or data_len < period or (period % sampling_rate != 0): + raise BrainFlowError('invalid period or data length shorter than period', BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value) + if noise_var_x < 0.0 or noise_var_y < 0.0 or noise_var_z < 0.0: + raise BrainFlowError('noise variances must be non-negative', BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value) + num_epochs = data_len // period + if num_epochs == 0: + raise BrainFlowError('data length is shorter than period', BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value) + output = numpy.zeros(num_epochs).astype(numpy.float64) + res = DataHandlerDLL.get_instance().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 != BrainFlowExitCodes.STATUS_OK.value: + raise BrainFlowError('unable to calculate activity index', res) + return output + @classmethod def get_nearest_power_of_two(cls, value: int) -> int: """calc nearest power of two diff --git a/python_package/examples/tests/activity_index.py b/python_package/examples/tests/activity_index.py new file mode 100644 index 000000000..2c50d46aa --- /dev/null +++ b/python_package/examples/tests/activity_index.py @@ -0,0 +1,163 @@ +import numpy as np +import sys +import os + +# add python_package to sys.path +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..'))) + +from brainflow.data_filter import DataFilter +from brainflow.exit_codes import BrainFlowError, BrainFlowExitCodes + + +def test_activity_index(): + print("Testing get_activity_index...") + fs = 50 + + # 1. Test Constant Signal (AI must be exactly 0.0) + accel_x = np.full(100, 1.0, dtype=np.float64) + accel_y = np.full(100, 2.0, dtype=np.float64) + accel_z = np.full(100, -3.0, dtype=np.float64) + + ai_constant = DataFilter.get_activity_index(accel_x, accel_y, accel_z, sampling_rate=fs) + print(f"Constant signal AI: {ai_constant}") + assert len(ai_constant) == 1 + assert np.isclose(ai_constant[0], 0.0), f"Expected 0.0, got {ai_constant[0]}" + + # 2. Test Stationary Noisy Signal with Baseline Noise Cancellation + # Stationary sensor produces noise variance. With baseline noise subtracted, AI should be 0.0. + np.random.seed(42) + noise_x = np.random.normal(0, 0.2, 200) # theoretical var ~ 0.04 + noise_y = np.random.normal(0, 0.2, 200) + noise_z = np.random.normal(0, 0.2, 200) + + sample_var_x = float(np.var(noise_x[:fs])) + sample_var_y = float(np.var(noise_y[:fs])) + sample_var_z = float(np.var(noise_z[:fs])) + + # If we pass the exact baseline noise variance of rest data, AI must be 0.0 + ai_rest = DataFilter.get_activity_index( + noise_x[:fs], noise_y[:fs], noise_z[:fs], + sampling_rate=fs, period=fs, + noise_var_x=sample_var_x, noise_var_y=sample_var_y, noise_var_z=sample_var_z + ) + print(f"Stationary noise cancelled AI: {ai_rest}") + assert np.isclose(ai_rest[0], 0.0) + + # If baseline noise exceeds signal variance, clamping to 0 must prevent negative values + ai_clamped = DataFilter.get_activity_index( + noise_x[:fs], noise_y[:fs], noise_z[:fs], + sampling_rate=fs, period=fs, + noise_var_x=1.0, noise_var_y=1.0, noise_var_z=1.0 + ) + assert np.isclose(ai_clamped[0], 0.0) + + # 3. Test Signal with Known Variance minus Noise Variance + # Sine wave with amplitude 1.0 has variance = 0.5 + t = np.linspace(0, 1.0, fs, endpoint=False) + sig_x = np.sin(2 * np.pi * 2 * t) # var = 0.5 + sig_y = np.cos(2 * np.pi * 2 * t) # var = 0.5 + sig_z = np.zeros(fs) # var = 0.0 + noise_val = 0.1 + + actual_var_x = float(np.var(sig_x)) + actual_var_y = float(np.var(sig_y)) + actual_var_z = float(np.var(sig_z)) + expected_ai = np.sqrt(max(0.0, ((actual_var_x - noise_val) + (actual_var_y - noise_val) + (actual_var_z - 0.0)) / 3.0)) + + ai_test = DataFilter.get_activity_index( + sig_x, sig_y, sig_z, sampling_rate=fs, period=fs, + noise_var_x=noise_val, noise_var_y=noise_val, noise_var_z=0.0 + ) + print(f"Known signal AI: {ai_test[0]}, Expected: {expected_ai}") + assert np.isclose(ai_test[0], expected_ai) + + # 4. Test Epoch Semantics (Bai et al. 2016 summing adjacent 1-second AIs) + # Second 1: low motion wave (amplitude 0.5, var = 0.125) + # Second 2: high motion wave (amplitude 2.0, var = 2.0) + t1 = np.linspace(0, 1.0, fs, endpoint=False) + s1_x = 0.5 * np.sin(2 * np.pi * 2 * t1) + s1_y = 0.5 * np.cos(2 * np.pi * 2 * t1) + s1_z = np.zeros(fs) + + t2 = np.linspace(0, 1.0, fs, endpoint=False) + s2_x = 2.0 * np.sin(2 * np.pi * 2 * t2) + s2_y = 2.0 * np.cos(2 * np.pi * 2 * t2) + s2_z = np.zeros(fs) + + # 1-second AIs individually + ai_sec1 = DataFilter.get_activity_index(s1_x, s1_y, s1_z, sampling_rate=fs, period=fs)[0] + ai_sec2 = DataFilter.get_activity_index(s2_x, s2_y, s2_z, sampling_rate=fs, period=fs)[0] + + # Combined 2-second signal + comb_x = np.concatenate([s1_x, s2_x]) + comb_y = np.concatenate([s1_y, s2_y]) + comb_z = np.concatenate([s1_z, s2_z]) + + # 2-second epoch AI should be EXACTLY ai_sec1 + ai_sec2 per Bai et al. 2016 + ai_2sec = DataFilter.get_activity_index(comb_x, comb_y, comb_z, sampling_rate=fs, period=2 * fs) + assert len(ai_2sec) == 1 + print(f"2-second epoch AI: {ai_2sec[0]}, Sum of 1s AIs: {ai_sec1 + ai_sec2}") + assert np.isclose(ai_2sec[0], ai_sec1 + ai_sec2) + + # Show that it is NOT simply the variance recomputed over the entire 2-second period + recomputed_full_var_ai = np.sqrt((np.var(comb_x) + np.var(comb_y) + np.var(comb_z)) / 3.0) + print(f"Recomputed full variance AI: {recomputed_full_var_ai}") + assert not np.isclose(ai_2sec[0], recomputed_full_var_ai) + + # 5. Test Multi-epoch AI calculation (4 seconds -> two 2-second epochs) + comb4_x = np.concatenate([comb_x, comb_x]) + comb4_y = np.concatenate([comb_y, comb_y]) + comb4_z = np.concatenate([comb_z, comb_z]) + ai_4sec = DataFilter.get_activity_index(comb4_x, comb4_y, comb4_z, sampling_rate=fs, period=2 * fs) + assert len(ai_4sec) == 2 + assert np.isclose(ai_4sec[0], ai_2sec[0]) + assert np.isclose(ai_4sec[1], ai_2sec[0]) + + # 6. Test Invalid Arguments + # Shape mismatch + try: + DataFilter.get_activity_index(np.zeros(10), np.zeros(5), np.zeros(10), sampling_rate=fs) + assert False, "Should have raised BrainFlowError for shape mismatch" + except BrainFlowError as e: + assert e.exit_code == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value + + # Empty inputs + try: + DataFilter.get_activity_index(np.array([]), np.array([]), np.array([]), sampling_rate=fs) + assert False, "Should have raised BrainFlowError for empty inputs" + except BrainFlowError as e: + assert e.exit_code == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value + + # Invalid sampling rate + try: + DataFilter.get_activity_index(np.zeros(100), np.zeros(100), np.zeros(100), sampling_rate=0) + assert False, "Should have raised BrainFlowError for zero sampling rate" + except BrainFlowError as e: + assert e.exit_code == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value + + # Period greater than data length + try: + DataFilter.get_activity_index(np.zeros(100), np.zeros(100), np.zeros(100), sampling_rate=fs, period=200) + assert False, "Should have raised BrainFlowError for period > data_len" + except BrainFlowError as e: + assert e.exit_code == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value + + # Period not an integer multiple of sampling rate + try: + DataFilter.get_activity_index(np.zeros(100), np.zeros(100), np.zeros(100), sampling_rate=fs, period=75) + assert False, "Should have raised BrainFlowError for period not multiple of fs" + except BrainFlowError as e: + assert e.exit_code == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value + + # Negative noise variance + try: + DataFilter.get_activity_index(np.zeros(100), np.zeros(100), np.zeros(100), sampling_rate=fs, noise_var_x=-0.5) + assert False, "Should have raised BrainFlowError for negative noise variance" + except BrainFlowError as e: + assert e.exit_code == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value + + print("All activity index tests passed successfully!") + + +if __name__ == '__main__': + test_activity_index() diff --git a/rust_package/brainflow/src/data_filter.rs b/rust_package/brainflow/src/data_filter.rs index 8bdb06ee6..09a14c3e0 100644 --- a/rust_package/brainflow/src/data_filter.rs +++ b/rust_package/brainflow/src/data_filter.rs @@ -859,6 +859,55 @@ where Ok(check_brainflow_exit_code(res)?) } +/// Calculate activity index from 3-axis accelerometer data using the Bai et al. (2016) formulation. +pub fn get_activity_index( + accel_x: &[f64], + accel_y: &[f64], + accel_z: &[f64], + sampling_rate: usize, + period: Option, + noise_var: Option<(f64, f64, f64)>, +) -> Result> { + if accel_x.len() != accel_y.len() || accel_x.len() != accel_z.len() { + return Err(Error::BrainFlowError(BrainFlowError::InvalidArgumentsError)); + } + let data_len = accel_x.len(); + if data_len == 0 || sampling_rate == 0 || data_len < sampling_rate { + return Err(Error::BrainFlowError(BrainFlowError::InvalidArgumentsError)); + } + let period = period.unwrap_or(data_len - (data_len % sampling_rate)); + if period < sampling_rate || period > data_len || period % sampling_rate != 0 { + return Err(Error::BrainFlowError(BrainFlowError::InvalidArgumentsError)); + } + let num_epochs = data_len / period; + if num_epochs == 0 { + return Err(Error::BrainFlowError(BrainFlowError::InvalidArgumentsError)); + } + let (noise_var_x, noise_var_y, noise_var_z) = noise_var.unwrap_or((0.0, 0.0, 0.0)); + if noise_var_x < 0.0 || noise_var_y < 0.0 || noise_var_z < 0.0 { + return Err(Error::BrainFlowError(BrainFlowError::InvalidArgumentsError)); + } + let mut output = Vec::::with_capacity(num_epochs); + let res = unsafe { + data_handler::get_activity_index( + accel_x.as_ptr() as *const c_double, + accel_y.as_ptr() as *const c_double, + accel_z.as_ptr() as *const c_double, + data_len as c_int, + sampling_rate as c_int, + period as c_int, + noise_var_x, + noise_var_y, + noise_var_z, + output.as_mut_ptr() as *mut c_double, + ) + }; + check_brainflow_exit_code(res)?; + + unsafe { output.set_len(num_epochs) }; + Ok(output) +} + /// Get DataFilter version. pub fn get_version() -> Result { const MAX_CHARS: usize = 64; diff --git a/rust_package/brainflow/src/ffi/data_handler.rs b/rust_package/brainflow/src/ffi/data_handler.rs index e466eefa0..d66ef000d 100644 --- a/rust_package/brainflow/src/ffi/data_handler.rs +++ b/rust_package/brainflow/src/ffi/data_handler.rs @@ -272,6 +272,20 @@ extern "C" { s_mat: *mut f64, ) -> ::std::os::raw::c_int; } +extern "C" { + pub fn get_activity_index( + accel_x: *const f64, + accel_y: *const f64, + accel_z: *const f64, + data_len: ::std::os::raw::c_int, + sampling_rate: ::std::os::raw::c_int, + period: ::std::os::raw::c_int, + noise_var_x: f64, + noise_var_y: f64, + noise_var_z: f64, + activity_index: *mut f64, + ) -> ::std::os::raw::c_int; +} extern "C" { pub fn set_log_level_data_handler(log_level: ::std::os::raw::c_int) -> ::std::os::raw::c_int; } diff --git a/src/data_handler/data_handler.cpp b/src/data_handler/data_handler.cpp index 2f177faad..5d6294a3c 100644 --- a/src/data_handler/data_handler.cpp +++ b/src/data_handler/data_handler.cpp @@ -1726,6 +1726,71 @@ int perform_ica (double *data, int rows, int cols, int num_components, double *w return res; } +int 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, double *output) +{ + if ((accel_x == NULL) || (accel_y == NULL) || (accel_z == NULL) || (output == NULL) || + (data_len <= 0) || (sampling_rate <= 0) || (period < sampling_rate) || + (data_len < period) || (period % sampling_rate != 0) || (noise_var_x < 0.0) || + (noise_var_y < 0.0) || (noise_var_z < 0.0)) + { + data_logger->error ("Invalid arguments for get_activity_index: accel_x {}, accel_y {}, " + "accel_z {}, output {}, data_len {}, sampling_rate {}, period {}, " + "noise_var_x {}, noise_var_y {}, noise_var_z {}", + (accel_x != NULL), (accel_y != NULL), (accel_z != NULL), (output != NULL), data_len, + sampling_rate, period, noise_var_x, noise_var_y, noise_var_z); + return (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR; + } + + int sec_per_epoch = period / sampling_rate; + int num_epochs = data_len / period; + + for (int epoch = 0; epoch < num_epochs; epoch++) + { + double epoch_ai = 0.0; + int epoch_start_sec = epoch * sec_per_epoch; + for (int sec = 0; sec < sec_per_epoch; sec++) + { + int sec_idx = epoch_start_sec + sec; + int start_pos = sec_idx * sampling_rate; + int end_pos = start_pos + sampling_rate; + + double mean_x = 0.0, mean_y = 0.0, mean_z = 0.0; + for (int i = start_pos; i < end_pos; i++) + { + mean_x += accel_x[i]; + mean_y += accel_y[i]; + mean_z += accel_z[i]; + } + mean_x /= sampling_rate; + mean_y /= sampling_rate; + mean_z /= sampling_rate; + + double var_x = 0.0, var_y = 0.0, var_z = 0.0; + for (int i = start_pos; i < end_pos; i++) + { + var_x += (accel_x[i] - mean_x) * (accel_x[i] - mean_x); + var_y += (accel_y[i] - mean_y) * (accel_y[i] - mean_y); + var_z += (accel_z[i] - mean_z) * (accel_z[i] - mean_z); + } + var_x /= sampling_rate; + var_y /= sampling_rate; + var_z /= sampling_rate; + + double adj_x = var_x - noise_var_x; + double adj_y = var_y - noise_var_y; + double adj_z = var_z - noise_var_z; + double mean_adj = (adj_x + adj_y + adj_z) / 3.0; + double ai_1sec = sqrt (std::max (0.0, mean_adj)); + epoch_ai += ai_1sec; + } + output[epoch] = epoch_ai; + } + + return (int)BrainFlowExitCodes::STATUS_OK; +} + int get_version_data_handler (char *version, int *num_chars, int max_chars) { strncpy (version, BRAINFLOW_VERSION_STRING, max_chars); diff --git a/src/data_handler/inc/data_handler.h b/src/data_handler/inc/data_handler.h index 2e50e284c..517623de2 100644 --- a/src/data_handler/inc/data_handler.h +++ b/src/data_handler/inc/data_handler.h @@ -69,6 +69,9 @@ extern "C" double *data, int data_len, int lag, double threshold, double influence, double *output); SHARED_EXPORT int CALLING_CONVENTION perform_ica (double *data, int rows, int cols, int num_components, double *w_mat, double *k_mat, double *a_mat, double *s_mat); + SHARED_EXPORT int CALLING_CONVENTION 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, double *activity_index); // logging methods SHARED_EXPORT int CALLING_CONVENTION set_log_level_data_handler (int log_level); diff --git a/swift_package/Sources/BrainFlow/DataFilter.swift b/swift_package/Sources/BrainFlow/DataFilter.swift index ec21ef5a9..89008de40 100644 --- a/swift_package/Sources/BrainFlow/DataFilter.swift +++ b/swift_package/Sources/BrainFlow/DataFilter.swift @@ -704,6 +704,44 @@ public enum DataFilter { BrainFlowArray.reshape_data_to_2d(num_rows: num_rows, num_cols: num_cols, linear_buffer: linear_buffer) } + public static func get_activity_index( + accel_x: [Double], + accel_y: [Double], + accel_z: [Double], + sampling_rate: Int, + period: Int = 0, + noise_var_x: Double = 0.0, + noise_var_y: Double = 0.0, + noise_var_z: Double = 0.0 + ) throws -> [Double] { + guard accel_x.count == accel_y.count, accel_x.count == accel_z.count else { throw invalidArguments("Array lengths must match") } + let dataLen = accel_x.count + guard dataLen > 0 else { throw invalidArguments("Input arrays must not be empty") } + guard sampling_rate > 0, dataLen >= sampling_rate else { throw invalidArguments("Invalid sampling rate or data length shorter than 1 second") } + let periodToUse = period <= 0 ? (dataLen - (dataLen % sampling_rate)) : period + guard periodToUse >= sampling_rate, dataLen >= periodToUse, periodToUse % sampling_rate == 0 else { + throw invalidArguments("Invalid period or data length is shorter than period") + } + guard noise_var_x >= 0.0, noise_var_y >= 0.0, noise_var_z >= 0.0 else { + throw invalidArguments("Noise variances must be non-negative") + } + let numEpochs = dataLen / periodToUse + guard numEpochs > 0 else { throw invalidArguments("Data length is shorter than period") } + var output = [Double](repeating: 0.0, count: numEpochs) + try accel_x.withUnsafeBufferPointer { xPtr in + try accel_y.withUnsafeBufferPointer { yPtr in + try accel_z.withUnsafeBufferPointer { zPtr in + try output.withUnsafeMutableBufferPointer { outPtr in + try DataFilterNative.withData { native in + try checkBrainFlowExitCode(native.get_activity_index(xPtr.baseAddress, yPtr.baseAddress, zPtr.baseAddress, CInt(dataLen), CInt(sampling_rate), CInt(periodToUse), noise_var_x, noise_var_y, noise_var_z, outPtr.baseAddress), "Failed to calculate activity index") + } + } + } + } + } + return output + } + private static func withMutableData(_ data: inout [Double], _ body: (UnsafeMutablePointer?, Int) throws -> T) throws -> T { try data.withUnsafeMutableBufferPointer { pointer in try body(pointer.baseAddress, pointer.count) @@ -753,6 +791,7 @@ final class DataFilterNative { let restore_data_from_wavelet_detailed_coeffs: @convention(c) (UnsafeMutablePointer?, CInt, CInt, CInt, CInt, UnsafeMutablePointer?) -> CInt let detect_peaks_z_score: @convention(c) (UnsafeMutablePointer?, CInt, CInt, Double, Double, UnsafeMutablePointer?) -> CInt let perform_ica: @convention(c) (UnsafeMutablePointer?, CInt, CInt, CInt, UnsafeMutablePointer?, UnsafeMutablePointer?, UnsafeMutablePointer?, UnsafeMutablePointer?) -> CInt + let get_activity_index: @convention(c) (UnsafePointer?, UnsafePointer?, UnsafePointer?, CInt, CInt, CInt, Double, Double, Double, UnsafeMutablePointer?) -> CInt let set_log_level_data_handler: @convention(c) (CInt) -> CInt let set_log_file_data_handler: @convention(c) (UnsafePointer?) -> CInt let log_message_data_handler: @convention(c) (CInt, UnsafeMutablePointer?) -> CInt @@ -805,6 +844,7 @@ final class DataFilterNative { restore_data_from_wavelet_detailed_coeffs = try library.symbol("restore_data_from_wavelet_detailed_coeffs", as: type(of: restore_data_from_wavelet_detailed_coeffs)) detect_peaks_z_score = try library.symbol("detect_peaks_z_score", as: type(of: detect_peaks_z_score)) perform_ica = try library.symbol("perform_ica", as: type(of: perform_ica)) + get_activity_index = try library.symbol("get_activity_index", as: type(of: get_activity_index)) set_log_level_data_handler = try library.symbol("set_log_level_data_handler", as: type(of: set_log_level_data_handler)) set_log_file_data_handler = try library.symbol("set_log_file_data_handler", as: type(of: set_log_file_data_handler)) log_message_data_handler = try library.symbol("log_message_data_handler", as: type(of: log_message_data_handler))