changeset 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.

--- a/cortical_model.py	Mon Feb 24 15:31:15 2014 +0000
+++ b/cortical_model.py	Mon Mar 31 19:22:23 2014 +0100
@@ -14,7 +14,6 @@
 External dependencies:
     > scipy
     > numpy
-    > copy
     > matplotlib
     
 Primary functions:
@@ -33,7 +32,6 @@
 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):
@@ -48,7 +46,7 @@
     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):
     
@@ -57,7 +55,7 @@
     
     loudness = 2*np.sum(sl,0)
     
-    return loudness
+    return sl, loudness
 
 def get_specific_loudness_noise(exc_i, analyse = "all", inc_loud_region = False):
     """
@@ -86,21 +84,19 @@
         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])
@@ -165,7 +161,7 @@
     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])
@@ -218,6 +214,7 @@
     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
@@ -658,8 +655,8 @@
         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")
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/loudnessRatios.py	Mon Mar 31 19:22:23 2014 +0100
@@ -0,0 +1,147 @@
+import cortical_model as cm
+import utils
+import numpy as np
+from os import chdir
+from os import listdir
+import matplotlib.pyplot as plt
+
+chdir("Stimuli")
+dirList = listdir(".")
+
+if(dirList[0] == ".DS_Store"):
+    dirList.pop(0)
+
+numEntries = len(dirList)
+ratio = np.zeros((3, numEntries))
+jazzRatio = np.zeros((3, 6))
+rockRatio = np.zeros((3, 5))
+funkRatio = np.zeros((3, 6))
+
+i = 0
+jazzi = 0
+funki = 0
+rocki = 0
+
+#get ratios
+for dir in dirList:
+    chdir(dir)
+    loudness = np.load("loudness.npy")
+    maxHat = loudness[0].max()
+    maxKick = loudness[1].max()
+    maxSnare = loudness[2].max()
+
+    ratio[0][i] = maxHat/maxKick
+    ratio[1][i] = maxHat/maxSnare
+    ratio[2][i] = maxKick/maxSnare
+    
+    if(dir.find("jazz") != -1):
+        jazzRatio[0][jazzi] = ratio[0][i]
+        jazzRatio[1][jazzi] = ratio[1][i]
+        jazzRatio[2][jazzi] = ratio[2][i]
+        jazzi = jazzi + 1
+    elif(dir.find("funk") != -1):
+        funkRatio[0][funki] = ratio[0][i]
+        funkRatio[1][funki] = ratio[1][i]
+        funkRatio[2][funki] = ratio[2][i]
+        funki = funki + 1
+    elif(dir.find("rock") != -1):
+        rockRatio[0][rocki] = ratio[0][i]
+        rockRatio[1][rocki] = ratio[1][i]
+        rockRatio[2][rocki] = ratio[2][i]
+        rocki = rocki + 1
+
+    chdir("../")
+
+    i=i+1
+
+chdir("..")
+
+def plot_all():
+    #calculate descriptive statistics
+    meanRatio = np.mean(ratio,1)
+    stdRatio = np.std(ratio,1)
+    seRatio = stdRatio/np.sqrt(numEntries)
+    rangeRatio = ratio.max(1)-ratio.min(1)
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+    ax.plot(meanRatio, 'xk')
+    ax.plot(meanRatio+stdRatio, '_b')
+    ax.plot(meanRatio-stdRatio, '_b')
+    ax.plot(meanRatio+seRatio, '_r')
+    ax.plot(meanRatio-seRatio, '_r')
+    #ax.plot(ratio, '.k')
+    plt.xlim([-1,3])
+    plt.xticks(np.array(range(5))-1,['','Hat:Kick','Hat:Snare','Kick:Snare',''])
+    plt.xlabel('Stimuli')
+    plt.ylabel('Loudness Ratio')
+
+    plt.show()
+
+def plot_indiv():
+    
+    #calculate descriptive statistics
+    meanJazz = np.mean(jazzRatio,1)
+    meanRock = np.mean(rockRatio,1)
+    meanFunk = np.mean(funkRatio,1)
+    stdJazz = np.std(jazzRatio,1)
+    stdRock = np.std(jazzRatio,1)
+    stdFunk = np.std(jazzRatio,1)
+    seJazz = stdJazz/np.sqrt(6)
+    seRock = stdRock/np.sqrt(5)
+    seFunk = stdFunk/np.sqrt(6)
+    #rangeRatio = ratio.max(1)-ratio.min(1)
+    
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+    ax.plot(meanJazz, 'or')
+    ax.plot(meanRock, 'og')
+    ax.plot(meanFunk, 'ob')
+    ax.legend(('Jazz', 'Rock', 'Funk'))
+    ax.plot(meanJazz+stdJazz, '_r')
+    ax.plot(meanRock+stdRock, '_g')
+    ax.plot(meanFunk+stdFunk, '_b')
+    ax.plot(meanJazz-stdJazz, '_r')
+    ax.plot(meanRock-stdRock, '_g')
+    ax.plot(meanFunk-stdFunk, '_b')
+    ax.plot(meanJazz+seJazz, '+r')
+    ax.plot(meanRock+seRock, '+g')
+    ax.plot(meanFunk+seFunk, '+b')
+    ax.plot(meanJazz-seJazz, '+r')
+    ax.plot(meanRock-seRock, '+g')
+    ax.plot(meanFunk-seFunk, '+b')
+    #ax.plot(ratio, '.k')
+    plt.xlim([-1,3])
+    plt.xticks(np.array(range(5))-1,['','Hat:Kick','Hat:Snare','Kick:Snare',''])
+    plt.xlabel('Stimuli')
+    plt.ylabel('Loudness Ratio')
+    
+    plt.show()
+
+plot_indiv()
+del dirList, maxHat, maxKick, maxSnare
+
+#plt.show()
+
+#chdir("../Stimuli2")
+#
+#for dir in dirList:
+#    chdir(dir)
+#    loudness = np.load("loudness.npy")
+#    maxOH = loudness[0].max()
+#    maxKick = loudness[1].max()
+#    maxSnare = loudness[2].max()
+#
+#    ratio = np.zeros((3))
+#    ratio[0] = maxOH/maxKick
+#    ratio[1] = maxOH/maxSnare
+#    ratio[2] = maxKick/maxSnare
+#
+#    if(ratio[0] > 16):
+#        print dir
+#
+#    np.save("ratio", ratio)
+#
+#    plt.plot(ratio, 'or')
+#
+#    chdir("../")
\ No newline at end of file

--- a/run.py	Mon Feb 24 15:31:15 2014 +0000
+++ b/run.py	Mon Mar 31 19:22:23 2014 +0100
@@ -2,17 +2,16 @@
 import numpy as np
 import utils
 
-t = np.array(range(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
+fs,hat = utils.wavread2('hat.wav')
+hat = sum(hat,0)
+fs,kick = utils.wavread2('kick.wav')
+kick = sum(kick,0)
+fs,snare = utils.wavread2('snare.wav')
+snare = sum(snare,0)
 
-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)
+drums = np.zeros((3,len(hat)))
+drums[0]=hat
+drums[1]=kick
+drums[2]=snare
 
-loudness = cm.calc_loudness_quiet(x,44100,90)
-
-utils.plot_waveform(loudness)
\ No newline at end of file
+loudness = cm.calc_loudness_noise(drums,fs,90, inc_loud_region = True)
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/runScript1.py	Mon Mar 31 19:22:23 2014 +0100
@@ -0,0 +1,62 @@
+import os
+import utils
+import cortical_model as cm
+import numpy as np
+
+os.chdir("Stimuli")
+dirList = os.listdir('.')
+if dirList[0] == '.DS_Store':
+    dirList.pop(0)
+
+print "Set 1"
+
+for dir in dirList:
+    print dir
+    
+    os.chdir(dir)
+    fs,hat = utils.wavread2('hat.wav')
+    fs,kick = utils.wavread2('kick.wav')
+    fs,snare = utils.wavread2('snare.wav')
+
+    drums = np.zeros((3,len(hat)))
+    drums[0]=hat
+    drums[1]=kick
+    drums[2]=snare
+    exc = np.zeros((3,39,len(hat)))
+
+    for i in range(3):
+        exc[i] = cm.get_excitation_i(drums[i],fs,90)
+
+    np.save('exc', exc)
+
+    os.chdir('..')
+
+kick = 0
+snare = 0
+hat = 0
+
+print "Set 2"
+
+os.chdir("../Stimuli2")
+
+for dir in dirList:
+    print dir
+    
+    os.chdir(dir)
+    fs,oh = utils.wavread2('oh.wav')
+    fs,kick = utils.wavread2('kick.wav')
+    fs,snare = utils.wavread2('snare.wav')
+    
+    drums = np.zeros((3,len(oh)))
+    drums[0]=oh
+    drums[1]=kick
+    drums[2]=snare
+    
+    exc = np.zeros((3,39,len(oh)))
+    
+    for i in range(3):
+        exc[i] = cm.get_excitation_i(drums[i],fs,90)
+    
+    np.save('exc', exc)
+    
+    os.chdir('..')
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testELC.py	Mon Mar 31 19:22:23 2014 +0100
@@ -0,0 +1,91 @@
+import cortical_model as cm
+import numpy as np
+import matplotlib.pyplot as plt
+import pdb
+import scipy.optimize as opt
+import scipy.interpolate as interp
+
+fArray = [20, 50, 100, 200, 400, 800, 1000, 2000, 3000, 4000, 5000, 6000, 10000, 12000, 15000]
+refLoudness = None
+refDB = None
+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
+
+def get_elc_range(verbose=False):
+    levels = range(10, 120, 10)
+    elc = np.zeros((len(levels), len(fArray)))
+    
+    for i in range(len(levels)):
+        if(verbose):
+            print "Level: " + str(levels[i])
+        elc[i] = get_elc(levels[i], verbose=verbose)
+        plot_elc(elc[i], show=False)
+
+    plt.show()
+
+    return elc
+
+def get_elc(SPL = 90, verbose=False):
+    fs = 44100
+    t = np.array(range(fs))
+    elc = np.zeros((len(fArray)))
+    set_refLoudness(fs, SPL)
+    
+    for i in range(len(fArray)):
+        if (fArray[i] == 1000):
+            elc[i] = SPL
+        else:
+            if(verbose):
+                print fArray[i]
+            input = np.sin(fArray[i] * t * 2 * np.pi / fs) * env
+            elc[i] = get_el(input, fs, refLoudness, SPL, verbose=verbose)
+
+    return elc
+
+def get_el(input, fs, refL, SPL = 90, verbose=False):
+    
+    x = opt.fsolve(getLoudnessDifference, SPL, (input, fs, refL, verbose), factor=0.1, xtol=0.00005)
+    
+    return x
+
+def getLoudnessDifference(SPL, input, fs, refL, verbose=False):
+
+    if(verbose):
+        print "SPL: ", SPL
+    compLoudness = rms(cm.calc_loudness_quiet(input, fs, SPL, rectify=False, inc_loud_region=True)[1])
+    difference = np.abs(refL - compLoudness)
+    if(verbose):
+        print "Difference: ", difference
+
+    return difference
+
+def rms(input):
+
+    output = np.sqrt(np.mean(input**2))
+
+    return output
+
+def set_refLoudness(fs, SPL = 90):
+    refDB = SPL
+    t = np.array(range(fs))
+
+    global refLoudness
+    input = np.sin(1000 * t * 2 * np.pi / fs) * env
+    refLoudness = rms(cm.calc_loudness_quiet(input, fs, SPL, rectify=False, inc_loud_region=True)[1])
+
+    return
+
+def plot_elc(cmOut, show = True):
+
+    x = np.array(range(min(fArray),max(fArray),5))
+    values = interp.interp1d(fArray, cmOut)
+    plt.gca().set_xscale('log')
+    plt.ylim((0,170))
+    plt.plot(fArray, cmOut, 'x', x, values(x))
+    if(show):
+        plt.show()
+
+    return
\ No newline at end of file

--- a/utils.py	Mon Feb 24 15:31:15 2014 +0000
+++ b/utils.py	Mon Mar 31 19:22:23 2014 +0100
@@ -13,9 +13,21 @@
 """
 
 import numpy as np
