view testELC.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.

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