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

"""
A utility module to assist cortical_model.py with general procedures, e.g., file reading and writing.

Packaged dependencies:
    * erb.dat
    * outMidFir.dat
    * tq.dat
    
External dependencies:
    * scipy
    * numpy
    * matplotlib
"""

import numpy as np
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():
	"""
	Loads and returns 39 ERB frequencies and bandwidths
    
    Parameters:
        NONE
        
    Returns:
        * fc (type: numpy array of floats) - a vector of length 39 with elements containing
            the centre frequencies of each ERB.
        * bw (type: numpy array of floats) - a vector of length 39 with elements containing
            the bandwidths of each ERB.
    
    Dependencies:
        * erb.dat (place in the same folder as utils.pyc)
	"""
	# read data from a text file
	data=np.loadtxt(gPath + "/erb.dat",delimiter=",")

	# load centre frequencies from the first column into fc
	fc=np.array(data[:,0])
	# load bandwidths from the second column into bw
	bw=np.array(data[:,1])

	return fc,bw

def load_outMidFir_coeff():
    """
    Loads and returns 39 ERB frequencies and bandwidths
        
    Parameters:
        NONE
        
    Returns:
        * b (type: numpy array of floats) - a vector of length 4097 containing the impulse
            response of the outer middle ear.
        
    Dependencies:
        * outMidFir.dat (place in the same folder as utils.pyc)
    """

    b=np.array(np.loadtxt(gPath + "/outMidFir.dat", delimiter=","))

    return b

def load_sl_parameters():
    """
    Loads the loudness parameters for each ERB fc into a tuple. A shortcut for load_tq_data(), load_A_data() and load_alpha_data().
        
    Parameters:
        NONE
        
    Returns:
        * tq (type: numpy array of floats) - a vector of length 39 containing the threshold excitation
        intensity in quietat each ERB centre frequency
        * A (type: numpy array of floats) - a vector of length 39 containing the parameter A for each ERB fc
        * alpha (type: numpy array of floats) - a vector of length 39 containing the parameter alpha for each ERB fc
        
    Dependencies:
        * tq.dat (place in the same folder as utils.pyc)
        * A.dat (place in the same folder as utils.pyc)
        * alpha.dat (place in the same folder as utils.pyc)
    """
    tq_dB = load_tqdB_data()
    A = load_A_data()
    alpha = load_alpha_data()

    return tq_dB, A, alpha

def load_tqdB_data():
    """
    Loads and returns the excitation threshold of quiet for each ERB fc.
    
    Parameters:
        NONE
    
    Returns:
        * tq (type: numpy array of floats) - a vector of length 39 containing the threshold excitation
        intensity in quietat each ERB centre frequency
    
    Dependencies:
        * tq.dat (place in the same folder as utils.pyc)
    """

    tq = np.array(np.loadtxt(gPath + "/tq_dB.dat",delimiter=","))

    return tq

def load_A_data():
    """
    Loads and returns the excitation A parameters for each ERB fc.
        
    Parameters:
        NONE
        
    Returns:
        * A (type: numpy array of floats) - a vector of length 39 containing the parameter A for each ERB fc
        
    Dependencies:
        * A.dat (place in the same folder as utils.pyc)
    """
    
    A = np.array(np.loadtxt(gPath + "/A.dat",delimiter=","))
    
    return A

def load_alpha_data():
    """
    Loads and returns the excitation alpha parameters for each ERB fc.
        
    Parameters:
        NONE
        
    Returns:
        * alpha (type: numpy array of floats) - a vector of length 39 containing the parameter alpha for each ERB fc
        
    Dependencies:
        * alpha.dat (place in the same folder as utils.pyc)
    """
    
    alpha = np.array(np.loadtxt(gPath + "/alpha.dat",delimiter=","))
    
    return alpha

def exp_sequence(start, stop, n, base=2):
    """
    Creates a linear sequence with n points starting from start and ending at stop. For each
    element in the sequence, i, the output sequence is 2**i, generating an exponential sequence
    with base 2.
    
    Parameters:
        * start (type: numerical int) - determines the first element of the sequence, i.e.,
            base**start. (Required)
        * stop (type: numerical int) - determines the last element of the sequence, i.e.,
            base**stop. (Required)
        * n (type = numerical int) - determines the number of elements in the sequence. (Required)
        * base (type: numerical) - determines the exponential base. (Optional; Default = 2)
    
    Returns:
        * seq - the exponential sequence
    """
    
    seq = [base**x for x in np.linspace(start,stop,n)]
    
    return seq

