Mercurial > hg > cm

--- a/cortical_model.py	Sun Feb 23 16:58:00 2014 +0000
+++ b/cortical_model.py	Sun Feb 23 22:50:48 2014 +0000
@@ -2,22 +2,30 @@
 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
+    * Simple frequency-modulation analysis model.
+    * Glasberg and Moore model using gammatone filters

 Packaged dependencies:
-    * utils.py and/or utils.pyc
-    * erb.dat
-    * outMidFir.dat
-    * tq.dat
+    > utils.py and/or utils.pyc
+    > erb.dat
+    > outMidFir.dat
+    > tq.dat

 External dependencies:
-    * scipy
-    * numpy
-    * copy
-    * matplotlib
+    > scipy
+    > numpy
+    > copy
+    > matplotlib
+
+Primary functions:
+    > calc_loudness_noise
+    > calc_loudness_quiet
+    > get_specific_loudness_noise
+    > get_specific_loudness_quiet
+    > get_excitation_i
+    > plot_excitation_response
+    > gammatone_filter
+    > holdsworth_gammatone_filter
 """

 import sys
@@ -28,24 +36,125 @@
 from copy import deepcopy
 import matplotlib.pyplot as plt

-def get_specific_loudness_noise(exc_i, inc_loud_region = False):
+def calc_loudness_noise(input, fs, SPL, rectify=False, inc_loud_region = False, verbose = False):

+    numTracks = np.shape(input)[0]
+    lenSig = np.shape(input)[-1]
+    exc_i = np.zeros((numTracks, 39, lenSig))
+
+    for i in range(numTracks):
+        exc_i[i] = get_excitation_i(input[i], fs, SPL, rectify = rectify, verbose = verbose)
+
+    sl = get_specific_loudness_noise(exc_i, inc_loud_region = inc_loud_region)
+    loudness = 2*np.sum(sl,1)

+    return loudness

-    return
+def calc_loudness_quiet(input, fs, SPL, rectify=False, inc_loud_region = False, verbose = False):
+
+    exc_i = get_excitation_i(input, fs, SPL, rectify = rectify, verbose = verbose)
+    sl = get_specific_loudness_noise(exc_i, inc_loud_region = inc_loud_region)
+
+    loudness = 2*np.sum(sl,0)
+
+    return loudness

-def get_specific_loudness(exc_i, inc_loud_region = False):
+def get_specific_loudness_noise(exc_i, analyse = "all", inc_loud_region = False):
     """
