Source code for mojo.analysis.acousticPressureScales

"""
acousticPressureScales.py
=========================

Module with methods to process 1D data sample with extension especially for acoustic applications.
"""

import math
import numpy as np

from .windowedFFT import WindowedFFT

# @todo: implement output in 1/3 octave bands


def get_a_filter(frequencies):
    """
    returns the A-weighting for given frequencies
    """
    freq_square = np.square(frequencies)  # pylint: disable=assignment-from-no-return
    a_filter = math.pow(12000, 2) * np.power(frequencies, 4) / ((freq_square + math.pow(20.6, 2)) *
                                                                (freq_square + math.pow(12200, 2)) * np.sqrt(freq_square + math.pow(107.7, 2)) *
                                                                np.sqrt(freq_square + math.pow(737.9, 2)))

    a_filter = np.ma.masked_less_equal(a_filter, 0., copy=False)
    a_filter = 20 * np.ma.log10(a_filter) + 2.

    return a_filter.filled(fill_value=0.)


def get_b_filter(frequencies):
    """
    returns the B-weighting for given frequencies
    """

    freq_square = np.square(frequencies)  # pylint: disable=assignment-from-no-return
    b_filter = math.pow(12000, 2) * np.power(frequencies, 3) / ((freq_square + math.pow(20.6, 2)) * (freq_square +
                                                                                                     math.pow(12200, 2)) * np.sqrt(freq_square + math.pow(158.5, 2)))

    b_filter = np.ma.masked_less_equal(b_filter, 0., copy=False)
    b_filter = 20 * np.ma.log10(b_filter) + 0.17

    return b_filter.filled(fill_value=0.)


def get_c_filter(frequencies):
    """
    returns the C-weighting for given frequencies
    """
    freq_square = np.square(frequencies)  # pylint: disable=assignment-from-no-return
    c_filter = math.pow(12000, 2) * freq_square / ((freq_square + math.pow(20.6, 2)) *
                                                   (freq_square + math.pow(12200, 2)))

    c_filter = np.ma.masked_less_equal(c_filter, 0., copy=False)
    c_filter = 20 * np.ma.log10(c_filter) + 0.06

    return c_filter.filled(fill_value=0.)


class AcousticPressureScales:
    """
    Class to process 1D data sample with extension especially for acoustic applications.
    """

    def __init__(self, pressure_time, d_time, pressure_reference=0.00002):

        self.__frequency = None

        self.__pressure_time = pressure_time
        self.__pressure_frq = None
        self.__pressure_frq_log = None

        self.__pref = pressure_reference

        self.__d_time = d_time
        self.__window_type = ""
        self.__window_size = 0
        self.__window_overlap = 0

        self.set_windowing_params("FLAT", self.__pressure_time.size, self.__window_overlap)

    def set_windowing_params(self, window_type, window_size, window_overlap):
        """
        set or change the parameter for the FFT
        """
        self.__window_type = window_type
        self.__window_size = window_size
        self.__window_overlap = window_overlap

    def get_autocorrelation(self):
        """
        calculate the auto correlation of a signal
        """
        signal_length = self.__pressure_time.size
        variance = np.var(self.__pressure_time)
        auto_corr = np.correlate(self.__pressure_time, self.__pressure_time, mode='full')[-signal_length:]
        auto_corr_norm = auto_corr / (variance * (np.arange(signal_length, 0, -1)))

        return auto_corr_norm

    def __perform_fft(self, pressure):
        """
        perform and store the FFT
        """
        if self.__pressure_frq is None:
            data_fft = WindowedFFT(pressure, self.__d_time, self.__window_type, self.__window_size)
            self.__frequency = data_fft.getSingleSidedFrequencies()
            self.__pressure_frq = data_fft.getSingleSidedSpectrum()

    def __compute_log_scale(self):
        """
        convert a equidistant linear frequency array into a equidistant log-scale array
        """
        if self.__pressure_frq_log is None:
            self.__pressure_frq_log = 10 * np.log10(np.square(np.abs(self.__pressure_frq) / self.__pref))

    def get_spl(self, filter_name=None):
        """
        Convert the pressure amplitudes to SPL based unit
        """
        self.__perform_fft(self.__pressure_time)
        self.__compute_log_scale()

        if filter_name is None:
            spl = self.__pressure_frq_log
        else:
            if filter_name.upper() == "A":
                filter_func = get_a_filter(self.__frequency)
            elif filter_name.upper() == "B":
                filter_func = get_b_filter(self.__frequency)
            elif filter_name.upper() == "C":
                filter_func = get_c_filter(self.__frequency)
            else:
                raise ValueError("No filter function defined with name '%s'." % filter_name)

            spl = filter_func * self.__pressure_frq_log

        return self.__frequency, spl

    def get_psd(self):
        """
        Convert the pressure amplitudes to PSD based unit
        """
        auto_corr = self.get_autocorrelation()
        self.__perform_fft(auto_corr)
        self.__compute_log_scale()

        return self.__frequency, self.__pressure_frq_log