changeset 4:29cd3e735c4c

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

--- a/cortical_model.py	Tue Feb 04 14:46:06 2014 +0000
+++ b/cortical_model.py	Sat Feb 22 13:53:55 2014 +0000
@@ -251,7 +251,8 @@
     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)
@@ -261,7 +262,7 @@
 
     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.
     
@@ -270,25 +271,28 @@
         * 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):
     """
@@ -334,7 +338,7 @@
 
     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.
     
@@ -347,11 +351,54 @@
             signal at each gammatone filter. The response at the nth gammatone filter can be accessed by y[n].
     """
     
+    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)
-    fc, bw = utils.load_erb_data()
-    y = decomposition(input,fs,fc,bw, verbose = verbose)
-    
-    return y
+    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):
     """
@@ -451,13 +498,49 @@
     
     return y
 
-def decomposition(input, fs, fc, bw, order=1024, verbose = False):
+def decomposition(input, b, a = None, 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.
+    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:
         * 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.
@@ -469,28 +552,22 @@
             (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].
+        * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap)
     """
-    
-    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
+
+    return b
 
 def exp_smoothing(fs, tc = 0.01):
     """
@@ -534,7 +611,7 @@
     
     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.
         

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/run.py	Sat Feb 22 13:53:55 2014 +0000
@@ -0,0 +1,20 @@
+import cortical_model as cm
+import numpy as np
+import utils
+
+t = np.array(range(44100))
+sine = np.sin(2*np.pi*t*1000/44100)
+sine1 = np.sin(2*np.pi*t*100/44100)
+x = np.append(sine, sine1)
+
+exc_i = cm.get_excitation_i(x, 44100, 90, verbose = True)
+sl = cm.get_specific_loudness(exc_i)
+
+for i in range(39):
+    print "ERB " + str(i)
+    utils.plot_waveform(sl[i])
+
+print "Loudness"
+loudness = np.sum(sl,0)
+
+utils.plot_waveform(loudness)
\ No newline at end of file

--- a/utils.py	Tue Feb 04 14:46:06 2014 +0000
+++ b/utils.py	Sat Feb 22 13:53:55 2014 +0000
@@ -233,8 +233,8 @@
         NONE
     """
 
-    if(xscale.lower() == "log"): plt.gca().set_xscale('log')
-    if(yscale.lower() == "log"): plt.gca().set_yscale('log')
+    plt.gca().set_xscale(xscale)
+    plt.gca().set_yscale(yscale)
     
     x = np.array(x)
     fftx = np.absolute(fft.fft(x))