def wavread(file):
    """
    Reads the audio data from a wav file and converts it to floating point ranging from -1 to 1.
    
    Parameters:
        * file (type: string) - the name of the wav file to read from. (Required)
        
    Returns:
        * fs (type: numerical) - the sampling frequency of the signal storied in file
        * data (type: numpy array of floats) - the data of the signal stored in file normalised
            to an amplitude range of -1 to 1.
    """

    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 
    writes the audio data to a wav file. (E.g., if 16-bit precision, then highest amplitude
    is equal to 2**16).
        
    Parameters:
        * file (type: string) - the name of the wav file to write to. (Required)
        * fs (type: numerical) - the sampling frequency of the signal. (Required)
        * data (type: array-like matrix of floats) - the signal data normalised from -1 to 1. The signal will be clipped if not
        in this range. (Required)
        * precision (type: numerical int)- the bit precision to store at. Can only be 16 or 32 bit, because that is
        all scipy allows.
        
    Returns:
        NONE
        
    TODO explore WAVE package to be allow for 24bit precision.
    """
    
    data[data>1] = 1
    data[data<-1] = -1
    
    if(precision == 16):
        dtype = np.int16
    elif(precision == 32):
        dtype = np.int32
    else:
        print "Error: precision can only be 16 or 32 bit due to scipy package."
        return

    data = nfloat_to_int(data)
    wav.write(file, fs, data)
    
    return

def plot_fft(x, xscale = 'log', yscale = 'log', show = True):
    """
    Plots the fft of signal x. If the figure is not shown, the current plot is held to allow other plots to
    be added to the same figure.
    
    Parameters:
        * x (type: array-like matrix of floats) - the signal to be analysed. (Required)
        * xscale (type: string) - the scale of the frequency axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
        * yscale (type: string) - the scale of the amplitude axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
        * show (type: boolean) - specifies whether the figure should be shown. If False, the current plot will be held so
        other plots can be added to the figure. (Optional; Default = True)
        
    Returns:
        NONE
    """

    plt.gca().set_xscale(xscale)
    plt.gca().set_yscale(yscale)
    
    x = np.array(x)
    fftx = np.absolute(fft.fft(x))
    plt.plot(range(np.shape(x)[0]),fftx)
    if(show):
        plt.show()
    else:
        plt.hold(True)

    return

def plot_waveform(x, show = True):
    """
    Plots the waveform of signal x. If the figure is not shown, the current plot is held to allow other plots to
    be added to the same figure.
        
    Parameters:
        * x (type: array-like matrix of floats) - the signal to be plotted. (Required)
        * show (type: boolean) - specifies whether the figure should be shown. If False, the current plot will be held so
        other plots can be added to the figure. (Optional; Default = True)
        
    Returns:
        NONE
    """
    
    x = np.array(x)
    plt.plot(x)
    
    if(show): plt.show()
    else: plt.hold(True)

    return

def int_to_nfloat(input, outputtype=np.float32):
    """
    Convert integer with to floating point with a range from -1 to 1.
        
    Parameters:
        * input (type: array-like matrix of ints) - a signed integer array. (Required)
        dtype : the output datatype. (Optional; Default = np.float32)
        
    Returns:
        * y (type: array-like matrix of floats) - a float array normalised to a
            range of -1 to 1.
    """
    
    input = np.array(input)
    assert input.dtype.kind == 'i', "'input' must be an array-like matrix of integers."
    outputtype = np.dtype(outputtype)
    inputdtype = np.dtype(type(input[0]))
    input = input.astype(outputtype)
    
    input[input > 0] = input[input > 0] / np.iinfo(inputdtype).max
    input[input < 0] = input[input < 0] / -np.iinfo(inputdtype).min
    
    y = input
    return y

def nfloat_to_int(input, type=np.int16):
    """
    Convert a float array with amplitude ranging from -1 to 1 to a unnormalised signed 
    integer array with specified precision.
        
    Parameters:
        * input (type: array-like matrix of floats) - a float array. (Required)
        dtype : the output datatype (also determines the precision). 
        (Optional; Default = np.int16)
        
    Returns:
        * y (type: array-like matrix of ints) - an unnormalised to a
        range of -1 to 1.
    """
    
    input = np.array(input)
    assert input.dtype.kind == 'f', "'input' must be an array of floats!"
    
    input[input > 0] = input[ input > 0 ] * np.iinfo(type).max
    input[input < 0] = input[ input < 0 ] * -np.iinfo(type).min
    
    y = input.astype(type)
    
    return y