Psd
===

Demonstration code for computing power spectral density using Psd().

Link to script: `feature_extraction/psd.py <https://gitlab.com/timhermans/nnsa_public/-/blob/main/examples/feature_extraction/psd.py>`_



.. code-block:: python

    import os
    import numpy as np
    import matplotlib.pyplot as plt
    
    from nnsa import Psd, print_object_summary, SUPPORTED_RESULT_FILE_TYPES, read_result_from_file, assert_equal, \
        TimeSeries
    
    plt.close('all')


Parameters. 



.. code-block:: python

    # Print the default parameters of Psd():
    print(Psd().default_parameters())
    
    # Descriptions of the parameters are documented in the default_parameters() code.
    
    # Create an instance of the Psd class with custom parameters, overruling some defaults:
    psd = Psd(window_length=2,
              df_max=0.1,
              detrend='constant')
    
    # See if the custom parameters were accepted:
    print('\nCustom parameters:')
    print(psd.parameters)


Main method: psd. 



.. code-block:: python

    # Now that we have initialized a Psd object with certain parameters, we can run the psd method, which computes the psd:
    # Create a random signal.
    signal = np.random.rand(10000)
    fs = 250
    result = psd.psd(signal, fs=fs, unit='uV', axis=-1, verbose=0)  # Note that unit, axis and verbose are optional args.
    
    # The returned object is a PsdResult object, which is a high-level interface for manipulating/visualizing/saving the
    # result:
    print('\nSummary of result object:')
    print_object_summary(result)


PsdResult attributes. 



.. code-block:: python

    # We could check whether the df does not exceed the df_max parameter that we specified:
    print('df_max: {}\ndf: {}'.format(psd.parameters['df_max'], result.df))
    
    # The main attribute is power_density, which contains the actual psd:
    print('power_density: {}'.format(result.power_density))
    
    # The frequencies corresponding to power_density can be constructed from df. This is done by the frequencies property:
    print('frequencies: {}'.format(result.frequencies))


PsdResult  methods. 



.. code-block:: python

    # Plot the power spectral density.
    result.plot(label='result')
    
    # Compute the power in a specific frequency band, with band limits f_low and f_high:
    f_low = 0
    f_high = result.frequencies[round(len(result.power_density)/2)]
    power = result.compute_power(f_low, f_high)
    print('Power {:.2f}-{:.2f} Hz = {:.2f} {}'.format(f_low, f_high, power, result.unit.replace('/Hz', '')))


Save the result. 



.. code-block:: python

    # 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 hdf5 (the file type to save is automatically detected from the extension of the filename).
    filename = 'temp_result.hdf5'
    result.save_to_file(filename)
    
    # We can reload the result back to a PsdResult 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.')


Psd of TimeSeries. 



.. code-block:: python

    # We can also directly call psd on a TimeSeries object, where we specify the parameters to the psd function as kwargs:
    ts = TimeSeries(signal=signal, fs=fs, label='Cz', unit='uV', info={'source': 'random'})
    result_ts = ts.psd(window_length=2,
                       df_max=0.1,
                       detrend='constant')
    
    # If we chose the same parameters, result and result_ts should be equal, except for the data_info object of the result,
    # which contains the info dictionary of the TimeSeries:
    print('result:\n', result)
    print('result_ts:\n', result_ts)
    
    result_ts.plot('--', label='result_ts')
    plt.legend()