annotate pyphon.py @ 11:ca98f5f26bcb tip

* More frequencies to verify toolbox with.

"""
A module used for auditory analysis.

Models currently implemented:
    * Simple frequency-modulation analysis model.
    * Glasberg and Moore model using gammatone filters
    
Packaged dependencies:
    > utils.py and/or utils.pyc
    > erb.dat
    > outMidFir.dat
    > tq.dat
    
External dependencies:
    > scipy
    > numpy
    > matplotlib
    > sys (Built-in)
    
Primary functions:
    > calc_loudness_quiet
    > calc_loudness_noise
    > get_excitation_i
    > plot_excitation_response
"""

import sys
import utils
import scipy.signal as sp
import scipy.fftpack as fft
import numpy as np
import matplotlib.pyplot as plt

def calc_loudness_noise(input, fs, SPL = None, inputType = 'signal', gammatoneLength=4096, rectify=False, inc_loud_region = True, verbose = False):
    """"
        Calculate the specifc loudness, instantaneous loudness, short term loudness and long term loudness of the input signal present in noise.
        
        Parameters:
        * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1 or excitation pattern. (Required)
        * fs (type: numerical) - sample frequency of the signal, input. (Required)
        * SPL (type: float) - the sound pressure level corresponding to 1 of the normalised range -1 to 1. (Optional if inputType='excitation' then doesn't do anything. Required otherwise)
        * inputType (type: string) - specify whether the input is the excitation pattern, or the input signal. Use 'signal' or 'excitation', respectively. (Optional; Default = 'signal')
        * gammaToneLength (type: int) - length of gammatone filters. (Optional; Default = 4096)
        * rectify (type = boolean) - rectify the input to model the direction shearing of the hairs in cochlea. Glasberg and Moore model doesn't consider this, so by default is false. (Optional; Default = False)
        * inc_loud_region (type: boolean) - Specifies whether to process the high levels differently to the mid levels. (Optional; Default = False)
        * Verbose (type: boolean) - print verbosely to stdout.
        
        Returns:
        * sl (type: numpy array of floats) - specific loudness
        * loudness (type: numpy array of floats) - instantaneous loudness
        * STL (type: numpy array of floats) - short term loudness
        * LTL (type: numpy array of floats) - long term loudness
        """
    
    
    if(inputType == 'signal' or inputType == 'sig'):
        if(SPL == None):
            raise ValueError("You need to provide the full scale sound pressure level!")
        
        numTracks = np.shape(input)[0]
        lenSig = np.shape(input)[-1]
        excitation = np.zeros((numTracks, 39, lenSig+gammatoneLength+4095))
    
        for i in range(numTracks):
            excitation[i] = get_excitation_i(input[i], fs, SPL, gammatoneLength=gammatoneLength, rectify = rectify, verbose = verbose)

    elif(inputType == 'excitation' or inputType == 'exc'):
        excitation = input
        numTracks = np.shape(excitation)[0]
    else:
        raise ValueError("Unknown inputType. Can be 'signal' (or 'sig') or 'excitation' (or 'exc')")

    del input
    
    sl = _get_specific_loudness_noise(excitation, inc_loud_region = inc_loud_region)
    loudness = 2*np.sum(sl,1)
    STL = np.zeros(np.shape(loudness))
    LTL = np.array(STL)
    
    for i in range(numTracks):
        STL[i] = temporal_integration(loudness[i], fs, 'STL')
        LTL[i] = temporal_integration(STL[i], fs, 'LTL')

    return sl, loudness, STL, LTL

def calc_loudness_quiet(input, fs, SPL = None, inputType = 'signal', gammatoneLength=4096, rectify = False, inc_loud_region = True, verbose = False):
    """"
    Calculate the specifc loudness, instantaneous loudness, short term loudness and long term loudness of the input signal.
        
        Parameters:
        * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1 or excitation pattern. (Required)
        * fs (type: numerical) - sample frequency of the signal, input. (Required)
        * SPL (type: float) - the sound pressure level corresponding to 1 of the normalised range -1 to 1. (Optional if inputType='excitation' then doesn't do anything. Required otherwise)
        * inputType (type: string) - specify whether the input is the excitation pattern, or the input signal. Use 'signal' or 'excitation', respectively. (Optional; Default = 'signal')
        * gammaToneLength (type: int) - length of gammatone filters. (Optional; Default = 4096)
        * rectify (type = boolean) - rectify the input to model the direction shearing of the hairs in cochlea. Glasberg and Moore model doesn't consider this, so by default is false. (Optional; Default = False)
        * inc_loud_region (type: boolean) - Specifies whether to process the high levels differently to the mid levels. (Optional; Default = False)
        * Verbose (type: boolean) - print verbosely to stdout.
        
        Returns:
        * sl (type: numpy array of floats) - specific loudness
        * loudness (type: numpy array of floats) - instantaneous loudness
        * STL (type: numpy array of floats) - short term loudness
        * LTL (type: numpy array of floats) - long term loudness
        """
    
    
    if(inputType == 'signal' or inputType == 'sig'):
        if(SPL == None):
            raise ValueError("You need to provide the full scale sound pressure level!")
        excitation = get_excitation_i(input, fs, SPL, gammatoneLength=gammatoneLength, rectify = rectify, verbose = verbose)
    elif(inputType == 'excitation' or inputType == 'exc'):
        excitation = input
    else:
        raise ValueError("Unknown inputType. Can be 'signal' (or 'sig') or 'excitation' (or 'exc')")

    del input

    sl = _get_specific_loudness_quiet(excitation, inc_loud_region = inc_loud_region)
    loudness = 2*np.sum(sl,0)
    STL = temporal_integration(loudness, fs, 'STL')
    LTL = temporal_integration(STL, fs, 'LTL')
    
    return sl, loudness, STL, LTL

