Power analysis

Demonstration code for performing a power analysis using PowerAnalysis().

Link to script: feature_extraction/power_analysis.py

import os
import matplotlib.pyplot as plt
import numpy as np

from nnsa import PowerAnalysis, EdfReader, print_object_summary, SUPPORTED_RESULT_FILE_TYPES, read_result_from_file, \
    assert_equal

plt.close('all')

Parameters.

# Print the default parameters of PowerAnalysis():
print(PowerAnalysis().default_parameters())

# Descriptions of the parameters are documented in the default_parameters() code.

# Create an instance of the PowerAnalysis class with custom parameters, overruling some defaults:
power_analysis = PowerAnalysis(artefact_criteria={'max_diff': 200},
                               psd={'window_length': 2.5},
                               segmentation={'segment_length': 15, 'overlap': 1})

# See if the custom parameters were accepted:
print('\nCustom parameters:')
print(power_analysis.parameters)

Main method: power_analysis.

# Now that we have initialized a PowerAnalysis object with certain parameters, we can run the power_analysis method,
# which performs the actual power analysis:
# Create 10 random signals to simulate 10 channel EEG data (amplitude 50).
fs = 250
data_matrix = np.random.rand(10, 100000)*100 - 50
result = power_analysis.power_analysis(data_matrix, fs=fs, unit='uV', verbose=1)  # Not that some args are optional.

# The returned object is a PowerAnalysisResult object, which is a high-level interface for
# manipulating/visualizing/saving the result:
print('\nSummary of result object:')
print_object_summary(result)

PowerAnalysisResult attributes.

# The main attribute is 'power', which contains the frequency power in certain bands:
print('result.power: {}'.format(result.power))

# power is a dictionary with the frequency bands as keys and its values are 2D arrays corresponding to
# (channels, segments).
freq_bands = list(result.power.keys())
print('Frequency bands: {}'.format(freq_bands))
print('power[\'{}\'].shape: {}'.format(freq_bands[0], result.power[freq_bands[0]].shape))

# The time (in seconds) corresponding to the segments in the data in power:
print('segment_times: {}'.format(result.segment_times))

PowerAnalysisResult methods.

# Print a summary of the power analysis:
result.summary()

# Plot the power for all channels as function of segments.
result.plot()

Save the result.

# We can save the result to any supported file. To list the supported file types/extensions:
print('Supported result file types: {}'.format(SUPPORTED_RESULT_FILE_TYPES))

# E.g. save as csv (the file type to save is automatically detected from the extension of the filename).
filename = 'temp_result.csv'
result.save_to_file(filename)

# We can reload the result back to a PowerAnalysisResult object:
result_loaded = read_result_from_file(filename)

# Remove temporary file.
os.remove(filename)

# Verify that the loaded object is equal to the original object:
assert_equal(actual=result_loaded, desired=result)  # Will raise an AssertionError if not equal.
print('Loaded result object is equal to original object.')

PowerAnalysis of EegDataset.

# Specify a file and open a reader to read the data (e.g. EdfReader to read an EDF(+) file).
filepath = 'C:/data_temp/test.edf'

# Read EEG channels in an EegDataset object.
with EdfReader(filepath) as r:
    ds = r.read_eeg_dataset(dtype=np.float32)  # We can specify the dtype, lower precision, means less memory usage.

# We can do the power analysis for multichannel EEG via the wrapper method 'power_analysis' of the EegDataset class.
# Specify the parameters to the PowerAnalysis class as kwargs. Preprocessing and power analysis can be done in one line:
result_ds = ds.reference('Cz').resample(fs_new=128).filter_saved_filter('bandpassfir_a').power_analysis(segmentation={'overlap': 0})

# A data info string is appended to power_result object, which summarizes the source and preprocessing of the data.
print('result_ds.data_info: {}'.format(result_ds.data_info))

# Note that the dtype of the result is the same as the dtype of the data (e.g. np.float32).
print('dtype result: {}'.format(result_ds.power[list(result_ds.power.keys())[0]].dtype))

# Print a summary.
result_ds.summary()

# Plot the power for all channels as function of segments.
plt.figure()
result_ds.plot()