c@0: """
c@0: A module used for auditory analysis.
c@0: 
c@0: Models currently implemented:
c@0:     * Frequency-modulation analysis model based on the human auditory system.
c@0: 
c@0: Model implementations in progress:
c@0:     * Glasberg and Moore model
c@0:     
c@0: Packaged dependencies:
c@0:     * utils.py and/or utils.pyc
c@0:     * erb.dat
c@0:     * outMidFir.dat
c@0:     
c@0: External dependencies:
c@0:     * scipy
c@0:     * numpy
c@0:     * copy
c@0:     * matplotlib
c@0: """
c@0: 
c@0: import utils
c@0: import scipy.signal as sp
c@0: import numpy as np
c@0: from copy import deepcopy
c@0: import matplotlib.pyplot as plt
c@0: 
c@0: def get_partial_loudness():
c@0:     """
c@0:     A function to calculate the partial loudness of an excitation pattern at given ERB.
c@0:     
c@0:     TO DO
c@0:     """
c@0:     
c@0:     
c@0:     return
c@0: 
c@0: def get_excitation_i(input, fs, SPL, rectify=False):
c@0:     """
c@0:     A function to calculate the excitation intensity of the input signal.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
c@0:         * fs (type: numerical) - sample frequency of the signal, input. (Required)
c@0:         * SPL (type: double) - the sound pressure level (SPL) at 0 dBFS, e.g., the SPL of a sine
c@0:             wave with peaks at amplitude 1 and troughs at amplitude -1. (Required)
c@0:         * rectify (type: boolean) - Specifies whether to include half wave rectification, modelling the direction
c@0:             of that the cochlear nerves vibrate.
c@0:             True to include, False to ignore. (Optional; Default = False)
c@0:         
c@0:     Returns:
c@0:         * gtfs (type: numpy array of floats) - array with size ((39,) + np.shape(input)) containing the excitation
c@0:             pattern (in sound intensity) for each ERB of input signal. The excitation pattern for the nth ERB can
c@0:             be accessed with gtfs[n].
c@0:     """
c@0: 
c@0:     input = np.array(input)
c@0:     inputOMFir = outMidFir(input)
c@0:     inputPa = v_Pascal(inputOMFir, SPL)
c@0:     gtfs = gamma_tone_filter(inputPa, fs)
c@0:     if (rectify):
c@0:         gtfs = half_rectification(gtfs)
c@0:     gtfs = pa_i(gtfs)
c@0:     gtfs = at_normalise(gtfs)
c@0: 
c@0:     return gtfs
c@0: 
c@0: def plot_excitation_tf(fs = 44100, outMidFilt = True, xscale = 'log', yscale = 'log'):
c@0:     """
c@0:     A function that plots the transfer function of the outer middle ear and each gammatone filter.
c@0:     
c@0:     Parameters:
c@0:         * fs (type: numerical) - the sampling frequency of the signal. (Optional; Default = 44100)
c@0:         * outMidFilt (type: boolean) - filter the signal by the outer and middle ear FIR filter. (Optional; Default = True)
c@0:         * xscale (type: string) - the scale of the frequency axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
c@0:         * yscale (type: string) - the scale of the amplitude axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
c@0:     
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with size ((39,np.ceil(fs))) containing the impulse response of the
c@0:             signal at each gammatone filter (and optionally, including the outer middle ear response). The response at the
c@0:             nth gammatone filter can be accessed by y[n].
c@0:     """
c@0:         
c@0:     input = np.zeros(np.ceil(fs))
c@0:     input[0] = 1
c@0:     if(outMidFilt): input = outMidFir(input)
c@0:     y = gamma_tone_filter(input, fs)
c@0:     
c@0:     for i in range(np.shape(y)[0]):
c@0:         utils.plot_fft(y[i],xscale, yscale, False)
c@0: 
c@0:     plt.show()
c@0: 
c@0:     return y
c@0: 
c@0: def get_modulation_i(input, fs):
c@0:     """
c@0:     A function to calculate the modulation intensity of the input intensity signal. The function implements a
c@0:     filter bank of bandpass filters with cut off frequencies ranging from 0.25 to 16 Hz.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - the input intensity signal.
c@0:             E.g., use get_excitation_i() to obtain excitation intensity and use as input.
c@0:         * fs (type: numerical) - sampling frequency of input signal
c@0:         
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with size ((10,) + np.shape(input)) containing the
c@0:         modulation intensity of the signal at each modulation filter. The modulation intensity for the nth filter can
c@0:         be accessed with y[n].
c@0:     """
c@0:     
c@0:     
c@0:     input = np.array(input)
c@0:     b = fir_antialias(fs)
c@0:     input_lp = sp.lfilter(b,(1),input_fr)
c@0:     input_ds = downsample(input_lp, fs)
c@0:     fc = np.array(utils.exp_sequence(-2,4,10))
c@0:     bw = fc/2
c@0:     y = decomposition(input_ds, fs, fc, bw)
c@0: 
c@0:     return y
c@0: 
c@0: def outMidFir(input):
c@0:     """
c@0:     A function to filter the input signal with the response of the outer and middle ear.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
c@0:         
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with dimensions equal to the input signal filtered by the response of
c@0:             the outer and middle ear.
c@0:     """
c@0:     
c@0:     input = np.array(input)
c@0:     b = utils.load_outMidFir_coeff()
c@0:     y = sp.lfilter(b, (1), input)
c@0: 
c@0:     return y
c@0: 
c@0: def gamma_tone_filter(input, fs):
c@0:     """
c@0:     A function to filter to decompose the input signal into 39 different gammatones filtered signals modelling the ERBs.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
c@0:         * fs (type: numerical) - sample frequency of the signal, input. (Required)
c@0:         
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with size ((39),np.shape(input)) containing the impulse response of the
c@0:             signal at each gammatone filter. The response at the nth gammatone filter can be accessed by y[n].
c@0:     """
c@0:     
c@0:     input = np.array(input)
c@0:     fc, bw = utils.load_erb_data()
c@0:     y = decomposition(input,fs,fc,bw)
c@0:     
c@0:     return y
c@0: 
c@0: def v_Pascal(input, SPL):
c@0:     """
c@0:     A function to convert a signal, normalised to an amplitude range of -1 to 1, to a signal represented in pressure (units: Pascal).
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
c@0:         * SPL (type: double) - the sound pressure level (SPL) at 0 dBFS, e.g., the SPL of a sine
c@0:             wave with peaks at amplitude 1 and troughs at amplitude -1. (Required)
c@0:             
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
c@0:             as a pressure signal.
c@0:     """
c@0:     
c@0:     input = np.array(input)
c@0:     y = np.sign(input)*(0.00002*10**(np.log10(np.abs(input))+dBFS/20))
c@0:     
c@0:     return y
c@0: 
c@0: def pa_i(input, C=406):
c@0:     """
c@0:     A function to convert a pressure signal (unit: Pascal) to an intensity signal.
c@0:         
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - pressure signal (unit: Pascal) (Required)
c@0:         * C (type: double) - the acoustic impedance of the air (Optional; Default = 406)
c@0:         
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
c@0:         as a pressure signal.
c@0:     """
c@0:     
c@0:     input = np.array(input)
c@0:     y = (input**2) / C
c@0:     
c@0:     return y
c@0: 
c@0: def at_normalise(input):
c@0:     """
c@0:     A function to normalise an intensity signal with the audibility threshold.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - intensity signal (unit: Pascal) (Required)
c@0:     
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input normalised
c@0:             with the audibility threshold.
c@0:     """
c@0:     
c@0:     input = np.array(input)
c@0:     y = input / 1*(10**12)
c@0:     
c@0:     return
c@0: 
c@0: def downsample(input, factor=100):
c@0:     """
c@0:     A function to downsample a signal, input, with sampling frequency, fs, by a downsample factor of factor.
c@0:     
c@0:     NOTE: It is advised to use the fir_antialias() function before downsampling to remove any high frequencies
c@0:         which would otherwise represented as low frequencies due to aliasing.
c@0:         
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - input signal. (Required)
c@0:         * factor - downsample factor (Optional; Default = 100)
c@0:         
c@0:     Returns:
c@0:         * output (type: numpy array of floats) - array with outer dimensions equivalent to the to the input, and
c@0:             inner dimension equal to np.floor(lenIn / factor) where lenIn is the length of the inner dimension.
c@0:     """
c@0:     
c@0:     input = np.array(input)
c@0:     shapeIn = np.shape(input)
c@0:     nDim = np.shape(shapeIn)
c@0:     lenIn = shapeIn[nDim[0]-1]
c@0:     lenOut = np.floor(lenIn / factor)
c@0:     n = np.linspace(0,lenIn,lenOut, endpoint=False).astype(np.int)
c@0:     output = input[...,n]
c@0: 
c@0:     return output
c@0: 
c@0: def antialias_fir(fs, fc=100, order=64):
c@0:     """
c@0:     A function which returns the b coefficients for a lowpass fir filter with specified requirements.
c@0:     Made specifically to remove aliasing when downsampling, but can be used for any application that
c@0:     requires a lowpass filter.
c@0:     
c@0:     Parameters:
c@0:         * fs (type: numerical) - sampling frequency of the signal to be filtered (Required)
c@0:         * fc (type: numerical) - unnormalised cut off frequency of the filter (Optional; Default = 100)
c@0:         * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 64)
c@0:         
c@0:     Returns:
c@0:         * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap)
c@0:     """
c@0:     
c@0:     nyquist = 0.5*fs
c@0:     fcNorm = fc/nyquist
c@0:     b = sp.firwin(order+1, fcNorm)
c@0:     
c@0:     return b
c@0: 
c@0: def half_rectification(input):
c@0:     """
c@0:     A function which performs a half-wave rectification on the input signal.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - input signal. (Required)
c@0:         
c@0:     Returns:
c@0:         * y (type: numpy array of floats) - an array with dimensions of input containing the half-wave rectification of
c@0:             input.
c@0:     """
c@0:     
c@0:     
c@0:     y = np.array(input)
c@0:     y[y<0] = 0
c@0:     
c@0:     return y
c@0: 
c@0: def decomposition(input, fs, fc, bw, order=1024, verbose = False):
c@0:     """
c@0:     A function to run the input filter through a bandpass filter bank of length equal to the length of fc. Each
c@0:     bandpass filter is designed by defining the centre frequencies, fc, and bandwidths, bw.
c@0:     
c@0:     Parameters:
c@0:         * input (type: array-like matrix of floats) - input signal. (Required)
c@0:         * fs (type: numerical) - the sampling frequency of the input signal. (Required)
c@0:         * fc (type: array-like vector of floats) - the centre off frequencies (unnormalised) of each bandpass filter.
c@0:             The length of this vector determines the number of filters in the bank.
c@0:         * bw (type: array-like vector of floats) - the bandwidths (unnormalised) of each bandpass filter. Must be equal
c@0:             to or more than the length of fc. If the length is more, all elements exceeding the length of fc - 1 will be
c@0:             ignored.
c@0:         * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 1024)
c@0:         * verbose (type: boolean) - determines whether to display the current subroutine/progress of the procedure.
c@0:             (Optional; Default = False)
c@0:         
c@0:     Returns:
c@0:     * y (type: numpy array of floats) - an array with inner dimensions equal to that of the input and outer dimension equal to
c@0:         the length of fc (i.e. the number of bandpass filters in the bank) containing the outputs to each filter. The output 
c@0:         signal of the nth filter can be accessed using y[n].
c@0:     """
c@0:     
c@0:     
c@0:     input = np.array(input)
c@0:     nFreqs = len(fc)
c@0:     shape = (nFreqs,) + np.shape(input)
c@0:     shape = shape[0:]
c@0:     y = np.zeros(shape)
c@0:     
c@0:     if(verbose): print "Running frequency decomposition."
c@0:     for i in range(nFreqs):
c@0:         if(verbose): print str(100.0*i/nFreqs) + "% complete."
c@0:         low = fc[i]-bw[i]/2;
c@0:         high = fc[i]+bw[i]/2;
c@0:         if(verbose): print "Low: " + str(low) + "; High: " + str(high)
c@0:         b = fir_bandpass(low, high, fs, order, verbose)
c@0:         x = deepcopy(input)
c@0:         y[i] = sp.lfilter(b,(1),x)
c@0:     
c@0:     return y
c@0: 
c@0: def fir_bandpass(low, high, fs, order = 1024, verbose = False):
c@0:     """
c@0:     A function which returns the b coefficients for a bandpass fir filter with specified requirements.
c@0:         
c@0:     Parameters:
c@0:         * low - the lower cutoff frequency of the filter. (Required)
c@0:         * high - the upper cutoff frequency of the filter. (Required)
c@0:         * fs (type: numerical) - sampling frequency of the signal to be filtered. (Required)
c@0:         * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 1024)
c@0:         * verbose (type: boolean) - determines whether to display the current progress (or info on the current subroutine)
c@0:         of the procedure. (Optional; Default = False)
c@0:         
c@0:     Returns:
c@0:         * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap).
c@0:     """
c@0:     
c@0:     nyquist = 0.5*fs
c@0:     lowNorm = low/nyquist
c@0:     highNorm = high/nyquist
c@0:     if(verbose): print "Low: " + str(lowNorm) + "; High: " + str(highNorm)
c@0:     b = sp.firwin(order+1, [lowNorm, highNorm], pass_zero=False)
c@0:     
c@0:     return b