def _get_specific_loudness_noise(exc_i, analyse = "all", inc_loud_region = True):
    """
        A function to calculate the specific loudness of the excitation patterns at each ERB when present in noise. Specific loudness is calculated for three regions of signal intensity, low level, mid level and high level. The high level processing is optional, and by default, ignored. When ignored the high level region is processed equivalent to that of the mid region.
        
        Parameters:
        * exc_i (type: array-like of floats) - an array of shape (39, lenSig) of the excitation intensity of the signal. The outer dimension determines the ERB channel. (Required)
        * inc_loud_region (type: boolean) - Specifies whether to process the high levels differently to the mid levels. (Optional; Default = False)
        """
    exc_i = np.array(exc_i)

    lenSig = np.shape(exc_i)[-1] #length signal
    numTracks = np.shape(exc_i)[0]
    tq, G, A, alpha = _get_specific_loudness_parameters(lenSig)
    K = _calc_k_parameter(lenSig)
    specific_loudness = exc_low = exc_mid = exc_high = np.zeros(np.shape(exc_i))

    exc_signoise = sum(exc_i, 0)
    signoise_gt_tl = _get_greater_than_tl(exc_signoise) #boolean array of elements greater than the loud threshold
    signoise_lt_tl = ~signoise_gt_tl
    
    if analyse.lower() == "all":
        loop = range(numTracks);
    elif analyse.lower() == "first":
        loop = range(1)
    
    for i in loop:
        
        exc_sig = exc_i[i]
        exc_noise = exc_signoise - exc_sig
        tn = K*exc_noise + tq

        less_than_tn = _get_less_than_tq(exc_sig, tn) #boolean array of elements less than or equal to the threshold of quiet
        greater_than_tn = ~less_than_tn #boolean array of elements greater than the threshold of quiet
        
        if inc_loud_region:
            c1 = signoise_lt_tl * greater_than_tn #case 1
            c2 = signoise_lt_tl * less_than_tn #case 2
            c3 = signoise_gt_tl * greater_than_tn #case 3
            c4 = signoise_gt_tl * less_than_tn #case 4
            specific_loudness[i, c3] = _get_specific_loudness_noise_c3(exc_sig[c3], exc_noise[c3], tq[c3], tn[c3], G[c3], A[c3], alpha[c3])
            specific_loudness[i, c4] = _get_specific_loudness_noise_c4(exc_sig[c4], exc_noise[c4], tq[c4], tn[c4], G[c4], A[c4], alpha[c4])
        else:
            c1 = greater_than_tn #case 1
            c2 = less_than_tn #case 2

        specific_loudness[i][c1] = _get_specific_loudness_noise_c1(exc_sig[c1], exc_noise[c1], tq[c1], tn[c1], G[c1], A[c1], alpha[c1])
        specific_loudness[i][c2] = _get_specific_loudness_noise_c2(exc_sig[c2], exc_noise[c2], tq[c2], tn[c2], G[c2], A[c2], alpha[c2])

    return specific_loudness

def _get_specific_loudness_noise_c1(exc_sig, exc_noise, tq, tn, G, A, alpha):
    """
    A function to calculate the specific loudness of a component of the signal for the case exc_noise < 10**10, and exc_quiet > threshold noise (c1).
    
    Parameters:
    * tq (type: numpy array of floats) - the threshold of quiet for each centre frequency of the ERBs
    * tn (type: numpy array of floats) - the threshold of quiet when present in noise for each centre frequency of the ERBs
    * G (type: numpy array of floats) - the frequency dependent gain values for each centre frequency of the ERBs
    * A (type: numpy array of floats) - the frequency dependent A values for each centre frequency of the ERBs
    * alpha (type: numpy array of floats) - the frequency dependent alpha values for each centre frequency of the ERBs
    
    Returns:
        * specific_loudness_c1 (type: array-like matrix of floats) - An array with dimensions equal to the exc_sig containing the
    specific loudness c1.
    """
    C = 0.046871 #constant parameter
    sl_c1 = C*(((exc_sig + exc_noise)*G + A)**alpha - A**alpha) - C*((((tn + exc_noise)*G + A)**alpha) - (tq*G + A)**alpha) * (tn/exc_sig)**0.3

    return sl_c1

