annotate equalLoudness.py @ 11:ca98f5f26bcb tip

* More frequencies to verify toolbox with.

import pyphon as pp
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, 1500, 2000, 3000, 4000, 5000, 6000, 10000, 12000, 15000]
#dArray = np.array([5, 10, 50, 100, 150, 200, 500, 1000])/1000.0
dArray = np.array([8, 16, 32, 64, 128, 256, 512, 1024])/1000.0

def get_elc_range(xtol=0.0000005, verbose=False):
    """
        Plot the prediction of the equal loudness contours of frequencies for a range of SPL levels from 10 to 120dB, intervals of 10dB.
        
        Parameters:
            * xtol (type: float) - x tolerance of fsolve
            * verbose (type: boolean) - I tell you what I know
        
        Returns:
            * elc (type: numpy array of floats) - SPL at which each frequency in fArray is equal to i*10 phons, where i is the first dimension index
        
        """
    levels = range(20, 100, 20)
    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], xtol=xtol, verbose=verbose)
        plot_elc(elc[i], show=False)

    plt.show()

    return elc

def get_elc(SPL = 90, xtol=0.0000005, verbose=False):
    """
        Plot the prediction of the equal loudness contours for specified level.
        
        Parameters:
        * SPL (type: float) - SPL of the reference, or the loudness level
        * xtol (type: float) - x tolerance of fsolve
        * verbose (type: boolean) - I tell you what I know if you give me True
        
        Returns:
        * elc (type: numpy array of floats) - SPL at which each frequency in fArray is equal to SPL phons
        """
    fs = 44100
    t = np.array(range(44100))
    elc = np.zeros((len(fArray)))
    env = get_cos2_env(44100, fs)
    
    refSig = np.sin(1000 * t * 2 * np.pi / fs) * env
    refSig = refSig/rms(refSig)
    refLoudness = pp.calc_loudness_quiet(refSig, fs, SPL, rectify=False, inc_loud_region=True)[2].max()
    
    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
            input = input / rms(input)
            elc[i] = get_el(input, fs, refLoudness, SPL, xtol=xtol, verbose=verbose)

    return elc

def get_temporal_descrimination(verbose=False):
    """
        Return the prediction of STL and peak of STL for tones with short duration.
        
        Parameters
        * verbose (type: boolean) - I tell you what I know if you give me True
        
        Returns:
            * STLMax (type: array-like of floats) - peak short term loudness for each duration
            * STL (type: array-like of floats) - short term loudness function for each duration
        """


    fs = 44100
    STLMax = np.zeros((len(dArray)))
    STL = []
    
    for i in range(len(dArray)):
        if(verbose):
            print dArray[i]
        length = np.floor(dArray[i]*fs)
        env = get_cos2_env(length, fs, 0.005*fs)
        input = np.sin(np.arange(length)*2*np.pi*1000/fs) * env
        input = input/rms(input)
        STL.append(pp.calc_loudness_quiet(input, fs, 60, inc_loud_region=True)[2])
        STLMax[i] = STL[i].max()

    return STLMax, STL

def get_cos2_env(lengthSig, fs, lengthEnv = None):
    """
        Return a cos**2 envelope.
        
        Parameters:
            * lengthSig (type: int) - length of signal to be enveloped
            * fs - (type: int) - sampling frequency
            * lengthEnv (type: int) - length of ramps
            
        Returns:
            * env (type: numpy array of floats)- envelope
        """
    
    if(lengthEnv == None):
        lengthEnv = np.floor(0.02*fs)
    
    env = np.ones((lengthSig))
    T = int(lengthEnv)
    tRamp = np.array(range(T))
    f = 1/(lengthEnv/fs*4)
    env[tRamp] = np.sin(2*np.pi*f*tRamp/fs)**2
    env[-tRamp] = np.sin(2*np.pi*f*tRamp/fs)**2
    
    return env

def get_el(input, fs, refL, SPL = 90, xtol=0.0000005, verbose=False):
    """
        Get the SPL of the input at equal loudness to the reference.
        
        Parameters:
        * input (type: array-like of floats) - input signal
        * fs (type: int) - sampling frequency
        * SPL (type: float) - the level to start searching at
        * xtol (type: float) - x tolerance of fsolve
        * verbose (type: boolean) - I tell you what I know if you give me True
        
        Returns:
            x (type: float) - sound pressure level at equal loudness
        """
    
    x = opt.fsolve(getLoudnessDifference, SPL, (input, fs, refL, verbose), factor=0.1, xtol=xtol)
    
    return x

def getLoudnessDifference(SPL, input, fs, refL, verbose=False):
    """
        Used for get_el fsolve
        
        Parameters:
        * SPL (type: float) - the level to set input to
        * input (type: array-like of floats) - input signal
        * fs (type: int) - sampling frequency
        
        Returns:
            difference (type: float) - difference in loudness
        """
    if(verbose):
        print "SPL: ", SPL
    input = input/rms(input)
    compLoudness = pp.calc_loudness_quiet(input, fs, SPL, rectify=False, inc_loud_region=True)[2].max()
    difference = np.abs(refL - compLoudness)
    if(verbose):
        print "Difference: ", difference

    return difference

def rms(input):
    """
        RMS of the input
        
        Parameters:
        * input (type: array-like of floats): input array
        
        Returns:
        * output (type: float): the rms
        """
    output = np.sqrt(np.mean(input**2, -1))

    return output

def plot_elc(elc, show = True):
    """
        Plot an equal loudness contour
        
        Parameters:
            * elc (type: array-like of floats) - equal loudness contour x, obtained from get_elc
            * show (type: boolean) - show the figure or turn hold on
        """
    plt.gca().set_xscale('log')
    plt.ylim((0,170))
    plt.plot(fArray, elc, 'x-')
    if(show):
        plt.show()
    else:
        plt.hold(True)

    return