Mercurial > hg > cm

changeset 7:0ee61aeb789a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
* Changed normalisation of gammatone_filter * Changed run.py to test calc_loudness_quiet()
author Carl Bussey <c.bussey@se10.qmul.ac.uk>
date Mon, 24 Feb 2014 15:31:15 +0000
parents 546dfcd45281
children d5693c7aa0ae
files cortical_model.py run.py
diffstat 2 files changed, 32 insertions(+), 39 deletions(-) [+]
line wrap: on
line diff
--- a/cortical_model.py	Sun Feb 23 22:50:48 2014 +0000
+++ b/cortical_model.py	Mon Feb 24 15:31:15 2014 +0000
@@ -18,8 +18,8 @@
     > matplotlib
     
 Primary functions:
-    > calc_loudness_noise
-    > calc_loudness_quiet
+    > calc_partial_loudness
+    > calc_loudness
     > get_specific_loudness_noise
     > get_specific_loudness_quiet
     > get_excitation_i
@@ -50,10 +50,10 @@
 
     return loudness
 
-def calc_loudness_quiet(input, fs, SPL, rectify=False, inc_loud_region = False, verbose = False):
+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)
+    sl = get_specific_loudness_quiet(exc_i, inc_loud_region = inc_loud_region)
     
     loudness = 2*np.sum(sl,0)
     
@@ -75,10 +75,10 @@
     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,)
+    if analyse.lower() == "all":
+        loop = range(numTracks);
+    elif analyse.lower() == "first":
+        loop = range(1)
     
     for i in loop:
         
@@ -91,19 +91,19 @@
         
         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
+            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
-            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])
+            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
-            c2 = less_than_tq
+            c1 = greater_than_tq #case 1
+            c2 = less_than_tq #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])
+        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
 
@@ -129,7 +129,7 @@
     """
     
     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)
+    sl_c2 = (C*((2*exc_sig)/(exc_sig + tn))**1.5) * (((tq*G + A)**alpha - A**alpha)/(((0.0 + 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
 
@@ -482,20 +482,16 @@
     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)
+        b[i] = t**(N-1) * np.exp(-2*np.pi*bw[i]*t) * np.cos(2*np.pi*fc[i]*t)
     
-    if normalise:
-        fftbmax = (np.abs(fft.fft(b, np.shape(b)[-1]))).max();
-        b = b/fftbmax;
-    
+        if normalise:
+            A = 1/(np.abs(fft.fft(b[i], 4*filterLength))).max()
+            b[i] *= A
+
     return b
 
 def holdsworth_gammatone_filter(input, fs):
--- a/run.py	Sun Feb 23 22:50:48 2014 +0000
+++ b/run.py	Mon Feb 24 15:31:15 2014 +0000
@@ -3,19 +3,16 @@
 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)
+env = np.ones((44100))
+#T = int(np.floor(0.02*44100))
+#tRamp = np.array(range(T))
+#env[tRamp] = np.sin(2*np.pi*12.5*tRamp/44100)**2
+#env[-tRamp] = np.sin(2*np.pi*12.5*tRamp/44100)**2
+
+sine = np.sin(2*np.pi*t*1000/44100) * env
+sine1 = np.sin(2*np.pi*t*100/44100) * env
 x = np.append(sine, sine1)
 
-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])
-
-print "Loudness"
-loudness = 2*np.sum(sl,0)
+loudness = cm.calc_loudness_quiet(x,44100,90)
 
 utils.plot_waveform(loudness)
\ No newline at end of file