def _get_specific_loudness_noise_c2(exc_sig, exc_noise, tq, tn, G, A, alpha):
    """
    A function to calculate the specific loudness of a component of the signal for the case exc_noise < 10**10, and exc_quiet < threshold noise (c2).
    
    Parameters:
        * tq (type: numpy array of floats) - the threshold of quiet for each centre frequency of the ERBs
        * tn (type: numpy array of floats) - the threshold of quiet when present in noise for each centre frequency of the ERBs
        * G (type: numpy array of floats) - the frequency dependent gain values for each centre frequency of the ERBs
        * A (type: numpy array of floats) - the frequency dependent A values for each centre frequency of the ERBs
    * alpha (type: numpy array of floats) - the frequency dependent alpha values for each centre frequency of the ERBs
    
    Returns:
        * specific_loudness_c2 (type: array-like matrix of floats) - An array with dimensions equal to the exc_sig containing the
    specific loudness c2.
    """
    
    C = 0.046871 #constant parameter
    sl_c2 = (C*((2*exc_sig)/(exc_sig + tn))**1.5) * (((tq*G + A)**alpha - A**alpha)/(((0.0 + tn + exc_noise)*G + A)**alpha - (exc_noise*G + A)**alpha)) * (((exc_sig + exc_noise)*G + A)**alpha - (exc_noise*G + A)**alpha)
    
    return sl_c2

def _get_specific_loudness_noise_c3(exc_sig, exc_noise, tq, tn, G, A, alpha):
    """
    A function to calculate the specific loudness of a component of the signal for the case exc_noise > 10**10, and exc_quiet > threshold noise (c3).
    
    Parameters:
    * tq (type: numpy array of floats) - the threshold of quiet for each centre frequency of the ERBs
    * tn (type: numpy array of floats) - the threshold of quiet when present in noise for each centre frequency of the ERBs
    * G (type: numpy array of floats) - the frequency dependent gain values for each centre frequency of the ERBs
    * A (type: numpy array of floats) - the frequency dependent A values for each centre frequency of the ERBs
    * alpha (type: numpy array of floats) - the frequency dependent alpha values for each centre frequency of the ERBs
    
    Returns:
        * specific_loudness_c3 (type: array-like matrix of floats) - An array with dimensions equal to the exc_sig containing the
    specific loudness c3.
    """
    C = 0.046871 #constant parameter
    C2 = C/((1.04 * 10**6)**0.5)
    sl_c3 = (C2*(exc_sig + exc_noise)**0.5) - C2*((tn + exc_noise)**0.5 - (tq*G + A)**alpha + A**alpha) * (tn/exc_sig)**0.3
    
    return sl_c3

def _get_specific_loudness_noise_c4(exc_sig, exc_noise, tq, tn, G, A, alpha):
    """
    A function to calculate the specific loudness of a component of the signal for the case exc_noise > 10**10, and exc_quiet < threshold noise (c4).
    
    Parameters:
        * tq (type: numpy array of floats) - the threshold of quiet for each centre frequency of the ERBs
        * tn (type: numpy array of floats) - the threshold of quiet when present in noise for each centre frequency of the ERBs
        * G (type: numpy array of floats) - the frequency dependent gain values for each centre frequency of the ERBs
        * A (type: numpy array of floats) - the frequency dependent A values for each centre frequency of the ERBs
        * alpha (type: numpy array of floats) - the frequency dependent alpha values for each centre frequency of the ERBs
        
    Returns:
        * specific_loudness_c4 (type: array-like matrix of floats) - An array with dimensions equal to the exc_sig containing the
    specific loudness c4.
    
    """
    C = 0.046871 #constant parameter
    sl_c4 = (C*((2*exc_sig)/(exc_sig + tn))**1.5) * (((tq*G + A)**alpha - A**alpha)/((tn + exc_noise)**0.5 - (exc_noise)**0.5)) * ((exc_sig + exc_noise)**0.5 - (exc_noise)**0.5)
    
    return sl_c4

def _get_specific_loudness_quiet(exc_i, inc_loud_region = True):
    """
    A function to calculate the specific loudness of the excitation patterns at each ERB. Specific loudness is calculated for three regions of signal intensity, low level, mid level and high level. The high level processing is optional, and by default, ignored. When ignored the high level region is processed equivalent to that of the mid region.
    
    Parameters:
        * exc_i (type: array-like of floats) - an array of shape (39, lenSig) of the excitation intensity of the signal. The outer dimension determines the ERB channel. (Required)
        * inc_loud_region (type: boolean) - Specifies whether to process the high levels differently to the mid levels. (Optional; Default = False)
        
    Returns:
        * specific_loudness (type: array-like matrix of floats) - specific loudness
    """
    exc_i = np.array(exc_i)
    
    lenSig = np.shape(exc_i)[-1] #length signal
    tq, G, A, alpha = _get_specific_loudness_parameters(lenSig)
    specific_loudness = exc_low = exc_mid = exc_high = np.zeros(np.shape(exc_i))
    
    less_than_tq = _get_less_than_tq(exc_i, tq) #boolean array of elements less than or equal to the threshold of quiet
    greater_than_tq = ~less_than_tq #boolean array of elements greater than the threshold of quiet
    greater_than_tl = _get_greater_than_tl(exc_i) #boolean array of elements greater than the loud threshold
    gttq_lttl = greater_than_tq * ~greater_than_tl #boolean array of elements greater than the threshold of quiet but less than the threshold of loud
    
    exc_low = exc_i[less_than_tq]
    specific_loudness[less_than_tq] = _get_specific_loudness_low(exc_low, G[less_than_tq], tq[less_than_tq], A[less_than_tq], alpha[less_than_tq])
    
    if(inc_loud_region):
        exc_mid = exc_i[gttq_lttl]
        exc_high = exc_i[greater_than_tl]
        specific_loudness[gttq_lttl] = _get_specific_loudness_mid(exc_mid, G[gttq_lttl], A[gttq_lttl], alpha[gttq_lttl])
        specific_loudness[greater_than_tl] = _get_specific_loudness_high(exc_high)
    else:
        exc_mid = exc_i[greater_than_tq]
        specific_loudness[greater_than_tq] = _get_specific_loudness_mid(exc_mid, G[greater_than_tq], A[greater_than_tq], alpha[greater_than_tq])
    
    return specific_loudness

