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 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 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 transformation_synthesis(inputFile, fs, hfreq, hmag, freqScaling = np.array([0, 2.0, 1, .3]), 
	freqStretching = np.array([0, 1, 1, 1.5]), timbrePreservation = 1, 
	timeScaling = np.array([0, .0, .671, .671, 1.978, 1.978+1.0])):
	# transform the analysis values returned by the analysis function and synthesize the sound
	# inputFile: name of input file
	# fs: sampling rate of input file	
	# tfreq, tmag: sinusoidal frequencies and magnitudes
	# freqScaling: frequency scaling factors, in time-value pairs
	# freqStretchig: frequency stretching factors, in time-value pairs
	# timbrePreservation: 1 preserves original timbre, 0 it does not
	# timeScaling: time scaling factors, in time-value pairs

	# size of fft used in synthesis
	Ns = 512

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

	# frequency scaling of the harmonics 
	yhfreq, yhmag = HT.harmonicFreqScaling(hfreq, hmag, freqScaling, freqStretching, timbrePreservation, fs)

	# time scale the sound
	yhfreq, yhmag = ST.sineTimeScaling(yhfreq, yhmag, timeScaling)

	# synthesis 
	y = SM.sineModelSynth(yhfreq, yhmag, np.array([]), Ns, H, fs)

	# write output sound 
	outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_harmonicModelTransformation.wav'
	UF.wavwrite(y,fs, outputFile)

	# --------- plotting --------------------

	# create figure to plot
	plt.figure(figsize=(12, 6))

	# frequency range to plot
	maxplotfreq = 15000.0

	plt.subplot(2,1,1)
	# plot the transformed sinusoidal frequencies
	tracks = yhfreq*np.less(yhfreq, maxplotfreq)
	tracks[tracks<=0] = np.nan
	numFrames = int(tracks[:,0].size)
	frmTime = H*np.arange(numFrames)/float(fs)
	plt.plot(frmTime, tracks, color='k')
	plt.title('transformed harmonic tracks')
	plt.autoscale(tight=True)

	# plot the output sound
	plt.subplot(2,1,2)
	plt.plot(np.arange(y.size)/float(fs), y)
	plt.axis([0, y.size/float(fs), min(y), max(y)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')
	plt.title('output sound: y')

	plt.tight_layout()
	plt.show()

def transformation_synthesis(inputFile, fs, tfreq, tmag, freqScaling = np.array([0, 2.0, 1, .3]), 
	timeScaling = np.array([0, .0, .671, .671, 1.978, 1.978+1.0])):
	"""
	Transform the analysis values returned by the analysis function and synthesize the sound
	inputFile: name of input file; fs: sampling rate of input file	
	tfreq, tmag: sinusoidal frequencies and magnitudes
	freqScaling: frequency scaling factors, in time-value pairs
	timeScaling: time scaling factors, in time-value pairs
	"""

	# size of fft used in synthesis
	Ns = 512

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

	# frequency scaling of the sinusoidal tracks 
	ytfreq = ST.sineFreqScaling(tfreq, freqScaling)

	# time scale the sinusoidal tracks 
	ytfreq, ytmag = ST.sineTimeScaling(ytfreq, tmag, timeScaling)

	# synthesis 
	y = SM.sineModelSynth(ytfreq, ytmag, np.array([]), Ns, H, fs)

	# write output sound 
	outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_sineModelTransformation.wav'
	UF.wavwrite(y,fs, outputFile)

	# create figure to plot
	plt.figure(figsize=(12, 6))

	# frequency range to plot
	maxplotfreq = 15000.0

	# plot the transformed sinusoidal frequencies
	if (ytfreq.shape[1] > 0):
		plt.subplot(2,1,1)
		tracks = np.copy(ytfreq)
		tracks = tracks*np.less(tracks, maxplotfreq)
		tracks[tracks<=0] = np.nan
		numFrames = int(tracks[:,0].size)
		frmTime = H*np.arange(numFrames)/float(fs)
		plt.plot(frmTime, tracks)
		plt.title('transformed sinusoidal tracks')
		plt.autoscale(tight=True)

	# plot the output sound
	plt.subplot(2,1,2)
	plt.plot(np.arange(y.size)/float(fs), y)
	plt.axis([0, y.size/float(fs), min(y), max(y)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')
	plt.title('output sound: y')

	plt.tight_layout()
	plt.show()

def main(inputFile='../../sounds/ocean.wav', H=256, stocf=.1):

	# ------- analysis parameters -------------------

	# inputFile: input sound file (monophonic with sampling rate of 44100)
	# H: hop size
	# stocf: decimation factor used for the stochastic approximation

	# --------- computation -----------------  

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

	# compute stochastic model                                          
	mYst = STM.stochasticModelAnal(x, H, stocf)             

	# synthesize sound from stochastic model
	y = STM.stochasticModelSynth(mYst, H)    

	outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_stochasticModel.wav'

	# write output sound
	UF.wavwrite(y, fs, outputFile)               

	# --------- plotting --------------------

	# create figure to plot
	plt.figure(figsize=(12, 9))

	# plot the input sound
	plt.subplot(3,1,1)
	plt.plot(np.arange(x.size)/float(fs), x)
	plt.axis([0, x.size/float(fs), min(x), max(x)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')
	plt.title('input sound: x')

	# plot stochastic representation
	plt.subplot(3,1,2)
	numFrames = int(mYst[:,0].size)
	frmTime = H*np.arange(numFrames)/float(fs)                             
	binFreq = np.arange(stocf*H)*float(fs)/(stocf*2*H)                      
	plt.pcolormesh(frmTime, binFreq, np.transpose(mYst))
	plt.autoscale(tight=True)
	plt.xlabel('time (sec)')
	plt.ylabel('frequency (Hz)')
	plt.title('stochastic approximation')

	# plot the output sound
	plt.subplot(3,1,3)
	plt.plot(np.arange(y.size)/float(fs), y)
	plt.axis([0, y.size/float(fs), min(y), max(y)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')

	plt.tight_layout()
	plt.show()

def downsampleAudio(inputFile, M):
    """
    Inputs:
        inputFile: file name of the wav file (including path)
        	M: downsampling factor (positive integer)
    """
    ## Your code here
    fs, x = wavread(inputFile)
    y = hopSamples(x, M)
    wavwrite(y, fs, 'test.wav')

def main(inputFile='../../sounds/ocean.wav', H=256, N=512, stocf=.1):
	"""
	inputFile: input sound file (monophonic with sampling rate of 44100)
	H: hop size, N: fft size
	stocf: decimation factor used for the stochastic approximation (bigger than 0, maximum 1)
	"""

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

	# compute stochastic model                                          
	stocEnv = STM.stochasticModelAnal(x, H, N, stocf)             

	# synthesize sound from stochastic model
	y = STM.stochasticModelSynth(stocEnv, H, N)    

	outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_stochasticModel.wav'

	# write output sound
	UF.wavwrite(y, fs, outputFile)               

	# create figure to plot
	plt.figure(figsize=(12, 9))

	# plot the input sound
	plt.subplot(3,1,1)
	plt.plot(np.arange(x.size)/float(fs), x)
	plt.axis([0, x.size/float(fs), min(x), max(x)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')
	plt.title('input sound: x')

	# plot stochastic representation
	plt.subplot(3,1,2)
	numFrames = int(stocEnv[:,0].size)
	frmTime = H*np.arange(numFrames)/float(fs)                             
	binFreq = np.arange(stocf*(N/2+1))*float(fs)/(stocf*N)                      
	plt.pcolormesh(frmTime, binFreq, np.transpose(stocEnv))
	plt.autoscale(tight=True)
	plt.xlabel('time (sec)')
	plt.ylabel('frequency (Hz)')
	plt.title('stochastic approximation')

	# plot the output sound
	plt.subplot(3,1,3)
	plt.plot(np.arange(y.size)/float(fs), y)
	plt.axis([0, y.size/float(fs), min(y), max(y)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')

	plt.tight_layout()
	plt.show(block=False)

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 main(inputFile='../../sounds/ocean.wav', H=256, N=512, stocf=.1):
	"""
	inputFile: input sound file (monophonic with sampling rate of 44100)
	H: hop size, N: fft size
	stocf: decimation factor used for the stochastic approximation (bigger than 0, maximum 1)
	"""

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

	# compute stochastic model                                          
	stocEnv = STM.stochasticModelAnal(x, H, N, stocf)             

	# synthesize sound from stochastic model
	y = STM.stochasticModelSynth(stocEnv, H, N)    

	outputFile = 'output_sounds/' + os.path.basename(inputFile)[:-4] + '_stochasticModel.wav'

	# write output sound
	UF.wavwrite(y, fs, outputFile)      
	return x, fs, stocEnv, y

def main(inputFile='../../sounds/vignesh.wav', window='blackman', M=1201, N=2048, t=-90, 
	minSineDur=0.1, nH=100, minf0=130, maxf0=300, f0et=7, harmDevSlope=0.01):
	"""
	Analysis and synthesis using the harmonic 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 could 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)

	# detect harmonics of input sound
	hfreq, hmag, hphase = HM.harmonicModelAnal(x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope, minSineDur)

	# synthesize the harmonics
	y = SM.sineModelSynth(hfreq, hmag, hphase, Ns, H, fs)  

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

	# write the sound resulting from harmonic analysis
	UF.wavwrite(y, fs, outputFile)
	return x,fs,hfreq,y

def main(inputFile='../../sounds/bendir.wav', window='hamming', M=2001, N=2048, t=-80, minSineDur=0.02, 
					maxnSines=150, freqDevOffset=10, freqDevSlope=0.001):
	"""
	Perform analysis/synthesis using the sinusoidal 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
	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)

	# analyze the sound with the sinusoidal model
	tfreq, tmag, tphase = SM.sineModelAnal(x, fs, w, N, H, t, maxnSines, minSineDur, freqDevOffset, freqDevSlope)

	# synthesize the output sound from the sinusoidal representation
	y = SM.sineModelSynth(tfreq, tmag, tphase, Ns, H, fs)

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

	# write the synthesized sound obtained from the sinusoidal synthesis
	UF.wavwrite(y, fs, outputFile)
	return x,fs,tfreq,y

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 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, stocf=0.1):
	"""
	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
	stocf: decimation factor used for the stochastic approximation
	"""

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

	# compute the harmonic plus stochastic model of the whole sound
	hfreq, hmag, hphase, stocEnv = HPS.hpsModelAnal(x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope, minSineDur, Ns, stocf)
		
	# synthesize a sound from the harmonic plus stochastic representation
	y, yh, yst = HPS.hpsModelSynth(hfreq, hmag, hphase, stocEnv, Ns, H, fs)

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

	# write sounds files for harmonics, stochastic, and the sum
	UF.wavwrite(yh, fs, outputFileSines)
	UF.wavwrite(yst, fs, outputFileStochastic)
	UF.wavwrite(y, fs, outputFile)
	return x, fs, hfreq, stocEnv, y

f0et=5
harmDevSlope=0.01
stocf=0.1

Ns = 512
H = 128

(fs, x) = UF.wavread(inputFile)
w = get_window(window, M)
hfreq, hmag, hphase, mYst = HPS.hpsModelAnal(x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope, minSineDur, Ns, stocf)
timeScaling = np.array([0, 0, 2.138, 2.138-1.5, 3.146, 3.146])
yhfreq, yhmag, ystocEnv = HPST.hpsTimeScale(hfreq, hmag, mYst, timeScaling)

y, yh, yst = HPS.hpsModelSynth(yhfreq, yhmag, np.array([]), ystocEnv, Ns, H, fs)

UF.wavwrite(y,fs, 'hps-transformation.wav')


plt.figure(figsize=(12, 9))

maxplotfreq = 14900.0

# plot the input sound
plt.subplot(4,1,1)
plt.plot(np.arange(x.size)/float(fs), x)
plt.axis([0, x.size/float(fs), min(x), max(x)])
plt.title('x (sax-phrase-short.wav')

# plot spectrogram stochastic compoment
plt.subplot(4,1,2)
numFrames = int(mYst[:,0].size)

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

	# create figure to plot
	plt.figure(figsize=(12, 9))

	# frequency range to plot
	maxplotfreq = 5000.0

	# plot the input sound
	plt.subplot(3,1,1)
	plt.plot(np.arange(x.size)/float(fs), x)
	plt.axis([0, x.size/float(fs), min(x), max(x)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')
	plt.title('input sound: x')

	# plot the magnitude spectrogram of residual
	plt.subplot(3,1,2)
	maxplotbin = int(N*maxplotfreq/fs)
	numFrames = int(mXr[:,0].size)
	frmTime = H*np.arange(numFrames)/float(fs)
	binFreq = np.arange(maxplotbin+1)*float(fs)/N
	plt.pcolormesh(frmTime, binFreq, np.transpose(mXr[:,:maxplotbin+1]))
	plt.autoscale(tight=True)

	# plot harmonic frequencies on residual spectrogram
	if (hfreq.shape[1] > 0):
		harms = hfreq*np.less(hfreq,maxplotfreq)
		harms[harms==0] = np.nan
		numFrames = int(harms[:,0].size)
		frmTime = H*np.arange(numFrames)/float(fs)
		plt.plot(frmTime, harms, color='k', ms=3, alpha=1)
		plt.xlabel('time(s)')
		plt.ylabel('frequency(Hz)')
		plt.autoscale(tight=True)
		plt.title('harmonics + residual spectrogram')

	# plot the output sound
	plt.subplot(3,1,3)
	plt.plot(np.arange(y.size)/float(fs), y)
	plt.axis([0, y.size/float(fs), min(y), max(y)])
	plt.ylabel('amplitude')
	plt.xlabel('time (sec)')
	plt.title('output sound: y')

	plt.tight_layout()
	plt.ion()
	plt.show()

binFreq = np.arange(maxplotbin+1)*float(fs)/N                         
plt.pcolormesh(frmTime, binFreq, np.transpose(mX[:,:maxplotbin+1]))
plt.autoscale(tight=True)

plt.subplot(4,1,3)
numFrames = int(ytfreq[:,0].size)
frmTime = H*np.arange(numFrames)/float(fs)
tracks = ytfreq*np.less(ytfreq, maxplotfreq)
tracks[tracks<=0] = np.nan
plt.plot(frmTime, tracks, color='k', lw=1)
plt.autoscale(tight=True)
plt.title('mY + time-scaled sine frequencies') 

maxplotbin = int(N*maxplotfreq/fs)
numFrames = int(mY[:,0].size)
frmTime = H*np.arange(numFrames)/float(fs)                             
binFreq = np.arange(maxplotbin+1)*float(fs)/N                         
plt.pcolormesh(frmTime, binFreq, np.transpose(mY[:,:maxplotbin+1]))
plt.autoscale(tight=True) 

plt.subplot(4,1,4)
plt.plot(np.arange(y.size)/float(fs), y, 'b')
plt.axis([0,y.size/float(fs),min(y),max(y)])
plt.title('y')    

plt.tight_layout()
UF.wavwrite(y, fs, 'mridangam-sineModelTimeScale.wav')
plt.savefig('sineModelTimeScale-mridangam.png')
plt.show()

def exploreSineModelMultiRes(inputFile='../../sounds/orchestra.wav'):
    """
    inputFile (string) = wav file including the path
    """
    fs, x = UF.wavread(inputFile)               # read input sound

    # First, let's check whether the new code returns same result as old one for mono-resolution case    
    
    verifySineModelMultiRes()

   # Let's find optimal parameters in a reasonable range 
  
    windows =['hanning', 'hamming', 'blackman', 'blackmanharris']

    best = Best()
    
    for k in range(5,80,5):
        m = k * 100 + 1                                  # Window size in samples
        for window in windows:                           # Window type
            for t in range(-90,-100,-10):                # Threshold
                for Ns in [512]:                         # size of fft used in synthesis
                    n = 2
                    while n < m: n = n * 2                           # size of fft used in analysis                   
                    for nPower in range(0,3):                        # try out the analysis window closest to window size, and some larger ones 
                        for nAdditionalResolutions in range(0,4):    # try out multi-resolution analysis windows
                            W = np.array([window])
                            M = np.array([m])
                            N = np.array([n])
                            B = np.array([ ])
                            T = np.array([t])

                            log_m = np.log(float(m))
                            log_n = np.log(float(n))
                            log_f = np.log(fs/2.0) 
                            log_step = np.log(2)
                            
                            executeStep = True
                            continueAddingResolutions = True
                            for additionalResolution in range(0,nAdditionalResolutions):
                                if continueAddingResolutions:
                                    scaledM = int(np.exp(log_m - log_step*(additionalResolution+1)))
                                    if scaledM % 2 == 0: scaledM = scaledM + 1
                                    scaledN = int(np.exp(log_n - log_step*(additionalResolution+1)))
                                    if scaledN < scaledM: scaledN = scaledM
                                    appropriateScaledN = 2
                                    while appropriateScaledN < scaledN: appropriateScaledN = appropriateScaledN * 2
                                    frequencyBoundary = np.exp(log_f - (log_step*(nAdditionalResolutions - additionalResolution)))
                                    if scaledM < Ns:
                                        continueAddingResolutions = False
                                        if additionalResolution == 0: executeStep = False
                                    else:
                                        W = np.append(W,window)
                                        M = np.append(M,scaledM)
                                        N = np.append(N,appropriateScaledN)
                                        B = np.append(B,frequencyBoundary)
                                        T = np.append(T,t)
                            if executeStep:
                                best.calculateAndUpdate(x, fs, Ns, W, M, N, B, T)
                        n = n * 2
                        
    print 'FILE:',inputFile
    print 'BEST:','diff =',best.diff,'for W =',best.W,', M =',best.M,', N =',best.N,', B =',best.B,', T =',best.T,', Ns =',best.Ns
    
    y_best = best.calculateAndUpdate(x, fs, best.Ns, best.W, best.M, best.N, best.B, best.T)
    outputFile = inputFile[:-4] + '_optimizedSineModel.wav'
    UF.wavwrite(y_best, fs, outputFile)

plt.subplot(311)
numFrames = int(mX[:,0].size)
frmTime = H1*np.arange(numFrames)/float(fs)                             
binFreq = fs*np.arange(N1*maxplotfreq/fs)/N1                       
plt.pcolormesh(frmTime, binFreq, np.transpose(mX[:,:int(N1*maxplotfreq/fs+1)])) 
plt.title('mX (orchestra.wav)')
plt.autoscale(tight=True)

plt.subplot(312)
numFrames = int(mX2[:,0].size)
frmTime = H1*np.arange(numFrames)/float(fs)  
                 
N = 2*mX2[0,:].size         
binFreq = fs*np.arange(N*maxplotfreq/fs)/N                       
plt.pcolormesh(frmTime, binFreq, np.transpose(mX2[:,:int(N*maxplotfreq/fs+1)]))
plt.title('mX2 (speech-male.wav)')
plt.autoscale(tight=True)

plt.subplot(313)
numFrames = int(mY[:,0].size)
frmTime = H1*np.arange(numFrames)/float(fs)                             
binFreq = fs*np.arange(N1*maxplotfreq/fs)/N1                       
plt.pcolormesh(frmTime, binFreq, np.transpose(mY[:,:int(N1*maxplotfreq/fs+1)])) 
plt.title('mY')
plt.autoscale(tight=True)

plt.tight_layout()
UF.wavwrite(y, fs, 'orchestra-speech-stftMorph.wav')
plt.savefig('stftMorph-orchestra.png')
plt.show()