comparison cortical_model.py @ 0:5609fd93e935

First commit.

+"""
+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
+    
+External dependencies:
+    * scipy
+    * numpy
+    * copy
+    * matplotlib
+"""
+
+import utils
+import scipy.signal as sp
+import numpy as np
+from copy import deepcopy
+import matplotlib.pyplot as plt
+
+def get_partial_loudness():
+    """
+    A function to calculate the partial loudness of an excitation pattern at given ERB.
+    
+    TO DO
+    """
+    
+    
+    return
+
+def get_excitation_i(input, fs, SPL, rectify=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)
+    if (rectify):
+        gtfs = half_rectification(gtfs)
+    gtfs = pa_i(gtfs)
+    gtfs = at_normalise(gtfs)
+
+    return gtfs
+
+def plot_excitation_tf(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].
+    """
+        
+    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 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):
+    """
+    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)
+    
+    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))+dBFS/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
+
+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