def _get_specific_loudness_parameters(lenSig = 1):
    """
    Loads and returns the specific loudness values. If lenSig is specified, the parameters are shaped equal to the excitation signal to allow for matrix elementwise operations between the parameters and signal. (Assumes excitation signal shape is [39, lenSig]).
    
    Parameters:
        * lenSig (type: numerical int) - the length of the excitation signal to be analysed. (Optional; Default = 1)
        
    Returns:
        * tq (type: numpy array of floats) - the threshold of quiet for each centre frequency of the ERBs
        * G (type: numpy array of floats) - the frequency dependent gain values for each centre frequency of the ERBs
        * A (type: numpy array of floats) - the frequency dependent A values for each centre frequency of the ERBs
        * alpha (type: numpy array of floats) - the frequency dependent alpha values for each centre frequency of the ERBs
    """

    tq_dB, A, alpha = utils.load_sl_parameters() #load tq, A and alpha parameters
    tq_dB = np.transpose(np.tile(tq_dB, (lenSig,1)))
    tq = 10**(tq_dB / 10)
    A = np.transpose(np.tile(A, (lenSig,1)))
    alpha = np.transpose(np.tile(alpha, (lenSig,1)))
    tq500_dB = 3.73 #threshold of quiet at 500Hz reference
    G = 10**((tq500_dB - tq_dB)/10) #gain parameters

    return tq, G, A, alpha

def _calc_k_parameter(lenSig = 1):
    """
    A function to return the K parameter of the loudness compression curve.
    
    Parameters:
        * lenSig (type: int): use if want to shape to the input data for matrix multiplication. (Optional; Default = 1)
    
    Returns:
        * K (type: numpy array of floats) - a parameter to the loudness compression curve
    """
    nerbs = 39
    erbs = np.array(range(1,nerbs+1))
    fc = (10**(erbs/21.366)-1)/0.004368
    
    fkRef = 10**np.append(np.append(np.array(1), np.linspace(np.log10(50), np.log10(500), 11)), np.linspace(3,4.5,16))
    kRef = np.append(np.append(13, np.linspace(13,-3,11)), np.ones((1,16))*-3)
    K = np.interp(fc,fkRef,kRef)
    K = 10**(K/10)
    K = np.transpose(np.tile(K, (lenSig,1)))

    return K

def _get_specific_loudness_low(exc_low, G, tq, A, alpha):
    """
    Returns the specific loudness of the low level parts of the signal. Use _get_specific_loudness_parameters() and specify lenSig
    to obtain the correct values for G, tq, A and alpha. Use _get_less_than_tq() to find the indexes of the samples with lower level
    excitations.
    
    EXAMPLE USE
    
    # exc_i is the excitation intensity
    lenSig = np.shape(exc_i)[-1] # find the length of the signal
    tq, G, A, alpha = _get_specific_loudness_parameters(lenSig) # get the shaped loudness parameters
    less_than_tq = _get_less_than_tq(exc_i, tq) # find which samples are lower level
    specific_loudness[less_than_tq] = _get_specific_loudness_low(exc_low[less_than_tq], G[less_than_tq], tq[less_than_tq], A[less_than_tq], alpha[less_than_tq])
    # only process the low level part of the signal
    
    Parameters:
        * exc_low (type: array-like matrix of floats) - the lower level (less than or equal to tq) excitation pattern for each ERB
            (Required)
        * G (type: array-like matrix of floats) - the frequency dependent loudness gain parameters for each ERB, must be same shape
            as exc_low. (Required)
        * tq (type: array-like matrix of floats) - the frequency dependent threshold of quiet parameters for each ERB. Must be same
            shape as exc_low. (Required)
        * A (type: array-like matrix of floats) - the frequency dependent A parameters for each ERB. Must be same shape as exc_low.
        * alpha (type: array-like matrix of floats) - the frequency dependent alpha parameters for each ERB. Must be same shape as
            exc_low. (Required)
            
    Returns:
        * specific_loudness_low (type: array-like matrix of floats) - An array with dimensions equal to the exc_low containing the
            specific loudness of the signal at levels below the threshold of quiet.
    """
    
    C = 0.046871 #constant parameter
    specific_loudness_low = C * ((2.0*exc_low)/(0.0+exc_low+tq))**1.5 * ((G * exc_low + A)**alpha - A**alpha)

    return specific_loudness_low