-import scipy.io.wavfile as wave
+import scipy.io.wavfile as wav
 import matplotlib.pyplot as plt
 import scipy.fftpack as fft
+import wave
+import struct
+from os import getcwd
+
+gPath = getcwd()
+
+def set_path_to_data(path):
+
+    global gPath
+    gPath = path
+
+    return
 
 def load_erb_data():
 	"""
@@ -34,7 +46,7 @@
         * erb.dat (place in the same folder as utils.pyc)
 	"""
 	# read data from a text file
-	data=np.loadtxt("erb.dat",delimiter=",")
+	data=np.loadtxt(gPath + "/erb.dat",delimiter=",")
 
 	# load centre frequencies from the first column into fc
 	fc=np.array(data[:,0])
@@ -58,7 +70,7 @@
         * outMidFir.dat (place in the same folder as utils.pyc)
     """
 
-    b=np.array(np.loadtxt("outMidFir.dat", delimiter=","))
+    b=np.array(np.loadtxt(gPath + "/outMidFir.dat", delimiter=","))
 
     return b
 
@@ -101,7 +113,7 @@
         * tq.dat (place in the same folder as utils.pyc)
     """
 
-    tq = np.array(np.loadtxt("tq_dB.dat",delimiter=","))
+    tq = np.array(np.loadtxt(gPath + "/tq_dB.dat",delimiter=","))
 
     return tq
 
