statistical_features.py

# -*- coding: utf-8 -*-
"""Statistical_features.ipynb

Automatically generated by Colaboratory.

"""

import pandas as pd
import numpy as np
import cupy as cp
from scipy.stats import skew, kurtosis, moment, gmean, hmean, trim_mean
from statsmodels.tsa.stattools import acf
import torch
import time

class StatisticalFeatures:
  def __init__(self,
              window_size,
              n_lags_auto_correlation=None,
              # moment_orders=None,
              trimmed_mean_thresholds=None,
              # higuchi_k_values=None,
              # tsallis_q_parameer=1,
              # renyi_alpha_parameter=2,
              # permutation_entropy_order=3,
              # permutation_entropy_delay=1,
              # svd_entropy_order=3,
              # svd_entropy_delay=1,
              ):

    self.window_size = window_size
    # self.tsallis_q_parameer = tsallis_q_parameer
    # self.renyi_alpha_parameter = renyi_alpha_parameter
    # self.permutation_entropy_order = permutation_entropy_order
    # self.permutation_entropy_delay = permutation_entropy_delay
    # self.svd_entropy_order = svd_entropy_order
    # self.svd_entropy_delay = svd_entropy_delay

    if n_lags_auto_correlation is None:
        self.n_lags_auto_correlation = int(min(10 * np.log10(window_size), window_size - 1))
    else:
        self.n_lags_auto_correlation = n_lags_auto_correlation

    # if moment_orders is None:
    #     self.moment_orders = [3, 4]
    # else:
    #     self.moment_orders = moment_orders

    if trimmed_mean_thresholds is None:
        self.trimmed_mean_thresholds = [0.1, 0.15, 0.2, 0.25, 0.3]
    else:
        self.trimmed_mean_thresholds = trimmed_mean_thresholds

    # if higuchi_k_values is None:
    #     self.higuchi_k_values = list({5, 10, 20, window_size // 5})
    # else:
    #     self.higuchi_k_values = list(higuchi_k_values)


  def calculate_mean(self, signal_gpu, signal_cpu):
    name = 'mean'
    start_time = time.time()
    mean = torch.mean(signal_gpu)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    mean_gpu = np.array([mean.cpu().numpy()]) # just convert to numpy array, but it is still calculated by GPU
    mean_cpu = np.array([np.mean(signal_cpu)]) # the feature compute by CPU/numpy
    return name, mean_gpu, mean_cpu, gpu_extraction_time

  def calculate_geometric_mean(self, signal_gpu, signal_cpu):
    name = 'geometric mean'
    start_time = time.time()
    log_mean = torch.mean(torch.log(signal_gpu))
    geometric_mean = torch.exp(log_mean)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    geometric_mean_gpu = np.array([geometric_mean.cpu().numpy()])
    geometric_mean_cpu = np.array([gmean(signal_cpu)])
    return name, geometric_mean_gpu, geometric_mean_cpu, gpu_extraction_time

  def calculate_harmonic_mean(self, signal_gpu, signal_cpu):
    name = 'harmonic mean'
    start_time = time.time()
    reciprocal = 1.0 / signal_gpu
    harmonic_mean = 1.0 / torch.mean(reciprocal)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    harmonic_mean_gpu = np.array([harmonic_mean.cpu().numpy()])
    harmonic_mean_cpu = np.array([hmean(signal_cpu)])
    return name, harmonic_mean_gpu, harmonic_mean_cpu, gpu_extraction_time

  def calculate_trimmed_means(self, signal_gpu, signal_cpu):
    # there are small discrepancies between values compute by GPU and CPU
    # probably because the way of sort in GPU computation cause it
    # in this method, we use cp to sort instead of torch because torch sort is not stable and should convert to cpu compute
    name = []
    trimmed_means_gpu = []
    trimmed_means_cpu = []
    extraction_times = []
    signal_gpu_cp = cp.array(signal_gpu)  # convert cpu signal to cp array, so that can be compute by GPU
    for proportiontocut in self.trimmed_mean_thresholds:
      name.append(f'trimmed_means_thresh_{str(proportiontocut)}')
      start_time = time.time()
      num_elements = len(signal_gpu_cp)
      trim_amount = int(num_elements * proportiontocut / 100)
      sorted_data = cp.sort(signal_gpu_cp)
      #sorted_data  = torch.sort(signal_gpu.cpu()).value
      trimmed_data = sorted_data[trim_amount:-trim_amount]
      trimmed_mean = cp.mean(trimmed_data)
      #trimmed_mean = torch.mean(trimmed_data)
      end_time = time.time()
      gpu_extraction_time = end_time - start_time
      trimmed_means_gpu.append(trimmed_mean.get())  # cp use get() instead of item()
      extraction_times.append(gpu_extraction_time)
    for proportiontocut in self.trimmed_mean_thresholds:
            trimmed_means_cpu.append(trim_mean(signal_cpu, proportiontocut=proportiontocut))
    return name, trimmed_means_gpu, trimmed_means_cpu, extraction_times

  def calculate_mean_abs(self, signal_gpu, signal_cpu):
    name = 'mean abs'
    start_time = time.time()
    mean_abs = torch.mean(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    mean_abs_gpu = np.array([mean_abs.cpu().numpy()])
    mean_abs_cpu = np.array([np.mean(np.abs(signal_cpu))])
    return name, mean_abs_gpu, mean_abs_cpu, gpu_extraction_time

  def calculate_geometric_mean_abs(self, signal_gpu, signal_cpu):
    name = 'geometric mean abs'
    start_time = time.time()
    log_mean = torch.mean(torch.log(torch.abs(signal_gpu)))
    geometric_mean_abs = torch.exp(log_mean)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    geometric_mean_abs_gpu = np.array([geometric_mean_abs.cpu().numpy()])
    geometric_mean_abs_cpu = np.array([gmean(np.abs(signal_cpu))])
    return name, geometric_mean_abs_gpu, geometric_mean_abs_cpu, gpu_extraction_time

  def calculate_harmonic_mean_abs(self, signal_gpu, signal_cpu):
    name = 'harmonic mean abs'
    start_time = time.time()
    reciprocal = 1.0 / torch.abs(signal_gpu)
    harmonic_mean_abs = 1.0 / torch.mean(reciprocal)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    harmonic_mean_abs_gpu = np.array([harmonic_mean_abs.cpu().numpy()])
    harmonic_mean_abs_cpu = np.array([hmean(np.abs(signal_cpu))])
    return name, harmonic_mean_abs_gpu, harmonic_mean_abs_cpu, gpu_extraction_time

  def calculate_trimmed_means_abs(self, signal_gpu, signal_cpu):
    name = []
    trimmed_means_abs_gpu = []
    trimmed_means_abs_cpu = []
    extraction_times = []
    signal_gpu = torch.abs(signal_gpu)
    signal_gpu_cp = cp.array(signal_gpu)  # convert cpu signal to cp array, so that can be compute by GPU
    for proportiontocut in self.trimmed_mean_thresholds:
      name.append(f'trimmed_means_abs_thresh_{str(proportiontocut)}')
      start_time = time.time()
      num_elements = len(signal_gpu_cp)
      trim_amount = int(num_elements * proportiontocut / 100)
      sorted_data = cp.sort(signal_gpu_cp)
      #sorted_data  = torch.sort(signal_gpu.cpu()).value
      trimmed_data = sorted_data[trim_amount:-trim_amount]
      trimmed_mean_abs = cp.mean(trimmed_data)
      #trimmed_mean = torch.mean(trimmed_data)
      end_time = time.time()
      gpu_extraction_time = end_time - start_time
      trimmed_means_abs_gpu.append(trimmed_mean_abs.get())  # cp use get() instead of item()
      extraction_times.append(gpu_extraction_time)
    for proportiontocut in self.trimmed_mean_thresholds:
            trimmed_means_abs_cpu.append(trim_mean(np.abs(signal_cpu), proportiontocut=proportiontocut))
    return name, trimmed_means_abs_gpu, trimmed_means_abs_cpu, extraction_times

  def calculate_std(self, signal_gpu, signal_cpu):
    name = 'standard deviation'
    start_time = time.time()
    std_ = torch.std(signal_gpu)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    std_gpu = np.array([std_.cpu().numpy()])
    std_cpu = np.array([np.std(signal_cpu)])
    return name, std_gpu, std_cpu, gpu_extraction_time

  def calculate_std_abs(self, signal_gpu, signal_cpu):
    name = 'standard deviation abs'
    start_time = time.time()
    std_abs = torch.std(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    std_abs_gpu = np.array([std_abs.cpu().numpy()])
    std_abs_cpu = np.array([np.std(np.abs(signal_cpu))])
    return name, std_abs_gpu, std_abs_cpu, gpu_extraction_time

  def calculate_skewness(self, signal_gpu, signal_cpu):
    name = 'skewness'
    start_time = time.time()
    mean_gpu = torch.mean(signal_gpu)
    std_dev_gpu = torch.std(signal_gpu)
    n = len(signal_gpu)
    skewness_ = ((n / ((n - 1) * (n - 2))) *
                  torch.sum(((signal_gpu - mean_gpu) / std_dev_gpu) ** 3))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    skewness_gpu = np.array([skewness_.cpu().numpy()])
    skewness_cpu = np.array([skew(signal_cpu)])
    return name, skewness_gpu, skewness_cpu, gpu_extraction_time

  def calculate_skewness_abs(self, signal_gpu, signal_cpu):
    name = 'skewness abs'
    start_time = time.time()
    mean_gpu = torch.mean(torch.abs(signal_gpu))
    std_dev_gpu = torch.std(torch.abs(signal_gpu))
    n = len(signal_gpu)
    skewness_abs = ((n / ((n - 1) * (n - 2))) *
                    torch.sum(((signal_gpu - mean_gpu) / std_dev_gpu) ** 3))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    skewness_abs_gpu = np.array([skewness_abs.cpu().numpy()])
    skewness_abs_cpu = np.array([skew(np.abs(signal_cpu))])
    return name, skewness_abs_gpu, skewness_abs_cpu, gpu_extraction_time

  def calculate_kurtosis(self, signal_gpu, signal_cpu):
    name = 'kurtosis'
    start_time = time.time()
    mean_gpu = torch.mean(signal_gpu)
    std_dev_gpu = torch.std(signal_gpu)
    n = len(signal_gpu)
    kurtosis_ = ((n * (n + 1) / ((n - 1) * (n - 2) * (n - 3))) *
                torch.sum(((signal_gpu - mean_gpu) / std_dev_gpu) ** 4)) - (3 * (n - 1) ** 2 / ((n - 2) * (n - 3)))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    kurtosis_gpu = np.array([kurtosis_.cpu().numpy()])
    kurtosis_cpu = np.array([kurtosis(signal_cpu)])
    return name, kurtosis_gpu, kurtosis_cpu, gpu_extraction_time

  def calculate_kurtosis_abs(self, signal_gpu, signal_cpu):
    name = 'kurtosis abs'
    start_time = time.time()
    mean_gpu = torch.mean(torch.abs(signal_gpu))
    std_dev_gpu = torch.std(torch.abs(signal_gpu))
    n = len(signal_gpu)
    kurtosis_abs = ((n * (n + 1) / ((n - 1) * (n - 2) * (n - 3))) *
                torch.sum(((signal_gpu - mean_gpu) / std_dev_gpu) ** 4)) - (3 * (n - 1) ** 2 / ((n - 2) * (n - 3)))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    kurtosis_abs_gpu = np.array([kurtosis_abs.cpu().numpy()])
    kurtosis_abs_cpu = np.array([kurtosis(np.abs(signal_cpu))])
    return name, kurtosis_abs_gpu, kurtosis_abs_cpu, gpu_extraction_time

  def calculate_median(self, signal_gpu, signal_cpu):
    name = 'median'
    start_time = time.time()
    median_ = torch.median(signal_gpu)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    median_gpu = np.array([median_.cpu().numpy()])
    median_cpu = np.array([np.median(signal_cpu)])
    return name, median_gpu, median_cpu, gpu_extraction_time

  def calculate_median_abs(self, signal_gpu, signal_cpu):
    name = 'median abs'
    start_time = time.time()
    median_abs = torch.median(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    median_abs_gpu = np.array([median_abs.cpu().numpy()])
    median_abs_cpu = np.array([np.median(np.abs(signal_cpu))])
    return name, median_abs_gpu, median_abs_cpu, gpu_extraction_time

  def calculate_max(self, signal_gpu, signal_cpu):
    name = 'max'
    start_time = time.time()
    max_ = torch.max(signal_gpu)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    max_gpu = np.array([max_.cpu().numpy()])
    max_cpu = np.array([np.max(signal_cpu)])
    return name, max_gpu, max_cpu, gpu_extraction_time

  def calculate_max_abs(self, signal_gpu, signal_cpu):
    name = 'max abs'
    start_time = time.time()
    max_abs = torch.max(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    max_abs_gpu = np.array([max_abs.cpu().numpy()])
    max_abs_cpu = np.array([np.max(np.abs(signal_cpu))])
    return name, max_abs_gpu, max_abs_cpu, gpu_extraction_time

  def calculate_min(self, signal_gpu, signal_cpu):
    name = 'min'
    start_time = time.time()
    min_ = torch.min(signal_gpu)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    min_gpu = np.array([min_.cpu().numpy()])
    min_cpu = np.array([np.min(signal_cpu)])
    return name, min_gpu, min_cpu, gpu_extraction_time

  def calculate_min_abs(self, signal_gpu, signal_cpu):
    name = 'min abs'
    start_time = time.time()
    min_abs = torch.min(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    min_abs_gpu = np.array([min_abs.cpu().numpy()])
    min_abs_cpu = np.array([np.min(np.abs(signal_cpu))])
    return name, min_abs_gpu, min_abs_cpu, gpu_extraction_time

  def calculate_var(self, signal_gpu, signal_cpu):
    name = 'variance'
    start_time = time.time()
    var_ = torch.var(signal_gpu)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    var_gpu = np.array([var_.cpu().numpy()])
    var_cpu = np.array([np.var(signal_cpu)])
    return name, var_gpu, var_cpu, gpu_extraction_time

  def calculate_var_abs(self, signal_gpu, signal_cpu):
    name = 'variance abs'
    start_time = time.time()
    var_abs = torch.var(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    var_abs_gpu = np.array([var_abs.cpu().numpy()])
    var_abs_cpu = np.array([np.var(np.abs(signal_cpu))])
    return name, var_abs_gpu, var_abs_cpu, gpu_extraction_time

  def calculate_cv(self, signal_gpu, signal_cpu):
    name = 'coefficient variance'
    start_time = time.time()
    mean_gpu = torch.mean(signal_gpu)
    std_dev_gpu = torch.std(signal_gpu)
    cv = std_dev_gpu / mean_gpu
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    cv_gpu = np.array([cv.cpu().numpy()])
    mean_cpu = np.mean(signal_cpu)
    std_dev_cpu = np.std(signal_cpu)
    cv_cpu = np.array([std_dev_cpu / mean_cpu])
    return name, cv_gpu, cv_cpu, gpu_extraction_time

  def calculate_iqr(self, signal_gpu, signal_cpu):
    name = 'inter-quartile range'
    start_time = time.time()
    sorted_data = torch.sort(signal_gpu).values
    q1_index = int(0.25 * (len(sorted_data) - 1))
    q3_index = int(0.75 * (len(sorted_data) - 1))
    q1_value = sorted_data[q1_index]
    q3_value = sorted_data[q3_index]
    iqr_ = q3_value - q1_value
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    iqr_gpu = np.array([iqr_.cpu().numpy()])
    iqr_cpu = np.array([np.percentile(signal_cpu, 75) - np.percentile(signal_cpu, 25)])
    return name, iqr_gpu, iqr_cpu, gpu_extraction_time

  def calculate_root_mean_square(self, signal_gpu, signal_cpu):
    name = 'root mean square'
    start_time = time.time()
    rms_ = torch.sqrt(torch.mean(signal_gpu ** 2))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    rms_gpu = np.array([rms_.cpu().numpy()])
    rms_cpu = np.array([np.sqrt(np.mean(signal_cpu ** 2))])
    return name, rms_gpu, rms_cpu, gpu_extraction_time

  def calculate_energy(self, signal_gpu, signal_cpu):
    name = 'energy'
    start_time = time.time()
    energy_ = torch.sum(signal_gpu ** 2)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    energy_gpu = np.array([energy_.cpu().numpy()])
    energy_cpu = np.array([np.sum(signal_cpu ** 2)])
    return name, energy_gpu, energy_cpu, gpu_extraction_time

  def calculate_log_energy(self, signal_gpu, signal_cpu):
    name = 'log energy'
    start_time = time.time()
    energy_log = torch.log(torch.sum(signal_gpu ** 2))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    energy_log_gpu = np.array([energy_log.cpu().numpy()])
    energy_log_cpu = np.array([np.log(np.sum(signal_cpu ** 2))])
    return name, energy_log_gpu, energy_log_cpu, gpu_extraction_time

  def calculate_entropy(self, signal_gpu, signal_cpu):
    name = 'entropy'
    start_time = time.time()
    values, counts = torch.unique(signal_gpu, return_counts=True)
    probabilities = counts / counts.sum()
    entropy_ = -torch.sum(probabilities * torch.log2(probabilities))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    entropy_gpu = np.array([entropy_.cpu().numpy()])
    values, counts = np.unique(signal_cpu, return_counts=True)
    probabilities = counts / counts.sum()
    entropy_cpu = np.array([-np.sum(probabilities * np.log2(probabilities))])
    return name, entropy_gpu, entropy_cpu, gpu_extraction_time

  def calculate_zero_crossings(self, signal_gpu, signal_cpu):
    name = 'zero crossings'
    start_time = time.time()
    # Compute sign changes: A sign change occurs where the product of adjacent elements is negative
    sign_changes = torch.mul(signal_gpu[:-1], signal_gpu[1:]) < 0
    zero_crossings = torch.sum(sign_changes)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    zero_crossings_gpu = np.array([zero_crossings.cpu().numpy()])
    zero_cross_diff = np.diff(np.signbit(signal_cpu))
    zero_crossings_cpu = np.array([zero_cross_diff.sum()])
    return name, zero_crossings_gpu, zero_crossings_cpu, gpu_extraction_time

  def calculate_crest_factor(self, signal_gpu, signal_cpu):
    name = 'crest factor'
    start_time = time.time()
    crest_factor = torch.max(torch.abs(signal_gpu)) / torch.sqrt(torch.mean(signal_gpu ** 2))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    crest_factor_gpu = np.array([crest_factor.cpu().numpy()])
    crest_factor_cpu = np.array([np.max(np.abs(signal_cpu)) / np.sqrt(np.mean(signal_cpu ** 2))])
    return name, crest_factor_gpu, crest_factor_cpu, gpu_extraction_time

  def calculate_mean_crossings(self, signal_gpu, signal_cpu):
    name = 'mean crossings'
    start_time = time.time()
    data_minus_mean = signal_gpu - torch.mean(signal_gpu)
    mean_crossings = torch.mul(data_minus_mean[:-1], data_minus_mean[1:]) < 0
    mean_crossings_count = torch.sum(mean_crossings)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    mean_crossings_gpu = np.array([mean_crossings_count.cpu().numpy()])
    mean_crossings_cpu = np.array([len(np.where(np.diff(np.sign(signal_cpu - np.mean(signal_cpu))))[0])])
    return name, mean_crossings_gpu, mean_crossings_cpu, gpu_extraction_time

  def calculate_mean_auto_correlation(self, signal_gpu, signal_cpu):
    name = 'mean auto correlation'
    n_lags = self.n_lags_auto_correlation
    #n = len(signal_gpu)
    start_time = time.time()
    autocorrelations = torch.zeros(n_lags).cuda()
    #mean_ = torch.mean(signal_gpu)
    for lag in range(1, n_lags + 1): #exclude the lag=0 here
      autocorrelations[lag - 1] = torch.mean((signal_gpu[:-lag] - torch.mean(signal_gpu[:-lag])) * (signal_gpu[lag:] - torch.mean(signal_gpu[lag:])))
      #autocorrelations[lag - 1] = torch.sum((signal_gpu[lag:] - mean_) * (signal_gpu[:-lag] - mean_)) / (n - lag)
    mean_autocorrelation = torch.mean(autocorrelations)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    mean_autocorrelation_gpu = np.array([mean_autocorrelation.cpu().numpy()])
    mean_autocorrelation_cpu = np.array([np.mean(acf(signal_cpu, nlags=self.n_lags_auto_correlation)[1:])])
    return name, mean_autocorrelation_gpu, mean_autocorrelation_cpu, gpu_extraction_time

  def calculate_median_abs_dev(self, signal_gpu, signal_cpu):
    name = 'median absolute deviation'
    start_time = time.time()
    median_ = torch.median(signal_gpu)
    abs_deviation = torch.abs(signal_gpu - median_)
    mad_ = torch.median(abs_deviation)
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    mad_gpu = np.array([mad_.cpu().numpy()])
    mad_cpu = np.array([np.median(np.abs(signal_cpu - np.median(signal_cpu)))])
    return name, mad_gpu, mad_cpu, gpu_extraction_time

  def calculate_magnitude_area(self, signal_gpu, signal_cpu):
    name = 'magnitude area'
    start_time = time.time()
    magnitude_area = torch.sum(torch.abs(signal_gpu))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    magnitude_area_gpu = np.array([magnitude_area.cpu().numpy()])
    magnitude_area_cpu = np.array([np.sum(np.abs(signal_cpu))])
    return name, magnitude_area_gpu, magnitude_area_cpu, gpu_extraction_time

  def calculate_avg_amplitude_change(self, signal_gpu, signal_cpu):
    name = 'average amplitude change'
    start_time = time.time()
    average_amplitude_change = torch.mean(torch.abs(torch.diff(signal_gpu)))
    end_time = time.time()
    gpu_extraction_time = end_time - start_time
    average_amplitude_change_gpu = np.array([average_amplitude_change.cpu().numpy()])
    average_amplitude_change_cpu = np.array([np.mean(np.abs(np.diff(signal_cpu)))])
    return name, average_amplitude_change_gpu, average_amplitude_change_cpu, gpu_extraction_time





  def calculate_statistial_features(self, signal_gpu, signal_cpu):

    feature_name = []
    feature_value_gpu = []
    feature_value_cpu = []
    feature_extraction_time = []

    # the following code can be write in a loop after finalize all features

    # mean
    name, value_gpu, value_cpu, time = self.calculate_mean(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    #geometric mean
    name, value_gpu, value_cpu, time = self.calculate_geometric_mean(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    #harmonic mean
    name, value_gpu, value_cpu, time = self.calculate_harmonic_mean(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    #trimmed mean
    name, value_gpu, value_cpu, time = self.calculate_trimmed_means(signal_gpu, signal_cpu)
    feature_name.extend(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.extend(time)

    # mean abs
    name, value_gpu, value_cpu, time = self.calculate_mean_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # geometric mean abs
    name, value_gpu, value_cpu, time = self.calculate_geometric_mean_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # harmonic mean abs
    name, value_gpu, value_cpu, time = self.calculate_harmonic_mean_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    #trimmed mean abs
    name, value_gpu, value_cpu, time = self.calculate_trimmed_means_abs(signal_gpu, signal_cpu)
    feature_name.extend(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.extend(time)

    # std
    name, value_gpu, value_cpu, time = self.calculate_std(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # std abs
    name, value_gpu, value_cpu, time = self.calculate_std_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # skewness
    name, value_gpu, value_cpu, time = self.calculate_skewness(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # skewness abs
    name, value_gpu, value_cpu, time = self.calculate_skewness_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # kurtosis
    name, value_gpu, value_cpu, time = self.calculate_kurtosis(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # kurtosis abs
    name, value_gpu, value_cpu, time = self.calculate_kurtosis_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # median
    name, value_gpu, value_cpu, time = self.calculate_median(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # median abs
    name, value_gpu, value_cpu, time = self.calculate_median_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # max
    name, value_gpu, value_cpu, time = self.calculate_max(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # max abs
    name, value_gpu, value_cpu, time = self.calculate_max_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # min
    name, value_gpu, value_cpu, time = self.calculate_min(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # min abs
    name, value_gpu, value_cpu, time = self.calculate_min_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # var
    name, value_gpu, value_cpu, time = self.calculate_var(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # var abs
    name, value_gpu, value_cpu, time = self.calculate_var_abs(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # coefficient of variation
    name, value_gpu, value_cpu, time = self.calculate_cv(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # inter-quartile range
    name, value_gpu, value_cpu, time = self.calculate_iqr(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # root mean square
    name, value_gpu, value_cpu, time = self.calculate_root_mean_square(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # energy
    name, value_gpu, value_cpu, time = self.calculate_energy(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # log energy
    name, value_gpu, value_cpu, time = self.calculate_log_energy(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # entropy
    name, value_gpu, value_cpu, time = self.calculate_entropy(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # zero crossings
    name, value_gpu, value_cpu, time = self.calculate_zero_crossings(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # crest factor
    name, value_gpu, value_cpu, time = self.calculate_crest_factor(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # mean crossings
    name, value_gpu, value_cpu, time = self.calculate_mean_crossings(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # mean auto correlation
    name, value_gpu, value_cpu, time = self.calculate_mean_auto_correlation(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # median absolute deviation
    name, value_gpu, value_cpu, time = self.calculate_median_abs_dev(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # magnitude area
    name, value_gpu, value_cpu, time = self.calculate_magnitude_area(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)

    # average amplitude area
    name, value_gpu, value_cpu, time = self.calculate_avg_amplitude_change(signal_gpu, signal_cpu)
    feature_name.append(name)
    feature_value_gpu.extend(value_gpu)
    feature_value_cpu.extend(value_cpu)
    feature_extraction_time.append(time)



    feature_information = pd.DataFrame({
          'feature_name': feature_name,
          'feature_value_gpu': feature_value_gpu,
          'feature_value_cpu': feature_value_cpu,
          'feature_extraction_time': feature_extraction_time
    })
    return feature_information

# using np to generate random signal data with 1000 data instances
np.random.seed(0)
signal_data = np.random.rand(1000)
signal_data = 10 * signal_data  #[0,10)

# Check if GPU (CUDA) is available
if torch.cuda.is_available():
    device = torch.device("cuda")
    print("Using GPU for computations.")
else:
    # Warning if GPU is not available
    print("Warning: GPU not available, please ensure that a GPU is available for this code to run properly.")
    # Raise an exception to stop the execution
    raise RuntimeError("GPU not available")

# in order to compare the feature compute from cpu and gpu, one use orignal data, one use tensor
data_cpu = signal_data

# Convert signal data as a PyTorch tensor
data_gpu = torch.tensor(signal_data, dtype=torch.float32).to(device)
statistical_feature = StatisticalFeatures(window_size = 2, trimmed_mean_thresholds=None)
statistical_feature.calculate_statistial_features(data_gpu, data_cpu)