def _get_specific_loudness_mid(exc_mid, G, A, alpha):
    """
    Returns the specific loudness of the mid level parts of the signal.
        
    EXAMPLE USE
        
    # exc_i is the excitation intensity
    lenSig = np.shape(exc_i)[-1] # find the length of the signal
    tq, G, A, alpha = _get_specific_loudness_parameters(lenSig) # get the shaped loudness parameters
    less_than_tq = _get_less_than_tq(exc_i, tq) # find which samples are lower level
    greater_than_tq = ~less_than_tq # find which samples are greater than the threshold of quiet
    greater_than_tl = _get_greater_than_tl(exc_i) # find which samples are greater than the loud threshold
    gttq_lttl = greater_than_tq[~greater_than_tl] # find which samples are mid level
    specific_loudness[gttq_lttl] = _get_specific_loudness_low(exc_low[gttq_lttl], G[gttq_lttl], tq[gttq_lttl], A[gttq_lttl], alpha[gttq_lttl])
    # only process the mid level part of the signal
    
    NOTE: The above is an example of use assuming the higher level processing IS NOT IGNORED. Use variable greater_than_tq if processing higher levels equivalent to the mid.
    
    
    Parameters:
        * exc_mid (type: array-like matrix of floats) - the mid level (larger than tq (and, optionally less than high level threshold))
            excitation pattern for each ERB (Required)
        * G (type: array-like matrix of floats) - the frequency dependent loudness gain parameters for each ERB, must be same shape as exc_low. (Required)
        * A (type: array-like matrix of floats) - the frequency dependent A parameters for each ERB. Must be same shape as exc_low.
        * alpha (type: array-like matrix of floats) - the frequency dependent alpha parameters for each ERB. Must be same shape as exc_low. (Required)
        
    Returns:
        * specific_loudness_mid (type: array-like matrix of floats) - An array with dimensions equal to the exc_mid containing the specific loudness of the signal at mid levels.
    """
    
    C = 0.046871 #constant parameter
    specific_loudness_mid = C * ((G * exc_mid + A)**alpha - (A**alpha))
    
    return specific_loudness_mid

def _get_specific_loudness_high(exc_high):
    """
    Returns the specific loudness of the high level parts of the signal.
        
    EXAMPLE USE
        
    # exc_i is the excitation intensity
    lenSig = np.shape(exc_i)[-1] # find the length of the signal
    tq, G, A, alpha = _get_specific_loudness_parameters(lenSig) # get the shaped loudness parameters
    greater_than_tl = _get_greater_than_tl(exc_i) # find which samples are greater than the loud threshold
    specific_loudness[greater_than_tl] = _get_specific_loudness_low(exc_low[greater_than_tl], G[greater_than_tl], tq[greater_than_tl], A[greater_than_tl], alpha[greater_than_tl])
    # only process the mid level part of the signal
        
    Parameters:
        * exc_high (type: array-like matrix of floats) - the high level (larger than the threshold of high level) excitation pattern
        for each ERB (Required)
        
    Returns:
        * specific_loudness_high (type: array-like matrix of floats) - An array with dimensions equal to the exc_high containing the specific loudness of the signal at high levels.
    """

    C = 0.046871 #constant parameter
    specific_loudness_high = C * (exc_high / (1.04 * 10**6))**0.5
    
    return specific_loudness_high

def _get_greater_than_tl(exc_i):
    """
    A function to return if each element of the excitation intensity is greater than the threshold of loud.
        
    Parameters:
        * exc_i (type: array-like matrix of floats) - the input excitation intensity
        
    Returns:
        * le_tq (type: array-like matrix of booleans) - a boolean array with dimensions equal to the input specifying if the excitation intensity is greater than the threshold of loud
    """
    
    lenSig = np.shape(exc_i)[-1]
    g_tl = exc_i[:,:]>np.transpose(np.tile(10**10,(lenSig,39)))
                 
    return g_tl

def _get_less_than_tq(exc_i, tq):
    """
    A function to return if each element of the excitation intensity is less than the threshold of quiet.
    
    Parameters:
        * exc_i (type: array-like matrix of floats) - the input excitation intensity
        * tq (type: array-like matrix of floats) - the threshold of quiet for each ERB.
            
    Returns:
        * le_tq (type: array-like matrix of booleans) - a boolean array with dimensions equal to the input specifying if the excitation intensity is less than the threshold of quiet
    """
    
    if (np.shape(exc_i)!=np.shape(tq)):
        np.transpose(np.tile(exc_i,(np.shape(exc_i)[-1],1)))
    
    le_tq = exc_i<=tq

    return le_tq