-    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.
+
+    """
+    exc_i = np.array(exc_i)
+
+    lenSig = np.shape(exc_i)[-1] #length signal
+    numTracks = np.shape(exc_i)[0]
+    tq, G, A, alpha = get_specific_loudness_parameters(lenSig)
+    K = calc_K_parameter(lenSig)
+    specific_loudness = exc_low = exc_mid = exc_high = np.zeros(np.shape(exc_i))
+
+    exc_signoise = sum(exc_i, 0)
+    signoise_gt_tl = get_greater_than_tl(exc_signoise) #boolean array of elements greater than the loud threshold
+    signoise_lt_tl = ~signoise_gt_tl
+
+    if(analyse.lower() == "all"):
+        loop = range(numTracks)
+    elif(analyse.lower() == "first")
+        loop = (0,)
+
+    for i in loop:
+
+        exc_sig = exc_i[i]
+        exc_noise = exc_signoise - exc_sig
+        tn = K*exc_noise + tq
+
+        less_than_tq = get_less_than_tq(exc_sig, tn) #boolean array of elements less than or equal to the threshold of quiet
+        greater_than_tq = ~less_than_tq #boolean array of elements greater than the threshold of quiet
+
+        if inc_loud_region:
+
+            c1 = signoise_lt_tl * greater_than_tn
+            c2 = signoise_lt_tl * less_than_tn
+            c3 = signoise_gt_tl * greater_than_tn
+            c4 = signoise_gt_tl * less_than_tn
+            specific_loudness[c3] = get_specific_loudness_noise_c3(exc_sig[c3], exc_noise[c3], tq[c3], tn[c3], G[c3], A[c3], alpha[c3])
+            specific_loudness[c4] = get_specific_loudness_noise_c4(exc_sig[c4], exc_noise[c4], tq[c4], tn[c4], G[c4], A[c4], alpha[c4])
+
+        else:
+            c1 = greater_than_tq
+            c2 = less_than_tq
+
+        specific_loudness[i, c1] = get_specific_loudness_noise_c1(exc_sig[c1], exc_noise[c1], tq[c1], tn[c1], G[c1], A[c1], alpha[c1])
+        specific_loudness[i, c2] = get_specific_loudness_noise_c2(exc_sig[c2], exc_noise[c2], tq[c2], tn[c2], G[c2], A[c2], alpha[c2])
+
+    return specific_loudness
+
+def get_specific_loudness_noise_c1(exc_sig, exc_noise, tq, tn, G, A, alpha):
+    """
+
+    """
+    C = 0.046871 #constant parameter
+    sl_c1 = C*(((exc_sig + exc_noise)*G + A)**alpha - A**alpha) - C*((((tn + exc_noise)*G + A)**alpha) - (tq*G + A)**alpha) * (tn/exc_sig)**0.3
+
+    return sl_c1
+
+def get_specific_loudness_noise_c2(exc_sig, exc_noise, tq, tn, G, A, alpha):
+    """
+    Parameters:
+            exc_sig
+            exc_noise
+            tq
+            tn
+            G
+            A
+            alpha
+    """
+
+    C = 0.046871 #constant parameter
+    sl_c2 = (C*((2*exc_sig)/(exc_sig + tn))**1.5) * (((tq*G + A)**alpha - A**alpha)/(((tn + exc_noise)*G + A)**alpha - (exc_noise*G + A)**alpha)) * (((exc_sig + exc_noise)*G + A)**alpha - (exc_noise*G + A)**alpha)
+
+    return sl_c2
+
+def get_specific_loudness_noise_c3(exc_sig, exc_noise, tq, tn, G, A, alpha):
+
+    C = 0.046871 #constant parameter
+    C2 = C/((1.04 * 10**6)**0.5)
+    sl_c3 = (C2*(exc_sig + exc_noise)**0.5) - C2*((tn + exc_noise)**0.5 - (tq*G + A)**alpha + A**alpha) * (tn/exc_sig)**0.3
+
+    return sl_c3
+
+def get_specific_loudness_noise_c4(exc_sig, exc_noise, tq, tn, G, A, alpha):
+
+    C = 0.046871 #constant parameter
+    sl_c4 = (C*((2*exc_sig)/(exc_sig + tn))**1.5) * (((tq*G + A)**alpha - A**alpha)/((tn + exc_noise)**0.5 - (exc_noise)**0.5)) * ((exc_sig + exc_noise)**0.5 - (exc_noise)**0.5)
+
+    return sl_c4
+
+def get_specific_loudness_quiet(exc_i, inc_loud_region = False):
+    """
+    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)
+        * 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)

@@ -74,8 +183,7 @@

 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]).
+    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)
@@ -97,13 +205,30 @@

     return tq, G, A, alpha

+def calc_K_parameter(lenSig = 1):
+    """
+
+    """
+
+
+    nerbs = 39
+    erbs = np.array(range(1,nerbs+1))
+    fc = (10**(erbs/21.366)-1)/0.004368
+
+    fkRef = 10**np.append(np.append(np.array(1), np.linspace(np.log10(50), np.log10(500), 11)), np.linspace(3,4.5,16))
+    kRef = np.append(np.append(13, np.linspace(13,-3,11)), np.ones((1,16))*-3)
+    K = np.interp(fc,fkRef,kRef)
+    K = np.transpose(np.tile(K, (lenSig,1)))
+
+    return K
+
 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.,
+    EXAMPLE USE

     # exc_i is the excitation intensity
     lenSig = np.shape(exc_i)[-1] # find the length of the signal
@@ -137,7 +262,7 @@
     """
     Returns the specific loudness of the mid level parts of the signal.

-    e.g.,
+    EXAMPLE USE

     # exc_i is the excitation intensity
     lenSig = np.shape(exc_i)[-1] # find the length of the signal
@@ -149,25 +274,21 @@
     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.
+    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)
+        * 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)
+        * 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.
+        * 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
+    C = 0.046871 #constant parameter
     specific_loudness_mid = C * ((G * exc_mid + A)**alpha - (A**alpha))

     return specific_loudness_mid
@@ -176,7 +297,7 @@
     """
     Returns the specific loudness of the high level parts of the signal.

-    e.g.,
+    EXAMPLE USE

     # exc_i is the excitation intensity
     lenSig = np.shape(exc_i)[-1] # find the length of the signal
@@ -190,11 +311,10 @@
         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.
+        * 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
+    C = 0.046871 #constant parameter
     specific_loudness_high = C * (exc_high / (1.04 * 10**6))**0.5

     return specific_loudness_high
