annotate cortical_model.py @ 3:2d80632482b3

* Added exponential smoothing filtering into get_excitation_i() + dependencies
* Broke up get_specific_loudness to allow user to return loudness parameters tq, A, G and alpha using get_specific_loudness_parameters().
* Commented all functions

"""
A module used for auditory analysis.

Models currently implemented:
    * Frequency-modulation analysis model based on the human auditory system.

Model implementations in progress:
    * Glasberg and Moore model
    
Packaged dependencies:
    * utils.py and/or utils.pyc
    * erb.dat
    * outMidFir.dat
    * tq.dat
    
External dependencies:
    * scipy
    * numpy
    * copy
    * matplotlib
"""

import utils
import scipy.signal as sp
import scipy.fftpack as fft
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt

def get_specific_loudness(exc_i, inc_loud_region = False):
    """
    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)
    
    TO DO
    """
    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 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.
    
    e.g.,
    
    # 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.047 #constant parameter
    specific_loudness_low = C * ((2*exc_low)/(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.
        
    e.g.,
        
    # 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.047 #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.
        
    e.g.,
        
    # 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.047 #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, 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)
    inputOMFir = outMidFir(input)
    inputPa = v_Pascal(inputOMFir, SPL)
    gtfs = gamma_tone_filter(inputPa, fs, verbose = verbose)
    if (rectify):
        gtfs = half_rectification(gtfs)
    gtfs = pa_i(gtfs)
    gtfs = at_normalise(gtfs)
    b,a = exp_smoothing(fs)
    gtfs = sp.lfilter(b,a,gtfs)

    return gtfs

def plot_excitation_response(input = 0, fs = 44100, outMidFilt = True, 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')
    
    Returns:
        * y (type: numpy array of floats) - array with size ((39,np.ceil(fs))) containing the impulse response of the
            signal at each gammatone filter (and optionally, including the outer middle ear response). The response at the
            nth gammatone filter can be accessed by y[n].
    """
    
    if ((np.shape(input)==())):
        input = np.zeros(np.ceil(fs))
        input[0] = 1
    if(outMidFilt): input = outMidFir(input)
    y = gamma_tone_filter(input, fs)
    
    for i in range(np.shape(y)[0]):
        utils.plot_fft(y[i],xscale, yscale, False)

    plt.show()

    return fft.fft(y)

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 = 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.lfilter(b, (1), input)

    return y

def gamma_tone_filter(input, fs, verbose = False):
    """
    A function to filter to decompose the input signal into 39 different gammatones filtered signals modelling the ERBs.
    
    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)
        
    Returns:
        * y (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 nth gammatone filter can be accessed by y[n].
    """
    
    input = np.array(input)
    fc, bw = utils.load_erb_data()
    y = decomposition(input,fs,fc,bw, verbose = verbose)
    
    return y

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.
    """
    
    input = np.array(input)
    y = np.sign(input)*(0.00002*10**(np.log10(np.abs(input))+SPL/20))
    
    return y

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 decomposition(input, fs, fc, bw, order=1024, verbose = False):
    """
    A function to run the input filter through a bandpass 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:
        * input (type: array-like matrix of floats) - input signal. (Required)
        * 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:
    * 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)
    nFreqs = len(fc)
    shape = (nFreqs,) + np.shape(input)
    shape = shape[0:]
    y = np.zeros(shape)
    
    if(verbose): print "Running frequency decomposition."
    for i in range(nFreqs):
        if(verbose): print str(100.0*i/nFreqs) + "% complete."
        low = fc[i]-bw[i]/2;
        high = fc[i]+bw[i]/2;
        if(verbose): print "Low: " + str(low) + "; High: " + str(high)
        b = fir_bandpass(low, high, fs, order, verbose)
        x = deepcopy(input)
        y[i] = sp.lfilter(b,(1),x)
    
    return y

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 = 1024, 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 the current progress (or 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): print "Low: " + str(lowNorm) + "; High: " + str(highNorm)
    b = sp.firwin(order+1, [lowNorm, highNorm], pass_zero=False)
    
    return b