def get_excitation_i(input, fs, SPL, gammatoneLength=4096, rectify=False, verbose = False):
    """
    A function to calculate the excitation intensity of the input signal.
    
    Parameters:
        * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
        * fs (type: numerical) - sample frequency of the signal, input. (Required)
        * SPL (type: double) - the sound pressure level (SPL) at 0 dBFS, e.g., the SPL of a sine wave with peaks at amplitude 1 and troughs at amplitude -1. (Required)
        * rectify (type: boolean) - Specifies whether to include half wave rectification, modelling the direction of that the cochlear nerves vibrate. True to include, False to ignore. (Optional; Default = False)
        
    Returns:
        * gtfs (type: numpy array of floats) - array with size ((39,) + np.shape(input)) containing the excitation pattern (in sound intensity) for each ERB of input signal. The excitation pattern for the nth ERB can be accessed with gtfs[n].
    """

    input = np.array(input)
    
    input = v_Pascal(input, SPL)

    excitation = outMidFir(input)
    
    b = gammatone_filter(fs, filterLength=gammatoneLength, normalise = True)
    excitation = filterbank_decomposition(excitation, b, verbose = verbose)
    
    if (rectify):
        excitation = half_rectification(excitation)
    
    excitation = pa_i(excitation)
    excitation = at_normalise(excitation)

    return excitation

def plot_excitation_response(input = None, fs = 44100, outMidFilt = True, nyquist=True, gammatone = True, gammatoneLength=4096, xscale = 'log', yscale = 'log'):
    """
    A function that plots the transfer function of the outer middle ear and each gammatone filter.
    
    Parameters:
        * fs (type: numerical) - the sampling frequency of the signal. (Optional; Default = 44100)
        * outMidFilt (type: boolean) - filter the signal by the outer and middle ear FIR filter. (Optional; Default = True)
        * xscale (type: string) - the scale of the frequency axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
        * yscale (type: string) - the scale of the amplitude axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
    """
    
    if input == None:
        input = np.zeros((np.ceil(fs)))
        input[0] = 1

    if(outMidFilt): input = outMidFir(input)
    if(gammatone):
        b = gammatone_filter(fs, filterLength=gammatoneLength, normalise = True)
        input = filterbank_decomposition(input, b)
        #input = holdsworth_gammatone_filter(input,fs)
        numPlot = range(np.shape(input)[0])
    else:
        numPlot = (0,)
        input = input.reshape(1,len(input))
        
    for i in numPlot:
        utils.plot_fft(input[i],xscale, yscale, nyquist=True, show=False)

    plt.show()

    return

def get_modulation_i(input, fs):
    """
    A function to calculate the modulation intensity of the input intensity signal. The function implements a filter bank of bandpass filters with cut off frequencies ranging from 0.25 to 16 Hz.
    
    Parameters:
        * input (type: array-like matrix of floats) - the input intensity signal.
            E.g., use get_excitation_i() to obtain excitation intensity and use as input.
        * fs (type: numerical) - sampling frequency of input signal
        
    Returns:
        * y (type: numpy array of floats) - array with size ((10,) + np.shape(input)) containing the modulation intensity of the signal at each modulation filter. The modulation intensity for the nth filter can be accessed with y[n].
    """
    
    input = np.array(input)
    b = fir_antialias(fs)
    input_lp = sp.lfilter(b,(1),input_fr)
    input_ds = downsample(input_lp, fs)
    fc = np.array(utils.exp_sequence(-2,4,10))
    bw = fc/2
    y = filterbank_decomposition(input_ds, fs, fc, bw)

    return y

def outMidFir(input):
    """
    A function to filter the input signal with the response of the outer and middle ear.
    
    Parameters:
        * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
        
    Returns:
        * y (type: numpy array of floats) - array with dimensions equal to the input signal filtered by the response of the outer and middle ear.
    """
    
    input = np.array(input)
    b = utils.load_outMidFir_coeff()
    y = sp.fftconvolve(input, b)
    
    return y

def gammatone_filter(fs, filterLength=4096, normalise = True):
    """
    A function to return the filter coefficients of 39 different gammatones with fc of ERBs 1 to 39.
    
    Parameters:
        * fs (type: numerical) - sample frequency of the signal, input. (Required)
        * filterLength (type: int) - length of the filters (Optional; Default = True)
        * normalise (type: boolean) - normalise peaks to united (Optional; Default = True)
        
    Returns:
        * b (type: numpy array of floats) - array with size ((39),np.shape(input)) containing the impulse response of the signal at each gammatone filter. The response at the bth gammatone filter can be accessed by y[n].
    """
    
    nerbs = 39
    erbs = np.array(range(1,nerbs+1))
    fc = (10**(erbs/21.366)-1)/0.004368
    bw = 24.673 * (1 + 0.004368*fc)
    N = 4
    t = 1.0*np.array(range(filterLength))/fs
    
    b = np.zeros((nerbs,filterLength))

    for i in range(39):
        b[i] = t**(N-1) * np.exp(-2*np.pi*bw[i]*t) * np.cos(2*np.pi*fc[i]*t)
    
        if normalise:
            A = 1/(np.abs(fft.fft(b[i], 4*filterLength))).max()
            b[i] *= A

    return b