@@ -207,8 +327,7 @@
         * 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]
@@ -225,8 +344,7 @@
         * tq (type: array-like matrix of floats) - the threshold of quiet for each ERB.

     Returns:
-        * le_tq (type: array-like matrix of booleans) - a boolean array with dimensions equal to the input
-            specifying if the excitation intensity is less than the threshold of quiet
+        * le_tq (type: array-like matrix of booleans) - a boolean array with dimensions equal to the input specifying if the excitation intensity is less than the threshold of quiet
     """

     if (np.shape(exc_i)!=np.shape(tq)):
@@ -243,16 +361,11 @@
     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)
+        * 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].
+        * 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)
@@ -261,7 +374,7 @@

     inputOMFir = outMidFir(inputPa)

-    b = gamma_tone_filter(fs, normalise = True)
+    b = gammatone_filter(fs, normalise = True)
     gtfs = decomposition(inputOMFir, b, verbose = verbose)

     if (rectify):
@@ -295,9 +408,9 @@

     if(outMidFilt): input = outMidFir(input)
     if(gammatone):
-        b = gamma_tone_filter(fs, normalise = True)
+        b = gammatone_filter(fs, normalise = True)
         input = decomposition(input, b)
-        #input = holdsworthGamma(input,fs)
+        #input = holdsworth_gammatone_filter(input,fs)
         numPlot = range(np.shape(input)[0])
     else:
         numPlot = (0,)
@@ -312,8 +425,7 @@

 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.
+    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.
@@ -321,9 +433,7 @@
         * 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].
+        * 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)
@@ -344,8 +454,7 @@
         * 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.
+        * 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)
@@ -354,17 +463,16 @@

     return y

-def gamma_tone_filter(fs, normalise = False):
+def gammatone_filter(fs, normalise = False):
     """
-    A function to filter to decompose the input signal into 39 different gammatones filtered signals modelling the ERBs.
+    A function to return the filter coefficients of 39 different gammatones with fc of ERBs 1 to 39.

     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].
+        * 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].
     """

     nerbs = 39
@@ -390,7 +498,7 @@

     return b

-def holdsworthGamma(input, fs):
+def holdsworth_gammatone_filter(input, fs):
     """

     """
@@ -432,6 +540,8 @@
     Returns:
         * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
             as a pressure signal.
+
+    NOTE If input signal contains one or more zeros, a divide by zero warning will be outputted. This warning can be ignored.
     """

     input = np.array(input)
@@ -547,14 +657,16 @@
     shape = shape[0:]
     y = np.zeros(shape)

-    if(verbose): sys.stdout.write("Running decomposition.\n")
+    if(verbose): sys.stdout.write("Running filterbank decomposition.\n")
     for i in range(nFilters):
         if(verbose):
             sys.stdout.write("\r" + str(int(np.round(100.0*i/nFilters))) + "% complete.")
             sys.stdout.flush()
         x = deepcopy(input)
         y[i] = sp.lfilter(b[i],a[i],x)
-
+
+    if(verbose): sys.stdout.write("\n")
+
     return y

 def get_Modulation_filter_bank_coefficients(fs, fc, bw, order = 1024, verbose = False):
@@ -588,7 +700,6 @@
             sys.stdout.flush()
         low = fc[i]-bw[i]/2;
         high = fc[i]+bw[i]/2;
-        if(verbose): sys.stdout.write("\r Low: " + str(low) + "; High: " + str(high))
         b = fir_bandpass(low, high, fs, order, verbose)

     return b
@@ -644,7 +755,7 @@
         * 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)
+        * verbose (type: boolean) - determines whether to display info on the current subroutine
         of the procedure. (Optional; Default = False)

     Returns:
--- a/run.py	Sun Feb 23 16:58:00 2014 +0000
+++ b/run.py	Sun Feb 23 22:50:48 2014 +0000
@@ -7,8 +7,10 @@
 sine1 = np.sin(2*np.pi*t*100/44100)
 x = np.append(sine, sine1)

-exc_i = cm.get_excitation_i(x, 44100, 60, verbose = True)
-sl = cm.get_specific_loudness(exc_i)
+exc_i = cm.get_excitation_i(x, 44100, 90, verbose = True)
+exc_i = np.reshape(exc_i, (1,)+np.shape(exc_i))
+sl = cm.get_specific_loudness_noise(exc_i)
+sl = np.reshape(sl, (39,88200))

 #for i in range(39):
 #    utils.plot_waveform(sl[i])