diff pyphon.py @ 9:ef62cca13f36

* Added STL and LTL functionality
* Changed cortical_model name to PyPhon because it's a funny pun!
* Added wavread2 to utils, lets you read 24bit wavefiles
* Added sound to utils, listen to your vector like in Matlab!
* changed testElc to equalLoudness

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pyphon.py	Tue Apr 22 00:55:19 2014 +0100
@@ -0,0 +1,927 @@
+"""
+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
\ No newline at end of file