def writeExampleFiles():
    """
    A convenience function: writes out example files, some of them with optimal parameters found by exploreSineModelMultiRes()
    """
    inputFile='../../sounds/orchestra.wav'
    fs, x = UF.wavread(inputFile)
    W = np.array(['blackmanharris'])
    M = np.array([1001])
    N = np.array([4096])
    B = np.array([ ])
    T = np.array([-90])
    Ns = 512
    best = Best()
    y = best.calculateAndUpdate(x, fs, Ns, W, M, N, B, T)
    outputFile = inputFile[:-4] + '_optimizedSineModel.wav'
    print '->',outputFile
    UF.wavwrite(y, fs, outputFile)

    inputFile='../../sounds/121061__thirsk__160-link-strings-2-mono.wav'
    fs, x = UF.wavread(inputFile)
    W = np.array(['hamming','hamming','hamming'])
    M = np.array([3001,1501,751])
    N = np.array([16384,8192,4096])
    B = np.array([2756.25,5512.5])
    T = np.array([-90,-90,-90])
    Ns = 512
    best = Best()
    y = best.calculateAndUpdate(x, fs, Ns, W, M, N, B, T)
    outputFile = inputFile[:-4] + '_optimizedSineModel.wav'
    print '->',outputFile
    UF.wavwrite(y, fs, outputFile)

    inputFile='../../sounds/orchestra.wav'
    fs, x = UF.wavread(inputFile)
    W = np.array(['hamming','hamming','hamming'])
    M = np.array([3001,1501,751])
    N = np.array([16384,8192,4096])
    B = np.array([2756.25,5512.5])
    T = np.array([-90,-90,-90])
    Ns = 512
    best = Best()
    y = best.calculateAndUpdate(x, fs, Ns, W, M, N, B, T)
    outputFile = inputFile[:-4] + '_nonOptimizedSineModel.wav'
    print '->',outputFile
    UF.wavwrite(y, fs, outputFile)

    inputFile='../../sounds/121061__thirsk__160-link-strings-2-mono.wav'
    fs, x = UF.wavread(inputFile)
    W = np.array(['blackmanharris'])
    M = np.array([1001])
    N = np.array([4096])
    B = np.array([ ])
    T = np.array([-90])
    Ns = 512
    best = Best()
    y = best.calculateAndUpdate(x, fs, Ns, W, M, N, B, T)
    outputFile = inputFile[:-4] + '_nonOptimizedSineModel.wav'
    print '->',outputFile
    UF.wavwrite(y, fs, outputFile)

def computeODF(inputFile, window, M, N, H):
    """
    Inputs:
            inputFile (string): input sound file (monophonic with sampling rate of 44100)
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                blackman, blackmanharris)
            M (integer): analysis window size (odd integer value)
            N (integer): fft size (power of two, bigger or equal than than M)
            H (integer): hop size for the STFT computation
    Output:
            The function should return a numpy array with two columns, where the first column is the ODF 
            computed on the low frequency band and the second column is the ODF computed on the high 
            frequency band.
            ODF[:,0]: ODF computed in band 0 < f < 3000 Hz 
            ODF[:,1]: ODF computed in band 3000 < f < 10000 Hz
    """
    
    ### your code here
    windowing = get_window(window, M)

    (fs, x) = UF.wavread(inputFile)

    mX, pX = stft.stftAnal(x, fs, windowing, N, H)

    bin0 = 1
    bin3000 = np.floor(3000.0*N/fs)
    bin10000 = np.floor(10000.0*N/fs)
    bin3000up = np.ceil(3000.0*N/fs)

    ODF = np.zeros((mX.shape[0], 2))

    prevODF3000 = 0.0
    prevODF10000 = 0.0

    for i in range(mX.shape[0]):
        env3000 = np.sum(np.square(10**(mX[i,1:bin3000+1] / 20)))
        env3000db = 10 * np.log10(env3000)

        odf3000 = env3000db - prevODF3000
        prevODF3000 = env3000db
 
        if odf3000 <= 0.0:
            odf3000 = 0.0

        ODF[i,0] = odf3000
        

        env10000 = np.sum(np.square(10**(mX[i,bin3000up:bin10000+1] / 20)))
        env10000db = 10 * np.log10(env10000)

        odf10000 = env10000db - prevODF10000
        prevODF10000 = env10000db

        if odf10000 <= 0.0:
            odf10000 = 0.0

        ODF[i,1] = odf10000

        
    return ODF

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    ## your code here

    def energy(X, k1, k2):
        X2 = np.power(X, 2)
        return np.sum(X2[k1:k2])

    fs, x = UF.wavread(inputFile)
    w = get_window(window, M)
    xsyn = stft.stft(x, fs, w, N, H)
    noise = np.subtract(xsyn, x)
    
    Esignal1 = energy(x, 0, len(x))
    Enoise1 = energy(noise, 0, len(noise))
    SNR1 = 10*np.log10(Esignal1/Enoise1)
    
    Esignal2 = energy(x, M+1, len(x)-M-1)
    Enoise2 = energy(noise, M+1, len(noise)-M-1)
    SNR2 = 10*np.log10(Esignal2/Enoise2)

    return SNR1, SNR2

