Mercurial > hg > cm

--- a/cortical_model.py	Sun Feb 02 22:48:18 2014 +0000
+++ b/cortical_model.py	Tue Feb 04 14:46:06 2014 +0000
@@ -29,18 +29,21 @@

 def get_specific_loudness(exc_i, inc_loud_region = False):
     """
-    A function to calculate the specific loudness of the excitation patterns at each ERB.
+    A function to calculate the specific loudness of the excitation patterns at each ERB. Specific loudness is calculated for three regions of
+    signal intensity, low level, mid level and high level. The high level processing is optional, and by default, ignored. When ignored the high level
+    region is processed equivalent to that of the mid region.
+
+    Parameters:
+        * exc_i (type: array-like of floats) - an array of shape (39, lenSig) of the excitation intensity of the signal. The outer dimension determines the
+            ERB channel. (Required)
+        * inc_loud_region (type: boolean) - Specifies whether to process the high levels differently to the mid levels. (Optional; Default = False)

     TO DO
     """
+    exc_i = np.array(exc_i)
+
     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
+    tq, G, A, alpha = get_specific_loudness_parameters(lenSig)
     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
@@ -62,7 +65,61 @@

     return specific_loudness

+def get_specific_loudness_parameters(lenSig = 1):
+    """
+    Loads and returns the specific loudness values. If lenSig is specified, the parameters are shaped equal to the excitation signal to allow for
+    matrix elementwise operations between the parameters and signal. (Assumes excitation signal shape is [39, lenSig]).
+
+    Parameters:
+        * lenSig (type: numerical int) - the length of the excitation signal to be analysed. (Optional; Default = 1)
+
+    Returns:
+        * tq (type: numpy array of floats) - the threshold of quiet for each centre frequency of the ERBs
+        * G (type: numpy array of floats) - the frequency dependent gain values for each centre frequency of the ERBs
+        * A (type: numpy array of floats) - the frequency dependent A values for each centre frequency of the ERBs
+        * alpha (type: numpy array of floats) - the frequency dependent alpha values for each centre frequency of the ERBs
+    """
+
+    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 parameters
+
+    return tq, G, A, alpha
+
 def get_specific_loudness_low(exc_low, G, tq, A, alpha):
+    """
+    Returns the specific loudness of the low level parts of the signal. Use get_specific_loudness_parameters() and specify lenSig
+    to obtain the correct values for G, tq, A and alpha. Use get_less_than_tq() to find the indexes of the samples with lower level
+    excitations.
+
+    e.g.,
+
+    # exc_i is the excitation intensity
+    lenSig = np.shape(exc_i)[-1] # find the length of the signal
+    tq, G, A, alpha = get_specific_loudness_parameters(lenSig) # get the shaped loudness parameters
+    less_than_tq = get_less_than_tq(exc_i, tq) # find which samples are lower level
+    specific_loudness[less_than_tq] = get_specific_loudness_low(exc_low[less_than_tq], G[less_than_tq], tq[less_than_tq], A[less_than_tq], alpha[less_than_tq])
+    # only process the low level part of the signal
+
+    Parameters:
+        * exc_low (type: array-like matrix of floats) - the lower level (less than or equal to tq) excitation pattern for each ERB
+            (Required)
+        * G (type: array-like matrix of floats) - the frequency dependent loudness gain parameters for each ERB, must be same shape
+            as exc_low. (Required)
+        * tq (type: array-like matrix of floats) - the frequency dependent threshold of quiet parameters for each ERB. Must be same
+            shape as exc_low. (Required)
+        * A (type: array-like matrix of floats) - the frequency dependent A parameters for each ERB. Must be same shape as exc_low.
+        * alpha (type: array-like matrix of floats) - the frequency dependent alpha parameters for each ERB. Must be same shape as
+            exc_low. (Required)
+
+    Returns:
+        * specific_loudness_low (type: array-like matrix of floats) - An array with dimensions equal to the exc_low containing the
+            specific loudness of the signal at levels below the threshold of quiet.
+    """

     C = 0.047 #constant parameter
     specific_loudness_low = C * ((2*exc_low)/(exc_low+tq))**1.5 * ((G * exc_low + A)**alpha - A**alpha)
@@ -70,6 +127,38 @@
     return specific_loudness_low

 def get_specific_loudness_mid(exc_mid, G, A, alpha):