def _holdsworth_gammatone_filter(input, fs):
    """
    Feed the input through a bank of gammatone filters, Holdsworth implementation. Not sure how great this works, don't use.
    
    Parameters:
        * input (type: array-like matrix of floats) - signal. (Required)
        * fs (type: numerical) - sample frequency of the signal, input. (Required)
        
    Returns:
        * out (type: numpy array of floats) - array with size ((39),np.shape(input)) containing the output of the signal at each gammatone filter.
    """
    input = np.array(input)
    input = input + 0j

    T = 1.0/fs

    ERBs = np.array(range(1,40))
    f0 = (10**(ERBs/21.4)-1)/4.37e-3

    inLen = len(input)
    b = 24.673 * (1 + 0.004368*f0)
    k = np.array(range(inLen)) + 0j
    out = np.zeros((39,inLen))

    for erb in range(39):
        
        zArr = input*np.exp(-2*np.pi*1j*f0[erb]*k*T)
        wArr = np.zeros((inLen+1))

        for i in range(1,inLen+1):
            wArr[i] = wArr[i-1] + (1 - np.exp(-2*np.pi*b[erb]*T))*(zArr[i-1] - wArr[i-1])

        out[erb] = (wArr[1:]*np.exp(2*np.pi*1j*f0[erb]*k*T)).real

    return out


def v_Pascal(input, SPL):
    """
    A function to convert a signal, normalised to an amplitude range of -1 to 1, to a signal represented in pressure (units: Pascal).
    
    Parameters:
        * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
        * SPL (type: double) - the sound pressure level (SPL) at 0 dBFS, e.g., the SPL of a sine
            wave with peaks at amplitude 1 and troughs at amplitude -1. (Required)
            
    Returns:
        * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
            as a pressure signal.
            
    NOTE If input signal contains one or more zeros, a divide by zero warning will be outputted. This warning can be ignored.
    """
    
    input = np.array(input)
    sigRMS = rms(input)
    y = input*0.00002*10**(SPL/20.0)
    
    return y

def rms(input):
    
    sigRMS = np.sqrt(np.mean(input**2))
    
    return sigRMS

def pa_i(input, C=406):
    """
    A function to convert a pressure signal (unit: Pascal) to an intensity signal.
        
    Parameters:
        * input (type: array-like matrix of floats) - pressure signal (unit: Pascal) (Required)
        * C (type: double) - the acoustic impedance of the air (Optional; Default = 406)
        
    Returns:
        * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
        as a pressure signal.
    """
    
    input = np.array(input)
    y = (input**2) / C
    
    return y

def at_normalise(input):
    """
    A function to normalise an intensity signal with the audibility threshold.
    
    Parameters:
        * input (type: array-like matrix of floats) - intensity signal (unit: Pascal) (Required)
    
    Returns:
        * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input normalised
            with the audibility threshold.
    """
    
    input = np.array(input)
    y = input / (1*(10**-12))
    
    return y

def downsample(input, factor=100):
    """
    A function to downsample a signal, input, with sampling frequency, fs, by a downsample factor of factor.
    
    NOTE: It is advised to use the fir_antialias() function before downsampling to remove any high frequencies
        which would otherwise represented as low frequencies due to aliasing.
        
    Parameters:
        * input (type: array-like matrix of floats) - input signal. (Required)
        * factor - downsample factor (Optional; Default = 100)
        
    Returns:
        * output (type: numpy array of floats) - array with outer dimensions equivalent to the to the input, and
            inner dimension equal to np.floor(lenIn / factor) where lenIn is the length of the inner dimension.
    """
    
    input = np.array(input)
    shapeIn = np.shape(input)
    nDim = np.shape(shapeIn)
    lenIn = shapeIn[nDim[0]-1]
    lenOut = np.floor(lenIn / factor)
    n = np.linspace(0,lenIn,lenOut, endpoint=False).astype(np.int)
    output = input[...,n]

    return output

def half_rectification(input):
    """
    A function which performs a half-wave rectification on the input signal.
    
    Parameters:
        * input (type: array-like matrix of floats) - input signal. (Required)
        
    Returns:
        * y (type: numpy array of floats) - an array with dimensions of input containing the half-wave rectification of
            input.
    """
    
    
    y = np.array(input)
    y[y<0] = 0
    
    return y

def filterbank_decomposition(input, b, a = None, verbose = False):
    """
    A function to run the input through a bandpass filter bank with parameters defined by the b and a coefficints.
    
    Parameters:
        * input (type: array-like matrix of floats) - input signal. (Required)
        * b (type: array-like matrix of floats) - the b coefficients of each filter in shape b[numOfFilters][numOfCoeffs]. (Required)
        * a (type: array-like matrix of floats) - the a coefficients of each filter in shape a[numOfFilters][numOfCoeffs].
            If not specified, the filter is treated as a FIR filter. (Optional; Default = 1)
        * verbose (type: boolean) - determines whether to display the current subroutine/progress of the procedure.
            (Optional; Default = False)
        
    Returns:
    * y (type: numpy array of floats) - an array with inner dimensions equal to that of the input and outer dimension equal to
        the length of fc (i.e. the number of bandpass filters in the bank) containing the outputs to each filter. The output 
        signal of the nth filter can be accessed using y[n].
    """
    
    input = np.array(input)
    bShape = np.shape(b)
    nFilters = bShape[0]
    lengthFilter = bShape[1]
    
    shape = (nFilters,) + (np.shape(input))
    shape = np.array(shape[0:])
    shape[-1]=shape[-1]+lengthFilter-1
    y = np.zeros((shape))
    
    if(verbose): sys.stdout.write("Running filterbank filterbank_decomposition.\n")
    for i in range(nFilters):
        if(verbose):
            sys.stdout.write("\r" + str(int(np.round(100.0*i/nFilters))) + "% complete.")
            sys.stdout.flush()
        x = np.array(input)
        y[i] = sp.fftconvolve(x,b[i])

    if(verbose): sys.stdout.write("\n")

    return y