def sineModelOriginalTest1(inputFile, M):
    
    print "\n\n\n###############  RUN THE ORIGINAL TEST (without multiresolution)  ###############\n"
    
    #M1 = 4095
    M1 = M
    
    print "M: "
    print M
    
    N1 = int(pow(2, np.ceil(np.log2(M1))))      # FFT Size, power of 2 larger than M

    print "N1: "
    print N1
    
    t = -80.0                                   # threshold

    fs, x = UF.wavread(inputFile)               # read input sound
    
    #print "Ploting \"x\""
    #plt.plot(x)
    
    window = 'blackman'                         # Window type
    w1 = get_window(window, M1)                 # compute analysis window
            
    return sineModelOriginal(x,fs,w1,N1,t)

def phaseFlux(inFile, window, M, bins, H, passes, th, inhibTh, inhibRel, plot=False):
    fs, x = UF.wavread(inFile)  # read file
    x = normalise(x)            # normalise
    if plot:
        plt.plot(-transient(x, 10))

    N = len(x)                  # length of file
    win = get_window(window, M) # create window
   

    # STFT:
    X = np.ndarray(shape=(N/H, bins), dtype='complex')
    for n in range(N/H):
        Xpart = x[n * H:n * H + M]
        if len(Xpart) < len(win):
            Xpart = zeropad(Xpart, len(win))
        X[n] = UF.fft(Xpart * win, bins)    # gets auto zerophased and padded

    '''
    bins:               0   1   2   3   .. fftSize
timefrms:      
            0       [   .   .   .   .   ..  .   
            1           .   .   .   .   ..  .   
            2           .   .   .   .   ..  .
            ..          ..  ..  ..  ..  ..  ..  ]
            N/H
    '''
    
    mX = np.abs(X)  # get magnitude
    mX = normalise(mX)
    pX = np.angle(X) # get phase
    pX = normalise(pX)

    pX_uw = nd_unwrapPhase(pX, 0); # unwrapPhase
    # pX_uw = pX 

    derv1 = nd_derivative(pX_uw, 0)
    derv2 = nd_derivative(derv1, 0)
    
    binmul = np.ndarray(shape=(N/H, bins), dtype='float')
    for n in range(N/H):
        for k in range(bins):
            binmul[n][k] = 1

    onsets = np.ndarray(shape=(passes, N), dtype='float')
    for p in range(passes):
        for n in range(N/H):
            val = np.sum(derv2[n] * mX[n])
            if val / bins > th:
                onsets[p][n*H] = val
                for k in range(bins):
                    if derv2[n][k] > inhibTh:
                        for m in range(inhibRel):
                            binmul[n+m][k] = m/inhibRel 
                            if n+m > N/H:
                                break;
    
            mX *= binmul;
    # onsets = np.transpose(onsets)
    return normalise(onsets)

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    ## your code here
    def energy(mag):
        e = np.sum((10 ** (mag / 20)) ** 2)
        return e
    
    (fs, x) = UF.wavread(inputFile)
    w = get_window(window, M)
    
    mX, pX = STFT.stftAnal(x, fs, w, N, H)
    y = STFT.stftSynth(mX, pX, M, H)
    n = x - y[:x.size]
    n2 = x[w.size:-w.size] - y[:x.size][w.size:-w.size]
    
    mN, pN = STFT.stftAnal(n, fs, w, N, H)
    mN2, pN2 = STFT.stftAnal(n2, fs, w, N, H)
    
    snr1 = 10 * np.log10(energy(mX) / energy(mN))
    snr2 = 10 * np.log10(energy(mX) / energy(mN2))
        
    return snr1, snr2

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    x = UF.wavread(inputFile)[1]
    w = get_window(window, M)

    xSynth = stft(x, 1.0, w, N, H)

    eSignal1 = sum(x**2)
    eNoise1 = sum((x-xSynth)**2)
    SNR1 = 10.0*np.log10(eSignal1/eNoise1)

    x2 = x[M:len(x)-M]
    xSynth2 = xSynth[M:len(xSynth)-M]

    eSignal2 = sum(x2**2)
    eNoise2 = sum((x2-xSynth2)**2)
    SNR2 = 10.0*np.log10(eSignal2/eNoise2)

    return (SNR1,SNR2)

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """

    #read from the file
    FS, x = UF.wavread(inputFile)

    w = get_window(window, M)
    #do a stft computation
    y = stft.stft(x, FS, w, N, H)

    #compute SNR over complete signal
    diff = y - x
    energy_signal = (y**2).sum()
    energy_noise = (diff**2).sum()
    SNR1 = 10 * np.log10(energy_signal/energy_noise)

    #compute SNR over sliced signal
    energy_signal_sliced = (y[M:-M]**2).sum()
    energy_noise_sliced = (diff[M:-M]**2).sum()
    SNR2 = 10 * np.log10(energy_signal_sliced/energy_noise_sliced)


    return (SNR1, SNR2)

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    ## your code here
    fs, x = UF.wavread(inputFile)

    w = get_window(window, M)

    y = stft.stft(x, fs, w, N, H)
    noise = np.array(x - y)

    E_x = np.sum( abs(x)**2 )
    E_noise = np.sum( abs(noise)**2 )

    E_xAfterM = np.sum( abs( x[M : x.size-M] )**2 )
    E_nAfterM = np.sum( abs( noise[M : x.size-M] )**2 )

    SNR1 = 10 * np.log10(E_x / E_noise)
    SNR2 = 10 * np.log10(E_xAfterM/E_nAfterM)

    return (SNR1, SNR2)

def computeODF(inputFile, window, M, N, H):
    """
    Inputs:
            inputFile (string): input sound file (monophonic with sampling rate of 44100)
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                blackman, blackmanharris)
            M (integer): analysis window size (odd integer value)
            N (integer): fft size (power of two, bigger or equal than than M)
            H (integer): hop size for the STFT computation
    Output:
            The function should return a numpy array with two columns, where the first column is the ODF 
            computed on the low frequency band and the second column is the ODF computed on the high 
            frequency band.
            ODF[:,0]: ODF computed in band 0 < f < 3000 Hz 
            ODF[:,1]: ODF computed in band 3000 < f < 10000 Hz
    """
    ### your code here
    fs, x = UF.wavread(inputFile)

    w = get_window(window, M)

    mX = stft.stftAnal(x, fs, w, N, H)[0]

    X = 10 ** (mX / 20.0)

    b3k = int(N*3000.0/fs)
    b10k = int(N*10000.0/fs)

    o3k = odf(X[:, 1:b3k+1])
    o10k = odf(X[:, b3k+1:b10k+1])

    return np.column_stack((o3k, o10k))

def computeModel(inputFile, B, M, window = 'hanning', t = -90):

    bands = range(len(B))

    fs, x = UF.wavread(inputFile)
    w = [get_window(window, M[i]) for i in bands]
    N = (2**np.ceil(np.log2(B))).astype(int)

    y_combined = SMMR.sineModelMultiRes(x, fs, w, N, t, B)

    #y, y_combined = SMMR.sineModelMultiRes_combined(x, fs, w, N, t, B)

    # output sound file name
    outputFileInputFile = 'output_sounds/' + os.path.basename(inputFile)
    #outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_sineModel.wav'
    outputFile_combined = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_sineModelMultiRes.wav'

    # write the synthesized sound obtained from the sinusoidal synthesis
    UF.wavwrite(x, fs, outputFileInputFile)
    #UF.wavwrite(y, fs, outputFile)
    UF.wavwrite(y_combined, fs, outputFile_combined)

    plt.figure()
    plt.plot(x)
    plt.plot(y_combined)
    plt.show()

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    ## your code here
    w = get_window(window, M)  # get the window
    (fs, x) = UF.wavread(inputFile)
    # x: input sound, w: analysis window, N: FFT size, H: hop size
    # returns y: output sound
    STFTX = stft.stft(x, fs, w, N, H)
    xoutput = np.arange(x.size)
    energynoise = 0
    energynoise2 = 0
    for i in range(0, x.size):
        energynoise += np.power(np.abs(x[i].real) - np.abs(STFTX[i].real), 2)
        if i > M and i < x.size - M:
            energynoise2 += np.power(np.abs(x[i].real) - np.abs(STFTX[i].real), 2)
    energysignal = 0
    energysignal2 = 0
    for i in range(0, x.size):
        energysignal += np.power(np.abs(x[i].real), 2)
        if i > M and i < x.size - M:
            energysignal2 += np.power(np.abs(x[i].real), 2)
    SNR1 = 10 * np.log10(energysignal / energynoise)
    SNR2 = 10 * np.log10(energysignal2 / energynoise2)
    return SNR1, SNR2

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    
    ## your code here
    
    fs, x = UF.wavread(inputFile)
    
    w = get_window(window, M)
    
    xrec = stft.stft(x, fs, w, N, H)
    
    eSignal = energy(x)
    eSignal_part = energy(x[M:-M])
    eNoise = energy(x-xrec)
    eNoise_part = energy((x-xrec)[M:-M])

    snr = 10 * np.log10(eSignal / eNoise)
    snr_part = 10 * np.log10(eSignal_part / eNoise_part)
    
    return snr, snr_part

def minFreqEstErr(inputFile, f):
    """
    Inputs:
            inputFile (string) = wav file including the path
            f (float) = frequency of the sinusoid present in the input audio signal (Hz)
    Output:
            fEst (float) = Estimated frequency of the sinusoid (Hz)
            M (int) = Window size
            N (int) = FFT size
    """
    # analysis parameters:
    window = 'blackman'
    t = -40
    fs, x = UF.wavread(inputFile)
    x_half = len(x) // 2
    f_error = np.inf
    k = 1
    while f_error > 0.05:  # Hz
        M = 100 * k + 1
        M2 = M // 2
        W = get_window(window, M)
        N = int(2 ** np.ceil(np.log2(M)))
        mX, pX = DFT.dftAnal(x[x_half - M2: x_half - M2 + M], W, N)
        ploc = UF.peakDetection(mX, t)
        iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)
        fEst = iploc * fs / N
        f_error = np.abs(f - fEst)
        k += 1
    return(fEst, M, N)

def sineODF(file='../../../../../audioDSP_course/assignments/sms-tools/sounds/piano.wav'):
    fs, x = UF.wavread(file)

    # set params:
    M = 1024    # window size
    H = int(M/3)     # hop size
    t = -80.0   #treshold (dB??)
    window = 'blackman' # window type
    fftSize = int(pow(2, np.ceil(np.log2(M))))  # size of FFT
    N = fftSize
    maxnSines = 10      # maximum simultaneous sines
    minSineDur = 0.1    # minimal duration of sines
    freqDevOffset = 30  # min(??) frequency deviation at 0Hz
    freqDevSlope = 0.001    # slope increase of min freq dev.


    w = get_window(window, M)    # get analysis window
    tStamps = genTimeStamps(len(x), M, fs, H)    # generate timestamp return?
    fTrackEst, mTrackEst, pTreckEst = SM.sineModelAnal(x, fs, w, fftSize, H, t, maxnSines, minSineDur, freqDevOffset, freqDevSlope)

    fTrackTrue = genTrueFreqTracks(tStamps) # get true freq. tracks

    # plotting:
    mX, pX = stft.stftAnal(x, fs, w, fftSize, H)
    maxplotfreq = 1500.0
    binFreq = fs*np.arange(N*maxplotfreq/fs)/N
    plt.pcolormesh(tStamps, binFreq, np.transpose(mX[:,:N*maxplotfreq/fs+1]),cmap = 'hot_r')
    # plt.plot(fTrackTrue, 'o-', color = 'c', linewidth=3.0)
    plt.plot(tStamps, fTrackEst, color = 'y', linewidth=2.0)
    # plt.legend(('True f1', 'True f2', 'Estimated f1', 'Estimated f2'))
    plt.xlabel('Time (s)')
    plt.ylabel('Frequency (Hz)')
    plt.autoscale(tight=True)
    return fTrackEst

def computeEngEnv(inputFile, window, M, N, H):
    """
    Inputs:
            inputFile (string): input sound file (monophonic with sampling rate of 44100)
            window (string): analysis window type (choice of rectangular, triangular, hanning, 
                hamming, blackman, blackmanharris)
            M (integer): analysis window size (odd positive integer)
            N (integer): FFT size (power of 2, such that N > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a numpy array engEnv with shape Kx2, K = Number of frames
            containing energy envelop of the signal in decibles (dB) scale
            engEnv[:,0]: Energy envelope in band 0 < f < 3000 Hz (in dB)
            engEnv[:,1]: Energy envelope in band 3000 < f < 10000 Hz (in dB)
    """
    
    (fs,x) = UF.wavread(inputFile)
    w = get_window(window, M)
    (xmX, xpX) = stft.stftAnal(x, fs, w, N, H)

    kLow1 = 0
    
    kLow2 = 0
    while (True):
	kLow2 += 1
	if( (kLow2 < N*(fLow2)/float(fs)) & (kLow2 > N*(fLow2)/float(fs) - 1.0 ) ):
	    break
    
    kHigh1 = 0
    while (True):
	kHigh1 += 1
	if( (kHigh1 < N*(fHigh1)/float(fs)) & (kHigh1 > N*(fHigh1)/float(fs) - 1.0 ) ):
	    break
    
    kHigh2 = 0
    while (True):
	kHigh2 += 1
	if( (kHigh2 < N*(fHigh2)/float(fs)) & (kHigh2 > N*(fHigh2)/float(fs) - 1.0 ) ):
	    break
    
    nHops = int(xmX.shape[0])
    out = np.zeros((nHops,2))
    
    i = 0
    while i < nHops:
        subxmX = xmX[i,:]
    
        subLowxmX = subxmX[kLow1+1:kLow2+1]
        subLowxmX = 10**(subLowxmX/20)
        eSignalLow = sum(subLowxmX**2)
        out[i,0] = 10.0*np.log10(eSignalLow)

        subHighxmX = subxmX[kHigh1+1:kHigh2+1]
        subHighxmX = 10**(subHighxmX/20)
        eSignalHigh = sum(subHighxmX**2)
        out[i,1] = 10.0*np.log10(eSignalHigh)
    
        i += 1 

    return out

def getJawaab(ipFile = '../dataset/testInputs/testInput_1.wav', ipulsePos = getPulsePosFromAnn('../dataset/testInputs/testInput_1.csv'), strokeModels = None, oFile = './tablaOutput.wav', randomFlag = 1):
    # If poolFeats are not built, give an error!
    if strokeModels == None:
        print "Train models first before calling getJawaab() ..."
        opulsePos = None
        strokeSeq = None
        oFile = None
        ts = None
    else:
        print "Getting jawaab..."
        pulsePeriod = np.median(np.diff(ipulsePos))
        print pulsePeriod
        fss, audioIn = UF.wavread(ipFile)
        if randomFlag == 1:
            strokeSeq, tStamps, opulsePos = genRandomComposition(pulsePeriod, pieceDur = len(audioIn)/params.Fs, strokeModels = strokeModels)
        else:
            invCmat = getInvCovarianceMatrix(strokeModels)
            strokeSeq, tStamps, opulsePos = genSimilarComposition(pulsePeriod, pieceDur = len(audioIn)/params.Fs, strokeModels = strokeModels, iAudioFile = ipFile, iPos = ipulsePos,invC = invCmat)
        print strokeSeq
        print tStamps
        print opulsePos
        if oFile != None:
            audio = genAudioFromStrokeSeq(strokeModels,strokeSeq,tStamps)
            audio = audio/(np.max(audio) + 0.01)
            UF.wavwrite(audio, params.Fs, oFile)
    return opulsePos, strokeSeq, tStamps, oFile

def main(inputFile = '../../sounds/piano.wav', window = 'hamming', M = 1024, N = 1024, H = 512):
	"""
	analysis/synthesis using the STFT
	inputFile: input sound file (monophonic with sampling rate of 44100)
	window: analysis window type (choice of rectangular, hanning, hamming, blackman, blackmanharris)	
	M: analysis window size 
	N: fft size (power of two, bigger or equal than M)  
	H: hop size (at least 1/2 of analysis window size to have good overlap-add)               
	"""

	# read input sound (monophonic with sampling rate of 44100)
	fs, x = UF.wavread(inputFile)

	# compute analysis window
	w = get_window(window, M)

	# compute the magnitude and phase spectrogram
	mX, pX = STFT.stftAnal(x, fs, w, N, H)
	 
	# perform the inverse stft
	y = STFT.stftSynth(mX, pX, M, H)

	# output sound file (monophonic with sampling rate of 44100)
	outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_stft.wav'   

	# write the sound resulting from the inverse stft
	UF.wavwrite(y, fs, outputFile)
	return x, fs, mX, pX, y

def computeEngEnv(inputFile, window, M, N, H):
    """
    Inputs:
            inputFile (string): input sound file (monophonic with sampling rate of 44100)
            window (string): analysis window type (choice of rectangular, triangular, hanning, 
                hamming, blackman, blackmanharris)
            M (integer): analysis window size (odd positive integer)
            N (integer): FFT size (power of 2, such that N > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a numpy array engEnv with shape Kx2, K = Number of frames
            containing energy envelop of the signal in decibles (dB) scale
            engEnv[:,0]: Energy envelope in band 0 < f < 3000 Hz (in dB)
            engEnv[:,1]: Energy envelope in band 3000 < f < 10000 Hz (in dB)
    """
    
    ### your code here
    fs,x = UF.wavread(inputFile)
    w = get_window(window,M)

    mX,pX = stft.stftAnal(x,w,N,H)
    mX = pow(10,mX/20.)
    
    band_energy = np.zeros((len(mX),2))
    for frm_idx in range(len(mX)):
        frm = mX[frm_idx]
        for k in range(len(frm)):
            cur_f = k*44100/N
            if cur_f > 0 and cur_f < 3000:
                band_energy[frm_idx,0] += (frm[k]*frm[k])
            elif cur_f > 3000 and cur_f < 10000:
                band_energy[frm_idx,1] += (frm[k]*frm[k])

    band_energy = 10.0*np.log10(band_energy)
    return band_energy

def computeEngEnv(inputFile, window, M, N, H):
    """
    Inputs:
            inputFile (string): input sound file (monophonic with sampling rate of 44100)
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                blackman, blackmanharris)
            M (integer): analysis window size (odd positive integer)
            N (integer): FFT size (power of 2, such that N > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a numpy array engEnv with shape Kx2, K = Number of frames
            containing energy envelop of the signal in decibles (dB) scale
            engEnv[:,0]: Energy envelope in band 0 < f < 3000 Hz (in dB)
            engEnv[:,1]: Energy envelope in band 3000 < f < 10000 Hz (in dB)
    """
    
    ### your code here
    def energy(mag):
        e = 10 * np.log10(np.sum((10 ** (mag / 20)) ** 2, axis=1))
        return e
    
    (fs, x) = UF.wavread(inputFile)
    border_bin = int(np.ceil(float(3000) * N / fs))
    max_bin = int(np.ceil(float(10000) * N / fs))
    w = get_window(window, M)
    
    mX, pX = STFT.stftAnal(x, fs, w, N, H)
    low = np.transpose(np.transpose(mX)[1:border_bin])
    high = np.transpose(np.transpose(mX)[border_bin:max_bin])
    
    e_low = energy(low)
    e_high = energy(high)
    
    envs = np.append([e_low], [e_high], axis=0)
    envs = np.transpose(envs)
    
    # draw graph
    plt.figure(1, figsize=(9.5, 6))

    plt.subplot(211)
    numFrames = mX.shape[0]
    frmTime = H*np.arange(numFrames)/float(fs)
    binFreq = np.arange(mX.shape[1])*float(fs)/N
    plt.pcolormesh(frmTime, binFreq, np.transpose(mX))
    plt.title('mX ({0}), M={1}, N={2}, H={3}'.format(inputFile, M, N, H))
    plt.autoscale(tight=True)
    
    plt.subplot(212)
    plt.plot(frmTime, e_low, color="blue", label="row")
    plt.plot(frmTime, e_high, color="red", label="high")
    plt.title('Energy of Envelopes')
    plt.autoscale(tight=True)

    plt.tight_layout()
    plt.show()
    
    return envs