+    """
+    Returns the specific loudness of the mid level parts of the signal.
+
+    e.g.,
+
+    # exc_i is the excitation intensity
+    lenSig = np.shape(exc_i)[-1] # find the length of the signal
+    tq, G, A, alpha = get_specific_loudness_parameters(lenSig) # get the shaped loudness parameters
+    less_than_tq = get_less_than_tq(exc_i, tq) # find which samples are lower level
+    greater_than_tq = ~less_than_tq # find which samples are greater than the threshold of quiet
+    greater_than_tl = get_greater_than_tl(exc_i) # find which samples are greater than the loud threshold
+    gttq_lttl = greater_than_tq[~greater_than_tl] # find which samples are mid level
+    specific_loudness[gttq_lttl] = get_specific_loudness_low(exc_low[gttq_lttl], G[gttq_lttl], tq[gttq_lttl], A[gttq_lttl], alpha[gttq_lttl])
+    # only process the mid level part of the signal
+
+    NOTE: The above is an example of use assuming the higher level processing IS NOT IGNORED. Use variable greater_than_tq if processing
+    higher levels equivalent to the mid.
+
+
+    Parameters:
+        * exc_mid (type: array-like matrix of floats) - the mid level (larger than tq (and, optionally less than high level threshold))
+            excitation pattern for each ERB (Required)
+        * G (type: array-like matrix of floats) - the frequency dependent loudness gain parameters for each ERB, must be same shape
+            as exc_low. (Required)
+        * A (type: array-like matrix of floats) - the frequency dependent A parameters for each ERB. Must be same shape as exc_low.
+        * alpha (type: array-like matrix of floats) - the frequency dependent alpha parameters for each ERB. Must be same shape as
+            exc_low. (Required)
+
+    Returns:
+        * specific_loudness_mid (type: array-like matrix of floats) - An array with dimensions equal to the exc_mid containing the
+            specific loudness of the signal at mid levels.
+    """

     C = 0.047 #constant parameter
     specific_loudness_mid = C * ((G * exc_mid + A)**alpha - A**alpha)
@@ -78,7 +167,24 @@

 def get_specific_loudness_high(exc_high):
     """
-
+    Returns the specific loudness of the high level parts of the signal.
+
+    e.g.,
+
+    # exc_i is the excitation intensity
+    lenSig = np.shape(exc_i)[-1] # find the length of the signal
+    tq, G, A, alpha = get_specific_loudness_parameters(lenSig) # get the shaped loudness parameters
+    greater_than_tl = get_greater_than_tl(exc_i) # find which samples are greater than the loud threshold
+    specific_loudness[greater_than_tl] = get_specific_loudness_low(exc_low[greater_than_tl], G[greater_than_tl], tq[greater_than_tl], A[greater_than_tl], alpha[greater_than_tl])
+    # only process the mid level part of the signal
+
+    Parameters:
+        * exc_high (type: array-like matrix of floats) - the high level (larger than the threshold of high level) excitation pattern
+        for each ERB (Required)
+
+    Returns:
+        * specific_loudness_high (type: array-like matrix of floats) - An array with dimensions equal to the exc_high containing the
+            specific loudness of the signal at high levels.
     """

     C = 0.047 #constant parameter
@@ -94,8 +200,8 @@
         * 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
+        * 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]
@@ -150,6 +256,8 @@
         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

@@ -325,27 +433,6 @@

     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.
@@ -405,6 +492,48 @@

     return y

+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.
+
+    Parameters:
+        * fs (type: numerical) - sampling frequency of the signal to be filtered (Required)
+        * tc (type: numerical) - the time constant of the filter. (Optional; Default = 0.01)
+
+    Returns:
+        * b (type: numerical) - the coefficient of x[n]. Equal to alpha.
+        * a (type: numpy array of floats) - an array of size 2 of the coefficients of y[n] (a[0] = 1) and y[n-1]
+    """
+
+    T = 1.0/fs
+    alpha = T/(tc + T)
+    b = [alpha]
+    a = [1, -(1-alpha)]
+
+    return b, a
+
+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 fir_bandpass(low, high, fs, order = 1024, verbose = False):
     """
     A function which returns the b coefficients for a bandpass fir filter with specified requirements.