Mercurial > hg > cm
changeset 5:d7b2784ff5a3
* Fixed some miscalculations from the way python does float division
* Made data files more precise
| author | Carl Bussey <c.bussey@se10.qmul.ac.uk> |
|---|---|
| date | Sun, 23 Feb 2014 16:58:00 +0000 |
| parents | 29cd3e735c4c |
| children | 546dfcd45281 |
| files | A.dat alpha.dat cortical_model.py run.py tq.dat tq_dB.dat utils.py |
| diffstat | 7 files changed, 57 insertions(+), 34 deletions(-) [+] |
line wrap: on
line diff
--- a/A.dat Sat Feb 22 13:53:55 2014 +0000 +++ b/A.dat Sun Feb 23 16:58:00 2014 +0000 @@ -1,1 +1,1 @@ -4.7209,4.7204,4.7169,4.7084,4.6964,4.6814,4.6658,4.6565,4.6494,4.6421,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387,4.6387 +4.72091216639346900052,4.72041046637031325162,4.71687810730626111422,4.70836341627149046474,4.69637854656524211805,4.68144965624411746319,4.66581035715129210928,4.65653474160862845821,4.64937797923400619737,4.64213180884578857643,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689,4.63869286742936459689
--- a/alpha.dat Sat Feb 22 13:53:55 2014 +0000 +++ b/alpha.dat Sun Feb 23 16:58:00 2014 +0000 @@ -1,1 +1,1 @@ -0.2,0.19999,0.19992,0.19975,0.19952,0.19921,0.19886,0.19863,0.19843,0.19818,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803,0.19803 +0.19999907005650591207,0.19998931739710795741,0.19992068067184415314,0.19975495561825465374,0.19951866936695022980,0.19921267556305100044,0.19886379821729613382,0.19863097123971321101,0.19842704941898242632,0.19817953122456480330,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337,0.19803309989228007337
--- a/cortical_model.py Sat Feb 22 13:53:55 2014 +0000 +++ b/cortical_model.py Sun Feb 23 16:58:00 2014 +0000 @@ -20,6 +20,7 @@ * matplotlib """ +import sys import utils import scipy.signal as sp import scipy.fftpack as fft @@ -27,6 +28,12 @@ from copy import deepcopy import matplotlib.pyplot as plt +def get_specific_loudness_noise(exc_i, inc_loud_region = False): + + + + return + def get_specific_loudness(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 @@ -121,8 +128,8 @@ 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) + C = 0.046871 #constant parameter + specific_loudness_low = C * ((2.0*exc_low)/(0.0+exc_low+tq))**1.5 * ((G * exc_low + A)**alpha - A**alpha) return specific_loudness_low @@ -161,7 +168,7 @@ """ C = 0.047 #constant parameter - specific_loudness_mid = C * ((G * exc_mid + A)**alpha - A**alpha) + specific_loudness_mid = C * ((G * exc_mid + A)**alpha - (A**alpha)) return specific_loudness_mid @@ -249,18 +256,27 @@ """ input = np.array(input) - inputOMFir = outMidFir(input) - inputPa = v_Pascal(inputOMFir, SPL) - b = gamma_tone_filter(fs) - gtfs = decomposition(inputPa, b, verbose = verbose) + + inputPa = v_Pascal(input, SPL) + + inputOMFir = outMidFir(inputPa) + + b = gamma_tone_filter(fs, normalise = True) + gtfs = decomposition(inputOMFir, b, verbose = verbose) + if (rectify): gtfs = half_rectification(gtfs) - gtfs = pa_i(gtfs) - gtfs = at_normalise(gtfs) - b,a = exp_smoothing(fs) - gtfs = sp.lfilter(b,a,gtfs) + + gtfsI = pa_i(gtfs) + gtfsN = at_normalise(gtfsI) - return gtfs + b1,a1 = exp_smoothing(fs) + gtfsF = np.zeros(np.shape(gtfsI)) + + for i in range(39): + gtfsF[i] = sp.lfilter(b1,a1,gtfsN[i]) + + return gtfsF def plot_excitation_response(input = None, fs = 44100, outMidFilt = True, gammatone = True, xscale = 'log', yscale = 'log'): """ @@ -279,7 +295,7 @@ if(outMidFilt): input = outMidFir(input) if(gammatone): - b = gamma_tone_filter(fs) + b = gamma_tone_filter(fs, normalise = True) input = decomposition(input, b) #input = holdsworthGamma(input,fs) numPlot = range(np.shape(input)[0]) @@ -338,7 +354,7 @@ return y -def gamma_tone_filter(fs): +def gamma_tone_filter(fs, normalise = False): """ A function to filter to decompose the input signal into 39 different gammatones filtered signals modelling the ERBs. @@ -368,6 +384,10 @@ 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) + if normalise: + fftbmax = (np.abs(fft.fft(b, np.shape(b)[-1]))).max(); + b = b/fftbmax; + return b def holdsworthGamma(input, fs): @@ -415,7 +435,7 @@ """ input = np.array(input) - y = np.sign(input)*(0.00002*10**(np.log10(np.abs(input))+SPL/20)) + y = np.sign(input)*(0.00002*10**(np.log10(np.abs(input))+SPL/20.0)) return y @@ -450,7 +470,7 @@ """ input = np.array(input) - y = input / 1*(10**12) + y = input / (1*(10**-12)) return y @@ -527,9 +547,11 @@ shape = shape[0:] y = np.zeros(shape) - if(verbose): print "Running decomposition." + if(verbose): sys.stdout.write("Running decomposition.\n") for i in range(nFilters): - if(verbose): print str(100.0*i/nFilters) + "% complete." + 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) @@ -558,13 +580,15 @@ nFreqs = len(fc) - if(verbose): print "Running frequency decomposition." + if(verbose): sys.stdout.write("Running frequency decomposition.") for i in range(nFreqs): - if(verbose): print str(100.0*i/nFreqs) + "% complete." + if(verbose): + sys.stdout.write("\r" + str(int(np.round(100.0*i/nFreqs))) + "% complete.") + sys.stdout.flush() low = fc[i]-bw[i]/2; high = fc[i]+bw[i]/2; - if(verbose): print "Low: " + str(low) + "; High: " + str(high) + if(verbose): sys.stdout.write("\r Low: " + str(low) + "; High: " + str(high)) b = fir_bandpass(low, high, fs, order, verbose) return b @@ -630,7 +654,7 @@ nyquist = 0.5*fs lowNorm = low/nyquist highNorm = high/nyquist - if(verbose): print "Low: " + str(lowNorm) + "; High: " + str(highNorm) + if(verbose): sys.stdout.write("Low: " + str(lowNorm) + "; High: " + str(highNorm)) b = sp.firwin(order+1, [lowNorm, highNorm], pass_zero=False) return b \ No newline at end of file
--- a/run.py Sat Feb 22 13:53:55 2014 +0000 +++ b/run.py Sun Feb 23 16:58:00 2014 +0000 @@ -7,14 +7,13 @@ sine1 = np.sin(2*np.pi*t*100/44100) x = np.append(sine, sine1) -exc_i = cm.get_excitation_i(x, 44100, 90, verbose = True) +exc_i = cm.get_excitation_i(x, 44100, 60, verbose = True) sl = cm.get_specific_loudness(exc_i) -for i in range(39): - print "ERB " + str(i) - utils.plot_waveform(sl[i]) +#for i in range(39): +# utils.plot_waveform(sl[i]) print "Loudness" -loudness = np.sum(sl,0) +loudness = 2*np.sum(sl,0) utils.plot_waveform(loudness) \ No newline at end of file
--- a/tq.dat Sat Feb 22 13:53:55 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1 +0,0 @@ -37.038,26.436,17.727,12.824,9.8853,7.7397,6.1369,5.3208,4.712,4.0702,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73,3.73
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tq_dB.dat Sun Feb 23 16:58:00 2014 +0000 @@ -0,0 +1,1 @@ +37.03830244398031368291,26.43583190974692698205,17.72694330081553815148,12.82375753482408597961,9.88527548511461162661,7.73969571489240770745,6.13691907262802072154,5.32077857252018304735,4.71196196585856785788,4.07020279955640518210,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224,3.72999999999999998224
--- a/utils.py Sat Feb 22 13:53:55 2014 +0000 +++ b/utils.py Sun Feb 23 16:58:00 2014 +0000 @@ -80,13 +80,13 @@ * A.dat (place in the same folder as utils.pyc) * alpha.dat (place in the same folder as utils.pyc) """ - tq = load_tq_data() + tq_dB = load_tqdB_data() A = load_A_data() alpha = load_alpha_data() - return tq, A, alpha + return tq_dB, A, alpha -def load_tq_data(): +def load_tqdB_data(): """ Loads and returns the excitation threshold of quiet for each ERB fc. @@ -101,7 +101,7 @@ * tq.dat (place in the same folder as utils.pyc) """ - tq = np.array(np.loadtxt("tq.dat",delimiter=",")) + tq = np.array(np.loadtxt("tq_dB.dat",delimiter=",")) return tq