comparison cortical_model.py @ 4:29cd3e735c4c

* Added run.py for testing
* Changed frequency decomposition to bank of gammatones

     """
 
     input = np.array(input)
     inputOMFir = outMidFir(input)
     inputPa = v_Pascal(inputOMFir, SPL)
-    gtfs = gamma_tone_filter(inputPa, fs, verbose = verbose)
+    b = gamma_tone_filter(fs)
+    gtfs = decomposition(inputPa, b, 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'):
+def plot_excitation_response(input = None, fs = 44100, outMidFilt = True, gammatone = 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))
+    """
+    
+    if input == None:
+        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)
+    if(gammatone):
+        b = gamma_tone_filter(fs)
+        input = decomposition(input, b)
+        #input = holdsworthGamma(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, False)
 
     plt.show()
 
-    return fft.fft(y)
+    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.

     b = utils.load_outMidFir_coeff()
     y = sp.lfilter(b, (1), input)
 
     return y
 
-def gamma_tone_filter(input, fs, verbose = False):
+def gamma_tone_filter(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)

     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
+    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
+    filterLength = 4096
+    t = 1.0*np.array(range(filterLength))/fs
+
+    gain=((1.019*bw*(2.0*np.pi)/float(fs))**4)/6.0
+    
+    PI = N * np.arctan(1)
+    
+    b = np.zeros((nerbs,filterLength))
+
+    for i in range(39):
+        b[i] = gain[i] * t**(N-1) * fs**(N-1) * np.exp(-2*np.pi*bw[i]*t) * np.cos(2*np.pi*fc[i]*t)
+    
+    return b
+
+def holdsworthGamma(input, fs):
+    """
+    
+    """
+    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).
     

     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
+def 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]
+    
+    if a == None:
+        a = np.ones(nFilters)
+    
+    shape = (nFilters,) + np.shape(input)
+    shape = shape[0:]
+    y = np.zeros(shape)
+    
+    if(verbose): print "Running decomposition."
+    for i in range(nFilters):
+        if(verbose): print str(100.0*i/nFilters) + "% complete."
+        x = deepcopy(input)
+        y[i] = sp.lfilter(b[i],a[i],x)
+    
+    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:
-        * input (type: array-like matrix of floats) - input signal. (Required)
+        
+    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

         * 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)
+        * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap)
+    """
+
+
     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
+
+    return b
 
 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.

     fcNorm = fc/nyquist
     b = sp.firwin(order+1, fcNorm)
     
     return b
 
-def fir_bandpass(low, high, fs, order = 1024, verbose = False):
+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)