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 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 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 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 plotSpectogramF0Segments(x, fs, w, N, H, f0, segments):
    """
    Code for plotting the f0 contour on top of the spectrogram
    """
    # frequency range to plot
    maxplotfreq = 1000.0    
    fontSize = 16

    fig = plt.figure()
    ax = fig.add_subplot(111)

    mX, pX = stft.stftAnal(x, fs, w, N, H)                      #using same params as used for analysis
    mX = np.transpose(mX[:,:int(N*(maxplotfreq/fs))+1])
    
    timeStamps = np.arange(mX.shape[1])*H/float(fs)                             
    binFreqs = np.arange(mX.shape[0])*fs/float(N)
    
    plt.pcolormesh(timeStamps, binFreqs, mX)
    plt.plot(timeStamps, f0, color = 'k', linewidth=5)

    for ii in range(segments.shape[0]):
        plt.plot(timeStamps[segments[ii,0]:segments[ii,1]], f0[segments[ii,0]:segments[ii,1]], color = '#A9E2F3', linewidth=1.5)        
    
    plt.autoscale(tight=True)
    plt.ylabel('Frequency (Hz)', fontsize = fontSize)
    plt.xlabel('Time (s)', fontsize = fontSize)
    plt.legend(('f0','segments'))
    
    xLim = ax.get_xlim()
    yLim = ax.get_ylim()
    ax.set_aspect((xLim[1]-xLim[0])/(2.0*(yLim[1]-yLim[0])))    
    plt.autoscale(tight=True) 
    plt.show()

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 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 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

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)
    """

    # read the input sound
    fs, signal = UF.wavread(inputFile)
    # compute window
    w = get_window(window, M)
    # compute the spectrum
    magnitude_frames, p = stft.stftAnal(signal, fs, w, N, H)

    # compute the boundaries for the energy bands
    k_3000 = 3000 * float(N) / fs
    k_10000 = 10000 * float(N) / fs

    # set up variables to hold the energy values
    # energy_low = 0
    # energy_high = 0

    # set up array to hold the energy values
    output_frame = []

    # loop through array and collect energy
    for frame in magnitude_frames:
        energy_low = 0
        energy_high = 0
        L = len(frame)
        for i in range(1, L):
            if i < k_3000:
                energy_low += (10 ** (frame[i] / 20)) ** 2
            elif i < k_10000 and i > k_3000:
                energy_high += (10 ** (frame[i] / 20)) ** 2

        # compute decibel value of energy
        energy_low = 10 * np.log10(energy_low)
        energy_high = 10 * np.log10(energy_high)

        output_frame.append([energy_low, energy_high])

    return np.array(output_frame)

def computeAndPlotF0(inputFile = '../sms-tools/sounds/piano.wav'):
    """
    Function to estimate fundamental frequency (f0) in an audio signal using TWM.
    Input:
        inputFile (string): wav file including the path    
    """
    window='hamming'
    M=2048
    N=2048
    H=256
    f0et=5.0
    t=-80
    minf0=100
    maxf0=300

    fs, x = UF.wavread(inputFile)                               #reading inputFile
    w  = get_window(window, M)                                  #obtaining analysis window    
    f0 = f0Detection(x, fs, w, N, H, t, minf0, maxf0, f0et)  #estimating F0

    ## Code for plotting the f0 contour on top of the spectrogram
    # frequency range to plot
    maxplotfreq = 500.0    
    fontSize = 16
    plot = 1

    fig = plt.figure()
    ax = fig.add_subplot(111)

    mX, pX = stft.stftAnal(x, w, N, H)                      #using same params as used for analysis
    mX = np.transpose(mX[:,:int(N*(maxplotfreq/fs))+1])
    
    timeStamps = np.arange(mX.shape[1])*H/float(fs)                             
    binFreqs = np.arange(mX.shape[0])*fs/float(N)
    
    plt.pcolormesh(timeStamps, binFreqs, mX)
    plt.plot(timeStamps, f0, color = 'k', linewidth=1.5)
    
    plt.autoscale(tight=True)
    plt.ylabel('Frequency (Hz)', fontsize = fontSize)
    plt.xlabel('Time (s)', fontsize = fontSize)
    plt.legend(('f0',))
    
    xLim = ax.get_xlim()
    yLim = ax.get_ylim()
    ax.set_aspect((xLim[1]-xLim[0])/(2.0*(yLim[1]-yLim[0])))    

    if plot == 1: 
        plt.autoscale(tight=True) 
        plt.show()
    elif plot == 2:                   #you can save the plot too!
        fig.tight_layout()
        fig.savefig('f0_over_Spectrogram.png', dpi=150, bbox_inches='tight')

def mainlobeTracker(inputFile = '../sms-tools/sounds/sines-440-602-hRange.wav'):
    """
    Input:
           inputFile (string): wav file including the path
    Output:
           window (string): The window type used for analysis
           t (float) = peak picking threshold (negative dB)
           tStamps (numpy array) = A Kx1 numpy array of time stamps at which the frequency components were estimated
           fTrackEst = A Kx2 numpy array of estimated frequency values, one row per time frame, one column per component
           fTrackTrue = A Kx2 numpy array of true frequency values, one row per time frame, one column per component
    """       
    # Analysis parameters: Modify values of the parameters marked XX
    window = 'blackman'                             # Window type
    t = -67                                               # threshold (negative dB)
    # window = blackman && t >= -67: Mean estimation error = [ 0.01060268  1.58192485] Hz
    # window = blackman harris && t >= -61: Mean estimation error = [ 0.01060268  1.58192485] Hz
    # ohers failed

    ### Go through the code below and understand it, do not modify anything ###   
    M = 2047                                             # Window size 
    N = 4096                                             # FFT Size
    H = 128                                              # Hop size in samples
    maxnSines = 2
    minSineDur = 0.02
    freqDevOffset = 10
    freqDevSlope = 0.001
    # read input sound
    fs, x = UF.wavread(inputFile)               
    w = get_window(window, M)                   # Compute analysis window
    tStamps = genTimeStamps(x.size, M, fs, H)   # Generate the tStamps to return
    # analyze the sound with the sinusoidal model
    fTrackEst, mTrackEst, pTrackEst = SM.sineModelAnal(x, fs, w, N, H, t, maxnSines, minSineDur, freqDevOffset, freqDevSlope)
    fTrackTrue = genTrueFreqTracks(tStamps)     # Generate the true frequency tracks
    tailF = 20                                 
    # Compute mean estimation error. 20 frames at the beginning and end not used to compute error
    meanErr = np.mean(np.abs(fTrackTrue[tailF:-tailF,:] - fTrackEst[tailF:-tailF,:]),axis=0)     
    print("Mean estimation error = " + str(meanErr) + ' Hz')      # Print the error to terminal
    # Plot the estimated and true frequency tracks
    mX, pX = stft.stftAnal(x, w, N, H)
    maxplotfreq = 900.0
    binFreq = fs * np.arange(N * maxplotfreq / fs) / N
    plt.pcolormesh(tStamps, binFreq, np.transpose(mX[:,:np.int(N * maxplotfreq / fs + 1)]), cmap='hot_r')
    plt.plot(tStamps,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.title('frequency detection: Window = ' + window + '& t = ' + str(t))
    plt.xlabel('Time (s)')
    plt.ylabel('Frequency (Hz)')
    plt.autoscale(tight=True)
    return window, float(t), tStamps, fTrackEst, fTrackTrue  # Output returned 

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, pX) = stft.stftAnal(x, fs, w, N, H)

    numFrames = int(mX[:,0].size)
    frmTime = H*np.arange(numFrames)/float(fs)
    binFreq = np.arange(N/2+1)*float(fs)/N

    cutoff1 = 3000
    cutoff2 = 10000

    cutoff_bucket1 = np.ceil(float(cutoff1) * N / fs)
    cutoff_bucket2 = np.ceil(float(cutoff2) * N / fs)

    low_band = mX[:,1:cutoff_bucket1]
    high_band = mX[:,cutoff_bucket1:cutoff_bucket2]

    E = np.zeros((numFrames, 2))
    E[:,0] = by_frame_energy(low_band)
    E[:,1] = by_frame_energy(high_band)

    O = np.zeros((numFrames, 2))
    O[1:,:] = E[1:,:] - E[:-1,:]

    # half wave rectification
    O[O<=0] = 0

    # plot_odf(mX, fs, inputFile, M, N, H, O)

    return O

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
    lowf_init = 0
    lowf_end = 3000
    highf_init = 3000
    highf_end = 10000

    w = get_window(window, M)
    fs, x = UF.wavread(inputFile)

    lowb_init = lowf_init * N / fs
    lowb_end = lowf_end * N / fs
    highb_init = highf_init * N / fs
    highb_end = highf_end * N / fs

    xmX, pmX = stft.stftAnal(x, fs, w, N, H)

    xmX_linear = 10**(xmX / 20)

    result_low  = 10 * np.log10 ( np.sum( abs( xmX_linear[:, 1 : lowb_end] )**2, 1 ) )
    result_high = 10 * np.log10 ( np.sum( abs( xmX_linear[:, highb_init + 1 : highb_end] )**2, 1 ) )


    frames = result_low.shape[0]
    result = np.array([result_low[0], result_high[0]])


    for i in range(1, frames):
        
        temp = np.array([result_low[i], result_high[i]])
        result = np.vstack( (result, temp) )

    return result
    

def chirpTracker(inputFile='../sms-tools/sounds/chirp-150-190-linear.wav'):
    """
    Input:
           inputFile (string) = wav file including the path
    Output:
           M (int) = Window length
           H (int) = hop size in samples
           tStamps (numpy array) = A Kx1 numpy array of time stamps at which the frequency components were estimated
           fTrackEst (numpy array) = A Kx2 numpy array of estimated frequency values, one row per time frame, one column per component
           fTrackTrue (numpy array) = A Kx2 numpy array of true frequency values, one row per time frame, one column per component
           K is the number of frames
    """
    # Analysis parameters: Modify values of the parameters marked XX
    M = 3298                                    # Window size in samples
    
    ### Go through the code below and understand it, do not modify anything ###    
    H = 128                                     # Hop size in samples
    N = int(pow(2, np.ceil(np.log2(M))))        # FFT Size, power of 2 larger than M
    t = -80.0                                   # threshold
    window = 'blackman'                         # Window type
    maxnSines = 2                               # Maximum number of sinusoids at any time frame
    minSineDur = 0.0                            # minimum duration set to zero to not do tracking
    freqDevOffset = 30                          # minimum frequency deviation at 0Hz
    freqDevSlope = 0.001                        # slope increase of minimum frequency deviation
    
    fs, x = UF.wavread(inputFile)               # read input sound
    w = get_window(window, M)                   # Compute analysis window
    tStamps = genTimeStamps(x.size, M, fs, H)   # Generate the tStamps to return
    # analyze the sound with the sinusoidal model
    fTrackEst, mTrackEst, pTrackEst = SM.sineModelAnal(x, fs, w, N, H, t, maxnSines, minSineDur, freqDevOffset, freqDevSlope)
    fTrackTrue = genTrueFreqTracks(tStamps)     # Generate the true frequency tracks
    tailF = 20                                 
    # Compute mean estimation error. 20 frames at the beginning and end not used to compute error
    meanErr = np.mean(np.abs(fTrackTrue[tailF:-tailF,:] - fTrackEst[tailF:-tailF,:]),axis=0)     
    print("Mean estimation error = " + str(meanErr) + ' Hz')      # Print the error to terminal    
    # Plot the estimated and true frequency tracks
    mX, pX = stft.stftAnal(x, w, N, H)  # stft from anal
    maxplotfreq = 1500.0
    binFreq = fs*np.arange(N*maxplotfreq/fs)/N
    plt.pcolormesh(tStamps, binFreq, np.transpose(mX[:,:int(N * maxplotfreq / fs + 1)]),cmap = 'hot_r')
    plt.plot(tStamps,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.title('True and estimated frequency, windowsize = ' + str(M))
    plt.xlabel('Time (s)')
    plt.ylabel('Frequency (Hz)')
    plt.autoscale(tight=True)
    plt.show()
    return M, H, tStamps, fTrackEst, fTrackTrue  # Output returned 

def mainlobeTracker(inputFile="../../sounds/sines-440-602-hRange.wav"):
    """
    Input:
           inputFile (string): wav file including the path
    Output:
           window (string): The window type used for analysis
           t (float) = peak picking threshold (negative dB)
           tStamps (numpy array) = A Kx1 numpy array of time stamps at which the frequency components were estimated
           fTrackEst = A Kx2 numpy array of estimated frequency values, one row per time frame, one column per component
           fTrackTrue = A Kx2 numpy array of true frequency values, one row per time frame, one column per component
    """
    # Analysis parameters: Modify values of the parameters marked XX
    window = "blackmanharris"  # Window type
    t = -93.0  # threshold (negative dB)

    ### Go through the code below and understand it, do not modify anything ###
    M = 2047  # Window size
    N = 4096  # FFT Size
    H = 128  # Hop size in samples
    maxnSines = 2
    minSineDur = 0.02
    freqDevOffset = 10
    freqDevSlope = 0.001
    # read input sound
    fs, x = UF.wavread(inputFile)
    w = get_window(window, M)  # Compute analysis window
    tStamps = genTimeStamps(x.size, M, fs, H)  # Generate the tStamps to return
    # analyze the sound with the sinusoidal model
    fTrackEst, mTrackEst, pTrackEst = SM.sineModelAnal(
        x, fs, w, N, H, t, maxnSines, minSineDur, freqDevOffset, freqDevSlope
    )
    fTrackTrue = genTrueFreqTracks(tStamps)  # Generate the true frequency tracks
    tailF = 10
    # Compute mean estimation error. 50 frames at the beginning and end not used to compute error
    meanErr = np.mean(np.abs(fTrackTrue[tailF:-tailF, :] - fTrackEst[tailF:-tailF, :]), axis=0)
    print "Mean estimation error = " + str(meanErr) + " Hz"  # Print the error to terminal
    # Plot the estimated and true frequency tracks
    mX, pX = stft.stftAnal(x, fs, w, N, H)
    maxplotfreq = 900.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(tStamps, 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 window, t, tStamps, fTrackEst, fTrackTrue  # Output returned

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,pX = stft.stftAnal(x,w,N,H)
    mX = pow(10,mX/20.)
    
    num_frames = len(mX)
    band_energy = np.zeros((len(mX),2))
    for frm_idx in range(num_frames):
        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)

    odf = np.zeros((num_frames,2))
    for frm_idx in range(1,num_frames):
        odf[frm_idx,0] = band_energy[frm_idx,0]-band_energy[frm_idx-1,0]
        odf[frm_idx,0] = 0 if  odf[frm_idx,0] < 0 else odf[frm_idx,0]
        odf[frm_idx,1] = band_energy[frm_idx,1]-band_energy[frm_idx-1,1]
        odf[frm_idx,1] = 0 if  odf[frm_idx,1] < 0 else odf[frm_idx,1]

    return odf

def main(inputFile='../../sounds/bendir.wav', window='hamming', M=2001, N=2048, t=-80, 
	minSineDur=0.02, maxnSines=150, freqDevOffset=10, freqDevSlope=0.001):
	"""
	inputFile: input sound file (monophonic with sampling rate of 44100)
	window: analysis window type (rectangular, hanning, hamming, blackman, blackmanharris)	
	M: analysis window size 
	N: fft size (power of two, bigger or equal than M)
	t: magnitude threshold of spectral peaks 
	minSineDur: minimum duration of sinusoidal tracks
	maxnSines: maximum number of parallel sinusoids
	freqDevOffset: frequency deviation allowed in the sinusoids from frame to frame at frequency 0   
	freqDevSlope: slope of the frequency deviation, higher frequencies have bigger deviation
	"""

	# size of fft used in synthesis
	Ns = 512

	# hop size (has to be 1/4 of Ns)
	H = 128

	# read input sound
	(fs, x) = UF.wavread(inputFile)

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

	# perform sinusoidal plus residual analysis
	tfreq, tmag, tphase, xr = SPR.sprModelAnal(x, fs, w, N, H, t, minSineDur, maxnSines, freqDevOffset, freqDevSlope)
		
	# compute spectrogram of residual
	mXr, pXr = STFT.stftAnal(xr, fs, w, N, H)

	# sum sinusoids and residual
	y, ys = SPR.sprModelSynth(tfreq, tmag, tphase, xr, Ns, H, fs)

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

	# write sounds files for sinusoidal, residual, and the sum
	UF.wavwrite(ys, fs, outputFileSines)
	UF.wavwrite(xr, fs, outputFileResidual)
	UF.wavwrite(y, fs, outputFile)
	return x, fs, mXr, tfreq, y

def main(inputFile='../../sounds/sax-phrase-short.wav', window='blackman', M=601, N=1024, t=-100, 
	minSineDur=0.1, nH=100, minf0=350, maxf0=700, f0et=5, harmDevSlope=0.01):
	"""
	Perform analysis/synthesis using the harmonic plus residual model
	inputFile: input sound file (monophonic with sampling rate of 44100)
	window: analysis window type (rectangular, hanning, hamming, blackman, blackmanharris)	
	M: analysis window size; N: fft size (power of two, bigger or equal than M)
	t: magnitude threshold of spectral peaks; minSineDur: minimum duration of sinusoidal tracks
	nH: maximum number of harmonics; minf0: minimum fundamental frequency in sound
	maxf0: maximum fundamental frequency in sound; f0et: maximum error accepted in f0 detection algorithm                                                                                            
	harmDevSlope: allowed deviation of harmonic tracks, higher harmonics have higher allowed deviation
	"""

	# size of fft used in synthesis
	Ns = 512

	# hop size (has to be 1/4 of Ns)
	H = 128

	# read input sound
	(fs, x) = UF.wavread(inputFile)

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

	# find harmonics and residual
	hfreq, hmag, hphase, xr = HPR.hprModelAnal(x, fs, w, N, H, t, minSineDur, nH, minf0, maxf0, f0et, harmDevSlope)
	  
	# compute spectrogram of residual
	mXr, pXr = STFT.stftAnal(xr, fs, w, N, H)
	  
	# synthesize hpr model
	y, yh = HPR.hprModelSynth(hfreq, hmag, hphase, xr, Ns, H, fs)

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

	# write sounds files for harmonics, residual, and the sum
	UF.wavwrite(yh, fs, outputFileSines)
	UF.wavwrite(xr, fs, outputFileResidual)
	UF.wavwrite(y, fs, outputFile)
	return x, fs, mXr,hfreq, 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
    w = get_window(window, M)         # get the window
    (fs,x) = UF.wavread(inputFile) 
    lowfreq = 3000.0 
    highfreq = 10000.0
    array_size = int( math.ceil (float(x.size) / H ) )#H = 128
    k1 =  math.ceil ( lowfreq / ( float(fs) / N ) ) 
    k2 =  math.ceil ( highfreq / ( float(fs) / N ) )
    energysignal = 0 
    energysign