def get_Modulation_filter_bank_coefficients(fs, fc, bw, order = 1024, verbose = False):
    """
    A function which returns the b and a coefficients of the modulation filter bank of length equal to the length of fc. Each
    bandpass filter is designed by defining the centre frequencies, fc, and bandwidths, bw.
        
    Parameters:
        * fs (type: numerical) - the sampling frequency of the input signal. (Required)
        * fc (type: array-like vector of floats) - the centre off frequencies (unnormalised) of each bandpass filter.
            The length of this vector determines the number of filters in the bank.
        * bw (type: array-like vector of floats) - the bandwidths (unnormalised) of each bandpass filter. Must be equal
            to or more than the length of fc. If the length is more, all elements exceeding the length of fc - 1 will be
            ignored.
        * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 1024)
        * verbose (type: boolean) - determines whether to display the current subroutine/progress of the procedure.
            (Optional; Default = False)
        
    Returns:
        * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap)
    """


    nFreqs = len(fc)

    if(verbose): sys.stdout.write("Running frequency filterbank_decomposition.")

    for i in range(nFreqs):
        if(verbose):
            sys.stdout.write("\r" + str(int(np.round(100.0*i/nFreqs))) + "% complete.")
            sys.stdout.flush()
        low = fc[i]-bw[i]/2;
        high = fc[i]+bw[i]/2;
        b = fir_bandpass(low, high, fs, order, verbose)

    return b

def temporal_integration(input, fs, type='STL'):
    """
    Convert from instantaneous loudness to STL (type='STL'), or STL to LTL (type='LTL').
    
    Parameters:
        * input (type: array-like matrix of floats) - signal. Should be instantaneous loudness if type='STL', or STL if type='LTL' (Required)
        * fs (type: numerical) - sample frequency of the signal, input. (Required)
        * type (type: string) - 'STL' or 'LTL'
        
    Returns:
        * Output (type: array-like matrix of floats) - output signal.
    """
    Ti = 1.0/fs
    prevOutput = 0;
    output = np.zeros(np.shape(input))
    if(type=='STL'):
        Ta = 0.0217
        Tr = 0.0495
    elif(type=='LTL'):
        Ta = 0.0995
        Tr = 1.9995
    else:
        str = "Type argument " + type + " is unknown"
        raise ValueError(str)

    a = 1 - np.exp(-Ti/Ta)
    r = 1 - np.exp(-Ti/Tr)

    for i in range(len(input)):
        if(input[i] > prevOutput):
            output[i] = a*input[i] + (1-a)*prevOutput
        else:
            output[i] = r*input[i] + (1-r)*prevOutput
        prevOutput = output[i]

    return output

def exp_smoothing(fs, tc = 0.01):
    """
    Designs an exponential filter, y[n] = alpha*x[n] - -(1-alpha)*y[n-1], where alpha = T/(tc+T), where T
    is the inverse of the sampling frequency and time constant, tc, is specified.
    
    Parameters:
        * fs (type: numerical) - sampling frequency of the signal to be filtered (Required)
        * tc (type: numerical) - the time constant of the filter. (Optional; Default = 0.01)
        
    Returns:
        * b (type: numerical) - the coefficient of x[n]. Equal to alpha.
        * a (type: numpy array of floats) - an array of size 2 of the coefficients of y[n] (a[0] = 1) and y[n-1]
    """
    
    T = 1.0/fs
    alpha = T/(tc + T)
    b = [alpha]
    a = [1, -(1-alpha)]
    
    return b, a

def antialias_fir(fs, fc=100, order=64):
    """
    A function which returns the b coefficients for a lowpass fir filter with specified requirements.
    Made specifically to remove aliasing when downsampling, but can be used for any application that
    requires a lowpass filter.
        
    Parameters:
        * fs (type: numerical) - sampling frequency of the signal to be filtered (Required)
        * fc (type: numerical) - unnormalised cut off frequency of the filter (Optional; Default = 100)
        * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 64)
        
    Returns:
        * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap)
    """
    
    nyquist = 0.5*fs
    fcNorm = fc/nyquist
    b = sp.firwin(order+1, fcNorm)
    
    return b

def fir_bandpass(low, high, fs, order = 4096, verbose = False):
    """
    A function which returns the b coefficients for a bandpass fir filter with specified requirements.
        
    Parameters:
        * low - the lower cutoff frequency of the filter. (Required)
        * high - the upper cutoff frequency of the filter. (Required)
        * fs (type: numerical) - sampling frequency of the signal to be filtered. (Required)
        * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 1024)
        * verbose (type: boolean) - determines whether to display info on the current subroutine
        of the procedure. (Optional; Default = False)
        
    Returns:
        * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap).
    """
    
    nyquist = 0.5*fs
    lowNorm = low/nyquist
    highNorm = high/nyquist
    if(verbose): sys.stdout.write("Low: " + str(lowNorm) + "; High: " + str(highNorm))
    b = sp.firwin(order+1, [lowNorm, highNorm], pass_zero=False)
    
    return b