@@ -119,7 +131,7 @@
         * A.dat (place in the same folder as utils.pyc)
     """
     
-    A = np.array(np.loadtxt("A.dat",delimiter=","))
+    A = np.array(np.loadtxt(gPath + "/A.dat",delimiter=","))
     
     return A
 
@@ -137,7 +149,7 @@
         * alpha.dat (place in the same folder as utils.pyc)
     """
     
-    alpha = np.array(np.loadtxt("alpha.dat",delimiter=","))
+    alpha = np.array(np.loadtxt(gPath + "/alpha.dat",delimiter=","))
     
     return alpha
 
@@ -176,11 +188,58 @@
             to an amplitude range of -1 to 1.
     """
 
-    fs,data = wave.read(file)
+    fs,data = wav.read(file)
     data = np.array(int_to_nfloat(data))
 
     return fs,data
 
+def wavread2(filename):
+
+    file = wave.open(filename)
+    nChans, sampwidth, fs, nFrames, comptype, compname = file.getparams()
+    
+    N = 8*sampwidth
+    
+    data = np.zeros((nChans, nFrames))
+    
+    if sampwidth == 2:
+        unpackChar = 'h'
+    elif sampwidth == 3:
+        unpackChar = 'hc'
+    elif sampwidth == 4:
+        unpackChar = 'i'
+
+    frame = [[] for i in range(nChans)]
+
+    
+    for i in range(nFrames):
+        bytes = file.readframes(1)
+        
+        if nChans > 1:
+            for j in range(nChans):
+                index = j*sampwidth
+                frame[j] = bytes[index:index+sampwidth]
+        else:
+            frame[0] = bytes
+
+        for j in range(nChans):
+            if sampwidth == 2:
+                data[j][i] = struct.unpack('<h', frame[j])[0]
+            elif sampwidth == 3:
+                data[j][i] = struct.unpack('<i', frame[j] + ('\0' if frame[j][2] < '\x80' else '\xff'))[0]
+            elif sampwidth == 4:
+                data[j][i] = struct.unpack('<i', frame[j])[0]
+
+    data[data<0] = data[data<0]/(2**N)
+    data[data>0] = data[data>0]/((2**N)-1)
+
+    if nChans == 1:
+        data = data.reshape((nFrames))
+
+    file.close()
+
+    return fs,data
+
 def wavwrite(file, fs, data, precision = 16):
     """
     Unnormalises the audio data to a specified integer precision and converts to an int, then 
@@ -213,7 +272,7 @@
         return
 
     data = nfloat_to_int(data)
-    wave.write(file, fs, data)
+    wav.write(file, fs, data)
     
     return