view cortical_model.py @ 1:dc43033a2c20

* Added .dat files containing parameters and filter coefficients
* In cortical_model: added get_specific_loudness() function and dependencies.
* In cortical_model & in utils: fixed bugs created from changing variable and function names before last submit
* In cortical_model: changed plot_excitation_tf (now named plot_excitation_response) to allow an input signal to be optionally given to be analysed.

TODO:
* In cortical_model: comment
* In cortical_model: add partial loudness function

"""
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.
    
    TO DO
    """
    lenSig = np.shape(exc_i)[-1] #length signal
    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 parameter
    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_low(exc_low, G, tq, A, alpha):
    
    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):
    
    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):
    """
    
    """

    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)

    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 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 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 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