comparison cortical_model.py @ 8:d5693c7aa0ae

* Excitation pattern extraction is much faster. Calculating partial loudness for three 2 second components is reduced from ~180s to ~25s.
* runScript1.py is a script to get loudness data of all Stimuli within a set data tree. Not useful for anyone else.
* testELC.py calculates the equal loudness contours. Can also be used to get the equal loudness of two sounds. When I get the time, this will be added to a loudness package.
* loudnessRatios performs some descriptive statistics on loudness data obtained from the data within the set file tree. Probably not useful for anybody else but the author.

     > tq.dat
     
 External dependencies:
     > scipy
     > numpy
-    > copy
     > matplotlib
     
 Primary functions:
     > calc_partial_loudness
     > calc_loudness

 import sys
 import utils
 import scipy.signal as sp
 import scipy.fftpack as fft
 import numpy as np
-from copy import deepcopy
 import matplotlib.pyplot as plt
 
 def calc_loudness_noise(input, fs, SPL, rectify=False, inc_loud_region = False, verbose = False):
 
     numTracks = np.shape(input)[0]

         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 sl, loudness
 
 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_quiet(exc_i, inc_loud_region = inc_loud_region)
     
     loudness = 2*np.sum(sl,0)
     
-    return loudness
+    return sl, loudness
 
 def get_specific_loudness_noise(exc_i, analyse = "all", inc_loud_region = False):
     """
         
     """

         
         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
+        less_than_tn = get_less_than_tq(exc_sig, tn) #boolean array of elements less than or equal to the threshold of quiet
+        greater_than_tn = ~less_than_tn #boolean array of elements greater than the threshold of quiet
         
         if inc_loud_region:
-            
             c1 = signoise_lt_tl * greater_than_tn #case 1
             c2 = signoise_lt_tl * less_than_tn #case 2
             c3 = signoise_gt_tl * greater_than_tn #case 3
-            c4 = signoise_gt_tl * less_than_tn
+            c4 = signoise_gt_tl * less_than_tn #case 4
             specific_loudness[i, c3] = get_specific_loudness_noise_c3(exc_sig[c3], exc_noise[c3], tq[c3], tn[c3], G[c3], A[c3], alpha[c3])
             specific_loudness[i, 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 #case 1
-            c2 = less_than_tq #case 2
+            c1 = greater_than_tn #case 1
+            c2 = less_than_tn #case 2
 
         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

     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
     greater_than_tq = ~less_than_tq #boolean array of elements greater than the threshold of quiet
     greater_than_tl = get_greater_than_tl(exc_i) #boolean array of elements greater than the loud threshold
-    gttq_lttl = greater_than_tq[~greater_than_tl] #boolean array of elements greater than the threshold of quiet but less than the threshold of loud
+    gttq_lttl = greater_than_tq * ~greater_than_tl #boolean array of elements greater than the threshold of quiet but less than the threshold of loud
     
     exc_low = exc_i[less_than_tq]
     specific_loudness[less_than_tq] = get_specific_loudness_low(exc_low, G[less_than_tq], tq[less_than_tq], A[less_than_tq], alpha[less_than_tq])
     
     if(inc_loud_region):

     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 = 10**(K/10)
     K = np.transpose(np.tile(K, (lenSig,1)))
 
     return K
 
 def get_specific_loudness_low(exc_low, G, tq, A, alpha):

     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)
+        x = np.array(input)
+        y[i] = sp.fftconvolve(x,b[i],'same')
 
     if(verbose): sys.stdout.write("\n")
 
     return y