diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..c10666e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+.idea
+*.pyc
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/analyzeMovieSound.py b/analyzeMovieSound.py
new file mode 100644
index 0000000..9005792
--- /dev/null
+++ b/analyzeMovieSound.py
@@ -0,0 +1,195 @@
+import os, sys, shutil, glob, numpy, csv, cPickle
+import scipy.io.wavfile as wavfile
+import audioBasicIO
+import audioTrainTest as aT
+import audioSegmentation as aS
+import matplotlib.pyplot as plt
+import scipy.spatial.distance
+minDuration = 7;
+
+def classifyFolderWrapper(inputFolder, modelType, modelName, outputMode=False):
+	if not os.path.isfile(modelName):
+		raise Exception("Input modelName not found!")
+
+	if modelType=='svm':
+		[Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadSVModel(modelName)
+	elif modelType=='knn':
+		[Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadKNNModel(modelName)
+
+	PsAll = numpy.zeros((len(classNames), ))	
+		
+	files = "*.wav"
+	if os.path.isdir(inputFolder):
+		strFilePattern = os.path.join(inputFolder, files)
+	else:
+		strFilePattern = inputFolder + files
+
+	wavFilesList = []
+	wavFilesList.extend(glob.glob(strFilePattern))
+	wavFilesList = sorted(wavFilesList)
+	if len(wavFilesList)==0:
+		print "No WAV files found!"
+		return 
+	
+	Results = []
+	for wavFile in wavFilesList:	
+		[Fs, x] = audioBasicIO.readAudioFile(wavFile)	
+		signalLength = x.shape[0] / float(Fs)
+		[Result, P, classNames] = aT.fileClassification(wavFile, modelName, modelType)					
+		PsAll += (numpy.array(P) * signalLength)		
+		Result = int(Result)
+		Results.append(Result)
+		if outputMode:
+			print "{0:s}\t{1:s}".format(wavFile,classNames[Result])
+	Results = numpy.array(Results)
+	
+	# print distribution of classes:
+	[Histogram, _] = numpy.histogram(Results, bins=numpy.arange(len(classNames)+1))
+	if outputMode:	
+		for i,h in enumerate(Histogram):
+			print "{0:20s}\t\t{1:d}".format(classNames[i], h)
+	PsAll = PsAll / numpy.sum(PsAll)
+
+
+	if outputMode:	
+		fig = plt.figure()
+		ax = fig.add_subplot(111)
+		plt.title("Classes percentage " + inputFolder.replace('Segments',''))
+		ax.axis((0, len(classNames)+1, 0, 1))
+		ax.set_xticks(numpy.array(range(len(classNames)+1)))
+		ax.set_xticklabels([" "] + classNames)
+		ax.bar(numpy.array(range(len(classNames)))+0.5, PsAll)
+		plt.show()
+	return classNames, PsAll
+
+def getMusicSegmentsFromFile(inputFile):	
+	modelType = "svm"
+	modelName = "data/svmMovies8classes"
+	
+	dirOutput = inputFile[0:-4] + "_musicSegments"
+	
+	if os.path.exists(dirOutput) and dirOutput!=".":
+		shutil.rmtree(dirOutput)	
+	os.makedirs(dirOutput)	
+	
+	[Fs, x] = audioBasicIO.readAudioFile(inputFile)	
+
+	if modelType=='svm':
+		[Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadSVModel(modelName)
+	elif modelType=='knn':
+		[Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadKNNModel(modelName)
+
+	flagsInd, classNames, acc, CM = aS.mtFileClassification(inputFile, modelName, modelType, plotResults = False, gtFile = "")
+	segs, classes = aS.flags2segs(flagsInd, mtStep)
+
+	for i, s in enumerate(segs):
+		if (classNames[int(classes[i])] == "Music") and (s[1] - s[0] >= minDuration):
+			strOut = "{0:s}{1:.3f}-{2:.3f}.wav".format(dirOutput+os.sep, s[0], s[1])	
+			wavfile.write( strOut, Fs, x[int(Fs*s[0]):int(Fs*s[1])])
+
+def analyzeDir(dirPath):
+	for i,f in enumerate(glob.glob(dirPath + os.sep + '*.wav')):				# for each WAV file					
+		getMusicSegmentsFromFile(f)	
+		[c, P]= classifyFolderWrapper(f[0:-4] + "_musicSegments", "svm", "data/svmMusicGenre8", False)
+		if i==0:
+			print "".ljust(100)+"\t",
+			for C in c:
+				print C.ljust(12)+"\t",
+			print
+		print f.ljust(100)+"\t",
+		for p in P:
+				print "{0:.2f}".format(p).ljust(12)+"\t",
+		print
+		
+def main(argv):	
+	
+	if argv[1]=="--file":
+		getMusicSegmentsFromFile(argv[2])	
+		classifyFolderWrapper(argv[2][0:-4] + "_musicSegments", "svm", "data/svmMusicGenre8", True)		
+		
+	elif argv[1]=="--dir":	
+		analyzeDir(argv[2])	
+		
+	elif argv[1]=="--sim":
+		csvFile = argv[2]
+		f = []
+		fileNames = []
+		with open(csvFile, 'rb') as csvfile:
+			spamreader = csv.reader(csvfile, delimiter='\t', quotechar='|')
+			for j,row in enumerate(spamreader):
+				if j>0:
+					ftemp = []
+					for i in range(1,9):
+						ftemp.append(float(row[i]))
+					f.append(ftemp)
+					R = row[0]
+					II = R.find(".wav");
+					fileNames.append(row[0][0:II])
+			f = numpy.array(f)
+
+			Sim = numpy.zeros((f.shape[0], f.shape[0]))
+			for i in range(f.shape[0]):
+				for j in range(f.shape[0]):	
+					Sim[i,j] = scipy.spatial.distance.cdist(numpy.reshape(f[i,:], (f.shape[1],1)).T, numpy.reshape(f[j,:], (f.shape[1],1)).T, 'cosine')
+								
+			Sim1 = numpy.reshape(Sim, (Sim.shape[0]*Sim.shape[1], 1))
+			plt.hist(Sim1)
+			plt.show()
+
+			fo = open(csvFile + "_simMatrix", "wb")
+			cPickle.dump(fileNames,  fo, protocol = cPickle.HIGHEST_PROTOCOL)
+			cPickle.dump(f, fo, protocol = cPickle.HIGHEST_PROTOCOL)			
+			cPickle.dump(Sim, fo, protocol = cPickle.HIGHEST_PROTOCOL)
+			fo.close()
+
+	elif argv[1]=="--loadsim":
+		try:
+			fo = open(argv[2], "rb")
+		except IOError:
+				print "didn't find file"
+				return
+		try:			
+			fileNames 	= cPickle.load(fo)
+			f 			= cPickle.load(fo)
+			Sim 		= cPickle.load(fo)
+		except:
+			fo.close()
+		fo.close()	
+		print fileNames
+		Sim1 = numpy.reshape(Sim, (Sim.shape[0]*Sim.shape[1], 1))
+		plt.hist(Sim1)
+		plt.show()
+
+	elif argv[1]=="--audio-event-dir":		
+		files = "*.wav"
+		inputFolder = argv[2]
+		if os.path.isdir(inputFolder):
+			strFilePattern = os.path.join(inputFolder, files)
+		else:
+			strFilePattern = inputFolder + files
+
+		wavFilesList = []
+		wavFilesList.extend(glob.glob(strFilePattern))
+		wavFilesList = sorted(wavFilesList)		
+		for i,w in enumerate(wavFilesList):			
+			[flagsInd, classesAll, acc, CM] = aS.mtFileClassification(w, "data/svmMovies8classes", "svm", False, '')
+			histTemp = numpy.zeros( (len(classesAll), ) )
+			for f in flagsInd:
+				histTemp[int(f)] += 1.0
+			histTemp /= histTemp.sum()
+			
+			if i==0:
+				print "".ljust(100)+"\t",
+				for C in classesAll:
+					print C.ljust(12)+"\t",
+				print
+			print w.ljust(100)+"\t",
+			for h in histTemp:				
+				print "{0:.2f}".format(h).ljust(12)+"\t",
+			print
+
+			
+	return 0
+	
+if __name__ == '__main__':
+	main(sys.argv)
diff --git a/audacityAnnotation2WAVs.py b/audacityAnnotation2WAVs.py
new file mode 100644
index 0000000..dd456f9
--- /dev/null
+++ b/audacityAnnotation2WAVs.py
@@ -0,0 +1,67 @@
+import glob
+import os
+import audioBasicIO
+import sys
+import csv
+import scipy.io.wavfile as wavfile
+
+
+def annotation2files(wavFile, csvFile):
+    '''
+        Break an audio stream to segments of interest, 
+        defined by a csv file
+        
+        - wavFile:    path to input wavfile
+        - csvFile:    path to csvFile of segment limits
+        
+        Input CSV file must be of the format <T1>\t<T2>\t<Label>
+    '''    
+    
+    [Fs, x] = audioBasicIO.readAudioFile(wavFile)
+    with open(csvFile, 'rb') as csvfile:
+        reader = csv.reader(csvfile, delimiter='\t', quotechar='|')
+        for j, row in enumerate(reader):
+            T1 = float(row[0].replace(",","."))
+            T2 = float(row[1].replace(",","."))            
+            label = "%s_%s_%.2f_%.2f.wav" % (wavFile, row[2], T1, T2)
+            label = label.replace(" ", "_")
+            xtemp = x[int(round(T1*Fs)):int(round(T2*Fs))]            
+            print T1, T2, label, xtemp.shape
+            wavfile.write(label, Fs, xtemp)  
+
+def main(argv):
+    if argv[1] == "-f":
+        wavFile = argv[2]
+        annotationFile = argv[3]
+        annotation2files(wavFile, annotationFile)
+    elif argv[1] == "-d":
+        inputFolder = argv[2]
+        types = ('*.txt', '*.csv')
+        annotationFilesList = []
+        for files in types:
+            annotationFilesList.extend(glob.glob(os.path.join(inputFolder, files)))
+        for anFile in annotationFilesList:
+            wavFile = os.path.splitext(anFile)[0] + ".wav"
+            if not os.path.isfile(wavFile):
+                wavFile = os.path.splitext(anFile)[0] + ".mp3"
+                if not os.path.isfile(wavFile):
+                    print "Audio file not found!"
+                    return
+            annotation2files(wavFile, anFile)
+
+
+if __name__ == '__main__':
+    # Used to extract a series of annotated WAV files based on (a) an audio file (mp3 or wav) and 
+    # (b) a segment annotation file e.g. a "label" file generated in audacity
+    #
+    # usage 1:
+    # python audacityAnnotation2WAVs.py -f <audiofilepath> <annotationfilepath>
+    # The <annotationfilepath> is actually a tab-seperated file where each line has the format <startTime>\t<entTime>\t<classLabel>
+    # The result of this process is a  series of WAV files with a file name <audiofilepath>_<startTime>_<endTime>_<classLabel>
+    # 
+    # usage 2:
+    # python audacityAnnotation2WAVs.py -d <annotationfolderpath>
+    # Same but searches all .txt and .csv annotation files. Audio files are supposed to be in the same path / filename with a WAV extension
+
+    main(sys.argv)
+    
\ No newline at end of file
diff --git a/audioAnalysis.py b/audioAnalysis.py
new file mode 100644
index 0000000..7e1cc7d
--- /dev/null
+++ b/audioAnalysis.py
@@ -0,0 +1,628 @@
+#!/usr/bin/env python2.7
+import argparse
+import os
+import numpy
+import glob
+import matplotlib.pyplot as plt
+import audioFeatureExtraction as aF
+import audioTrainTest as aT
+import audioSegmentation as aS
+import audioVisualization as aV
+import audioBasicIO
+import scipy.io.wavfile as wavfile
+import matplotlib.patches
+
+
+def dirMp3toWavWrapper(directory, samplerate, channels):
+    if not os.path.isdir(directory):
+        raise Exception("Input path not found!")
+
+    useMp3TagsAsNames = True
+    audioBasicIO.convertDirMP3ToWav(directory, samplerate, channels,
+                                    useMp3TagsAsNames)
+
+
+def dirWAVChangeFs(directory, samplerate, channels):
+    if not os.path.isdir(directory):
+        raise Exception("Input path not found!")
+
+    audioBasicIO.convertFsDirWavToWav(directory, samplerate, channels)
+
+
+def featureExtractionFileWrapper(wavFileName, outFile, mtWin, mtStep,
+                                 stWin, stStep):
+    if not os.path.isfile(wavFileName):
+        raise Exception("Input audio file not found!")
+
+    aF.mtFeatureExtractionToFile(wavFileName, mtWin, mtStep, stWin,
+                                 stStep, outFile, True, True, True)
+
+
+def beatExtractionWrapper(wavFileName, plot):
+    if not os.path.isfile(wavFileName):
+        raise Exception("Input audio file not found!")
+    [Fs, x] = audioBasicIO.readAudioFile(wavFileName)
+    F = aF.stFeatureExtraction(x, Fs, 0.050 * Fs, 0.050 * Fs)
+    BPM, ratio = aF.beatExtraction(F, 0.050, plot)
+    print "Beat: {0:d} bpm ".format(int(BPM))
+    print "Ratio: {0:.2f} ".format(ratio)
+
+
+def featureExtractionDirWrapper(directory, mtWin, mtStep, stWin, stStep):
+    if not os.path.isdir(directory):
+        raise Exception("Input path not found!")
+    aF.mtFeatureExtractionToFileDir(directory, mtWin, mtStep, stWin,
+                                    stStep, True, True, True)
+
+
+def featureVisualizationDirWrapper(directory):
+    if not os.path.isdir(directory):
+        raise Exception("Input folder not found!")
+    aV.visualizeFeaturesFolder(directory, "pca", "")
+    #aV.visualizeFeaturesFolder(directory, "lda", "artist")
+
+
+def fileSpectrogramWrapper(wavFileName):
+    if not os.path.isfile(wavFileName):
+        raise Exception("Input audio file not found!")
+    [Fs, x] = audioBasicIO.readAudioFile(wavFileName)
+    x = audioBasicIO.stereo2mono(x)
+    specgram, TimeAxis, FreqAxis = aF.stSpectogram(x, Fs, round(Fs * 0.040),
+                                                   round(Fs * 0.040), True)
+
+
+def fileChromagramWrapper(wavFileName):
+    if not os.path.isfile(wavFileName):
+        raise Exception("Input audio file not found!")
+    [Fs, x] = audioBasicIO.readAudioFile(wavFileName)
+    x = audioBasicIO.stereo2mono(x)
+    specgram, TimeAxis, FreqAxis = aF.stChromagram(x, Fs, round(Fs * 0.040),
+                                                   round(Fs * 0.040), True)
+
+
+def trainClassifierWrapper(method, beatFeatures, directories, modelName):
+    if len(directories) < 2:
+        raise Exception("At least 2 directories are needed")
+    aT.featureAndTrain(directories, 1, 1, aT.shortTermWindow, aT.shortTermStep,
+                       method.lower(), modelName, computeBEAT=beatFeatures)
+
+
+def trainRegressionWrapper(method, beatFeatures, dirName, modelName):
+    aT.featureAndTrainRegression(dirName, 1, 1, aT.shortTermWindow,
+                                 aT.shortTermStep, method.lower(), modelName,
+                                 computeBEAT=beatFeatures)
+
+
+def classifyFileWrapper(inputFile, modelType, modelName):
+    if not os.path.isfile(modelName):
+        raise Exception("Input modelName not found!")
+    if not os.path.isfile(inputFile):
+        raise Exception("Input audio file not found!")
+
+    [Result, P, classNames] = aT.fileClassification(inputFile, modelName,
+                                                    modelType)
+    print "{0:s}\t{1:s}".format("Class", "Probability")
+    for i, c in enumerate(classNames):
+        print "{0:s}\t{1:.2f}".format(c, P[i])
+    print "Winner class: " + classNames[int(Result)]
+
+
+def regressionFileWrapper(inputFile, modelType, modelName):
+    if not os.path.isfile(inputFile):
+        raise Exception("Input audio file not found!")
+
+    R, regressionNames = aT.fileRegression(inputFile, modelName, modelType)
+    for i in range(len(R)):
+        print "{0:s}\t{1:.3f}".format(regressionNames[i], R[i])
+
+
+def classifyFolderWrapper(inputFolder, modelType, modelName, outputMode=False):
+    if not os.path.isfile(modelName):
+        raise Exception("Input modelName not found!")
+
+    types = ('*.wav', '*.aif',  '*.aiff', '*.mp3')
+    wavFilesList = []
+    for files in types:
+        wavFilesList.extend(glob.glob(os.path.join(inputFolder, files)))
+    wavFilesList = sorted(wavFilesList)
+    if len(wavFilesList) == 0:
+        print "No WAV files found!"
+        return
+
+    Results = []
+    for wavFile in wavFilesList:
+        [Result, P, classNames] = aT.fileClassification(wavFile, modelName,
+                                                        modelType)
+        Result = int(Result)
+        Results.append(Result)
+        if outputMode:
+            print "{0:s}\t{1:s}".format(wavFile, classNames[Result])
+    Results = numpy.array(Results)
+
+    # print distribution of classes:
+    [Histogram, _] = numpy.histogram(Results,
+                                     bins=numpy.arange(len(classNames) + 1))
+    for i, h in enumerate(Histogram):
+        print "{0:20s}\t\t{1:d}".format(classNames[i], h)
+
+
+def regressionFolderWrapper(inputFolder, modelType, modelName):
+    files = "*.wav"
+    if os.path.isdir(inputFolder):
+        strFilePattern = os.path.join(inputFolder, files)
+    else:
+        strFilePattern = inputFolder + files
+
+    wavFilesList = []
+    wavFilesList.extend(glob.glob(strFilePattern))
+    wavFilesList = sorted(wavFilesList)
+    if len(wavFilesList) == 0:
+        print "No WAV files found!"
+        return
+    Results = []
+    for wavFile in wavFilesList:
+        R, regressionNames = aT.fileRegression(wavFile, modelName, modelType)
+        Results.append(R)
+    Results = numpy.array(Results)
+
+    for i, r in enumerate(regressionNames):
+        [Histogram, bins] = numpy.histogram(Results[:, i])
+        centers = (bins[0:-1] + bins[1::]) / 2.0        
+        plt.subplot(len(regressionNames), 1, i + 1)
+        plt.plot(centers, Histogram)
+        plt.title(r)
+    plt.show()
+
+
+def trainHMMsegmenter_fromfile(wavFile, gtFile, hmmModelName, mtWin, mtStep):
+    if not os.path.isfile(wavFile):
+        print "Error: wavfile does not exist!"
+        return
+    if not os.path.isfile(gtFile):
+        print "Error: groundtruth does not exist!"
+        return
+
+    aS.trainHMM_fromFile(wavFile, gtFile, hmmModelName, mtWin, mtStep)
+
+
+def trainHMMsegmenter_fromdir(directory, hmmModelName, mtWin, mtStep):
+    if not os.path.isdir(directory):
+        raise Exception("Input folder not found!")
+    aS.trainHMM_fromDir(directory, hmmModelName, mtWin, mtStep)
+
+
+def segmentclassifyFileWrapper(inputWavFile, modelName, modelType):
+    if not os.path.isfile(modelName):
+        raise Exception("Input modelName not found!")
+    if not os.path.isfile(inputWavFile):
+        raise Exception("Input audio file not found!")    
+    gtFile = ""
+    if inputWavFile[-4::]==".wav":
+        gtFile = inputWavFile.replace(".wav", ".segments")
+    if inputWavFile[-4::]==".mp3":
+        gtFile = inputWavFile.replace(".mp3", ".segments")    
+    aS.mtFileClassification(inputWavFile, modelName, modelType, True, gtFile)
+
+
+def segmentclassifyFileWrapperHMM(wavFile, hmmModelName):
+    gtFile = wavFile.replace(".wav", ".segments")
+    aS.hmmSegmentation(wavFile, hmmModelName, PLOT=True, gtFileName=gtFile)
+
+
+def segmentationEvaluation(dirName, modelName, methodName):
+    aS.evaluateSegmentationClassificationDir(dirName, modelName, methodName)
+
+
+def silenceRemovalWrapper(inputFile, smoothingWindow, weight):
+    if not os.path.isfile(inputFile):
+        raise Exception("Input audio file not found!")
+
+    [Fs, x] = audioBasicIO.readAudioFile(inputFile)
+    segmentLimits = aS.silenceRemoval(x, Fs, 0.05, 0.05,
+                                      smoothingWindow, weight, True)
+    for i, s in enumerate(segmentLimits):
+        strOut = "{0:s}_{1:.3f}-{2:.3f}.wav".format(inputFile[0:-4], s[0], s[1])
+        wavfile.write(strOut, Fs, x[int(Fs * s[0]):int(Fs * s[1])])
+
+
+def speakerDiarizationWrapper(inputFile, numSpeakers, useLDA):
+    if useLDA:
+        aS.speakerDiarization(inputFile, numSpeakers, PLOT=True)
+    else:
+        aS.speakerDiarization(inputFile, numSpeakers, LDAdim=0, PLOT=True)
+
+
+def thumbnailWrapper(inputFile, thumbnailWrapperSize):
+    stWindow = 1.0
+    stStep = 1.0
+    if not os.path.isfile(inputFile):
+        raise Exception("Input audio file not found!")
+
+    [Fs, x] = audioBasicIO.readAudioFile(inputFile)
+    if Fs == -1:    # could not read file
+        return
+
+    [A1, A2, B1, B2, Smatrix] = aS.musicThumbnailing(x, Fs, stWindow, stStep,
+                                                     thumbnailWrapperSize)
+
+    # write thumbnailWrappers to WAV files:
+    if inputFile.endswith(".wav"):
+        thumbnailWrapperFileName1 = inputFile.replace(".wav", "_thumb1.wav")
+        thumbnailWrapperFileName2 = inputFile.replace(".wav", "_thumb2.wav")
+    if inputFile.endswith(".mp3"):
+        thumbnailWrapperFileName1 = inputFile.replace(".mp3", "_thumb1.mp3")
+        thumbnailWrapperFileName2 = inputFile.replace(".mp3", "_thumb2.mp3")        
+    wavfile.write(thumbnailWrapperFileName1, Fs, x[int(Fs * A1):int(Fs * A2)])
+    wavfile.write(thumbnailWrapperFileName2, Fs, x[int(Fs * B1):int(Fs * B2)])
+    print "1st thumbnailWrapper (stored in file {0:s}): {1:4.1f}sec" \
+          " -- {2:4.1f}sec".format(thumbnailWrapperFileName1, A1, A2)
+    print "2nd thumbnailWrapper (stored in file {0:s}): {1:4.1f}sec" \
+          " -- {2:4.1f}sec".format(thumbnailWrapperFileName2, B1, B2)
+
+    # Plot self-similarity matrix:
+    fig = plt.figure()
+    ax = fig.add_subplot(111, aspect="auto")
+    plt.imshow(Smatrix)
+    # Plot best-similarity diagonal:
+    Xcenter = (A1 / stStep + A2 / stStep) / 2.0
+    Ycenter = (B1 / stStep + B2 / stStep) / 2.0
+
+    e1 = matplotlib.patches.Ellipse((Ycenter, Xcenter),
+                                    thumbnailWrapperSize * 1.4, 3, angle=45,
+                                    linewidth=3, fill=False)
+    ax.add_patch(e1)
+
+    plt.plot([B1, Smatrix.shape[0]], [A1, A1], color="k",
+             linestyle="--", linewidth=2)
+    plt.plot([B2, Smatrix.shape[0]], [A2, A2], color="k",
+             linestyle="--", linewidth=2)
+    plt.plot([B1, B1], [A1, Smatrix.shape[0]], color="k",
+             linestyle="--", linewidth=2)
+    plt.plot([B2, B2], [A2, Smatrix.shape[0]], color="k",
+             linestyle="--", linewidth=2)
+
+    plt.xlim([0, Smatrix.shape[0]])
+    plt.ylim([Smatrix.shape[1], 0])
+
+    ax.yaxis.set_label_position("right")
+    ax.yaxis.tick_right()
+
+    plt.xlabel("frame no")
+    plt.ylabel("frame no")
+    plt.title("Self-similarity matrix")
+
+    plt.show()
+
+
+def parse_arguments():
+    parser = argparse.ArgumentParser(description="A demonstration script "
+                                                 "for pyAudioAnalysis library")
+    tasks = parser.add_subparsers(
+        title="subcommands", description="available tasks",
+        dest="task", metavar="")
+
+    dirMp3Wav = tasks.add_parser("dirMp3toWav",
+                                 help="Convert all .mp3 files in a directory "
+                                      "to .wav format")
+    dirMp3Wav.add_argument("-i", "--input", required=True, help="Input folder")
+    dirMp3Wav.add_argument("-r", "--rate", type=int,
+                           choices=[8000, 16000, 32000, 44100], required=True,
+                           help="Samplerate of generated WAV files")
+    dirMp3Wav.add_argument("-c", "--channels", type=int, choices=[1, 2],
+                           required=True,
+                           help="Audio channels of generated WAV files")
+
+    dirWavRes = tasks.add_parser("dirWavResample",
+                                 help="Change samplerate of .wav "
+                                      "files in a directory")
+    dirWavRes.add_argument("-i", "--input", required=True, help="Input folder")
+    dirWavRes.add_argument("-r", "--rate", type=int,
+                           choices=[8000, 16000, 32000, 44100], required=True,
+                           help="Samplerate of generated WAV files")
+    dirWavRes.add_argument("-c", "--channels", type=int, choices=[1, 2],
+                           required=True,
+                           help="Audio channels of generated WAV files")
+
+    featExt = tasks.add_parser("featureExtractionFile",
+                               help="Extract audio features from file")
+    featExt.add_argument("-i", "--input", required=True,
+                         help="Input audio file")
+    featExt.add_argument("-o", "--output", required=True,
+                         help="Output file")
+    featExt.add_argument("-mw", "--mtwin", type=float,
+                         required=True, help="Mid-term window size")
+    featExt.add_argument("-ms", "--mtstep", type=float,
+                         required=True, help="Mid-term window step")
+    featExt.add_argument("-sw", "--stwin", type=float,
+                         default=0.050, help="Short-term window size")
+    featExt.add_argument("-ss", "--ststep", type=float,
+                         default=0.050, help="Short-term window step")
+
+    beat = tasks.add_parser("beatExtraction",
+                            help="Compute beat features of an audio file")
+    beat.add_argument("-i", "--input", required=True, help="Input audio file")
+    beat.add_argument("--plot", action="store_true", help="Generate plot")
+
+    featExtDir = tasks.add_parser("featureExtractionDir",
+                                  help="Extract audio features "
+                                       "from files in a folder")
+    featExtDir.add_argument("-i", "--input", required=True,
+                            help="Input directory")
+    featExtDir.add_argument("-mw", "--mtwin", type=float, required=True,
+                            help="Mid-term window size")
+    featExtDir.add_argument("-ms", "--mtstep", type=float, required=True,
+                            help="Mid-term window step")
+    featExtDir.add_argument("-sw", "--stwin", type=float, default=0.050,
+                            help="Short-term window size")
+    featExtDir.add_argument("-ss", "--ststep", type=float, default=0.050,
+                            help="Short-term window step")
+
+    featVis = tasks.add_parser("featureVisualization")
+    featVis.add_argument("-i", "--input", required=True, help="Input directory")
+
+    spectro = tasks.add_parser("fileSpectrogram")
+    spectro.add_argument("-i", "--input", required=True,
+                         help="Input audio file")
+
+    chroma = tasks.add_parser("fileChromagram")
+    chroma.add_argument("-i", "--input", required=True, help="Input audio file")
+
+    trainClass = tasks.add_parser("trainClassifier",
+                                  help="Train an SVM or KNN classifier")
+    trainClass.add_argument("-i", "--input", nargs="+",
+                            required=True, help="Input directories")
+    trainClass.add_argument("--method",
+                            choices=["svm", "svm_rbf", "knn", "randomforest",
+                                     "gradientboosting","extratrees"],
+                            required=True, help="Classifier type")
+    trainClass.add_argument("--beat", action="store_true",
+                            help="Compute beat features")
+    trainClass.add_argument("-o", "--output", required=True,
+                            help="Generated classifier filename")
+
+    trainReg = tasks.add_parser("trainRegression")
+    trainReg.add_argument("-i", "--input", required=True,
+                          help="Input directory")
+    trainReg.add_argument("--method", choices=["svm", "randomforest","svm_rbf"],
+                          required=True, help="Classifier type")
+    trainReg.add_argument("--beat", action="store_true",
+                          help="Compute beat features")
+    trainReg.add_argument("-o", "--output", required=True,
+                          help="Generated classifier filename")
+
+    classFile = tasks.add_parser("classifyFile",
+                                 help="Classify a file using an "
+                                      "existing classifier")
+    classFile.add_argument("-i", "--input", required=True,
+                           help="Input audio file")
+    classFile.add_argument("--model", choices=["svm", "svm_rbf", "knn",
+                                               "randomforest",
+                                               "gradientboosting",
+                                               "extratrees"],
+                           required=True, help="Classifier type (svm or knn or"
+                                               " randomforest or "
+                                               "gradientboosting or "
+                                               "extratrees)")
+    classFile.add_argument("--classifier", required=True,
+                           help="Classifier to use (path)")
+
+    trainHMM = tasks.add_parser("trainHMMsegmenter_fromfile",
+                                help="Train an HMM from file + annotation data")
+    trainHMM.add_argument("-i", "--input", required=True,
+                          help="Input audio file")
+    trainHMM.add_argument("--ground", required=True,
+                          help="Ground truth path (segments CSV file)")
+    trainHMM.add_argument("-o", "--output", required=True,
+                          help="HMM model name (path)")
+    trainHMM.add_argument("-mw", "--mtwin", type=float, required=True,
+                          help="Mid-term window size")
+    trainHMM.add_argument("-ms", "--mtstep", type=float, required=True,
+                          help="Mid-term window step")
+
+    trainHMMDir = tasks.add_parser("trainHMMsegmenter_fromdir",
+                                   help="Train an HMM from file + annotation "
+                                        "data stored in a directory (batch)")
+    trainHMMDir.add_argument("-i", "--input", required=True,
+                             help="Input audio folder")
+    trainHMMDir.add_argument("-o", "--output", required=True,
+                             help="HMM model name (path)")
+    trainHMMDir.add_argument("-mw", "--mtwin", type=float, required=True,
+                             help="Mid-term window size")
+    trainHMMDir.add_argument("-ms", "--mtstep", type=float, required=True,
+                             help="Mid-term window step")
+
+    segmentClassifyFile = tasks.add_parser("segmentClassifyFile",
+                                           help="Segmentation - classification "
+                                                "of a WAV file given a trained "
+                                                "SVM or kNN")
+    segmentClassifyFile.add_argument("-i", "--input", required=True,
+                                     help="Input audio file")
+    segmentClassifyFile.add_argument("--model",
+                                     choices=["svm", "svm_rbf", "knn",
+                                              "randomforest","gradientboosting",
+                                              "extratrees"],
+                                     required=True, help="Model type")
+    segmentClassifyFile.add_argument("--modelName", required=True,
+                                     help="Model path")
+
+    segmentClassifyFileHMM = tasks.add_parser("segmentClassifyFileHMM",
+                                              help="Segmentation - "
+                                                   "classification of a WAV "
+                                                   "file given a trained HMM")
+    segmentClassifyFileHMM.add_argument("-i", "--input", required=True,
+                                        help="Input audio file")
+    segmentClassifyFileHMM.add_argument("--hmm", required=True,
+                                        help="HMM Model to use (path)")
+
+    segmentationEvaluation = tasks.add_parser("segmentationEvaluation", help=
+                                              "Segmentation - classification "
+                                              "evaluation for a list of WAV "
+                                              "files and CSV ground-truth "
+                                              "stored in a folder")
+    segmentationEvaluation.add_argument("-i", "--input", required=True,
+                                        help="Input audio folder")
+    segmentationEvaluation.add_argument("--model",
+                                        choices=["svm", "knn", "hmm"],
+                                        required=True, help="Model type")
+    segmentationEvaluation.add_argument("--modelName", required=True,
+                                        help="Model path")
+
+    regFile = tasks.add_parser("regressionFile")
+    regFile.add_argument("-i", "--input", required=True,
+                         help="Input audio file")
+    regFile.add_argument("--model", choices=["svm", "svm_rbf","randomforest"],
+                         required=True, help="Regression type")
+    regFile.add_argument("--regression", required=True,
+                         help="Regression model to use")
+
+    classFolder = tasks.add_parser("classifyFolder")
+    classFolder.add_argument("-i", "--input", required=True,
+                             help="Input folder")
+    classFolder.add_argument("--model", choices=["svm", "svm_rbf", "knn",
+                                                 "randomforest",
+                                                 "gradientboosting",
+                                                 "extratrees"],
+                             required=True, help="Classifier type")
+    classFolder.add_argument("--classifier", required=True,
+                             help="Classifier to use (filename)")
+    classFolder.add_argument("--details", action="store_true",
+                             help="Plot details (otherwise only "
+                                  "counts per class are shown)")
+
+    regFolder = tasks.add_parser("regressionFolder")
+    regFolder.add_argument("-i", "--input", required=True, help="Input folder")
+    regFolder.add_argument("--model", choices=["svm", "knn"],
+                           required=True, help="Classifier type")
+    regFolder.add_argument("--regression", required=True,
+                           help="Regression model to use")
+
+    silrem = tasks.add_parser("silenceRemoval",
+                              help="Remove silence segments from a recording")
+    silrem.add_argument("-i", "--input", required=True, help="input audio file")
+    silrem.add_argument("-s", "--smoothing", type=float, default=1.0,
+                        help="smoothing window size in seconds.")
+    silrem.add_argument("-w", "--weight", type=float, default=0.5,
+                        help="weight factor in (0, 1)")
+
+    spkrDir = tasks.add_parser("speakerDiarization")
+    spkrDir.add_argument("-i", "--input", required=True,
+                         help="Input audio file")
+    spkrDir.add_argument("-n", "--num", type=int, required=True,
+                         help="Number of speakers")
+    spkrDir.add_argument("--flsd", action="store_true",
+                         help="Enable FLsD method")
+
+    speakerDiarizationScriptEval = tasks.add_parser("speakerDiarizationScriptEval",
+                                                    help="Train an SVM or KNN "
+                                                         "classifier")
+    speakerDiarizationScriptEval.add_argument("-i", "--input", required=True,
+                                              help="Input directory")
+    speakerDiarizationScriptEval.add_argument("--LDAs", type=int, nargs="+",
+                                              required=True,
+                                              help="List FLsD params")
+
+    thumb = tasks.add_parser("thumbnail",
+                             help="Generate a thumbnailWrapper "
+                                  "for an audio file")
+    thumb.add_argument("-i", "--input", required=True, help="input audio file")
+    thumb.add_argument("-s", "--size",  default=10.0,  type=float,
+                       help="thumbnailWrapper size in seconds.")
+
+    return parser.parse_args()
+
+
+def trainHMMsegmenter_fromdir(directory, hmmModelName, mtWin, mtStep):
+    if not os.path.isdir(directory):
+        raise Exception("Input folder not found!")
+    aS.trainHMM_fromDir(directory, hmmModelName, mtWin, mtStep)
+
+
+def segmentclassifyFileWrapperHMM(wavFile, hmmModelName):
+    gtFile = wavFile.replace(".wav", ".segments")
+    aS.hmmSegmentation(wavFile, hmmModelName, PLOT=True, gtFileName=gtFile)
+
+
+if __name__ == "__main__":
+    args = parse_arguments()
+
+
+    if args.task == "dirMp3toWav":
+        # Convert mp3 to wav (batch - folder)
+        dirMp3toWavWrapper(args.input, args.rate, args.channels)
+    elif args.task == "dirWavResample":
+        # Convert Fs for a list of wavs stored in a folder
+        dirWAVChangeFs(args.input, args.rate, args.channels)
+    elif args.task == "featureExtractionFile":
+        # Feature extraction for WAV file
+        featureExtractionFileWrapper(args.input, args.output, args.mtwin,
+                                     args.mtstep, args.stwin, args.ststep)
+    elif args.task == "featureExtractionDir":
+        # Feature extraction for all WAV files stored in a folder
+        featureExtractionDirWrapper(args.input, args.mtwin, args.mtstep,
+                                    args.stwin, args.ststep)
+    elif args.task == "fileSpectrogram":
+        # Extract spectrogram from a WAV file
+        fileSpectrogramWrapper(args.input)
+    elif args.task == "fileChromagram":
+        # Extract chromagram from a WAV file
+        fileChromagramWrapper(args.input)
+    elif args.task == "featureVisualization":
+        # Visualize the content of a list of WAV files stored in a folder
+        featureVisualizationDirWrapper(args.input)
+    elif args.task == "beatExtraction":
+        # Extract bpm from file
+        beatExtractionWrapper(args.input, args.plot)
+    elif args.task == "trainClassifier":
+        # Train classifier from data (organized in folders)
+        trainClassifierWrapper(args.method, args.beat, args.input, args.output)
+    elif args.task == "trainRegression":
+        # Train a regression model from data (organized in
+        # a single folder, while ground-truth is provided in a CSV)
+        trainRegressionWrapper(args.method, args.beat, args.input, args.output)
+    elif args.task == "classifyFile":
+        # Apply audio classifier on audio file
+        classifyFileWrapper(args.input, args.model, args.classifier)
+    elif args.task == "trainHMMsegmenter_fromfile":
+        # Train an hmm segmenter-classifier from WAV file + annotation
+        trainHMMsegmenter_fromfile(args.input, args.ground, args.output,
+                                   args.mtwin, args.mtstep)
+    elif args.task == "trainHMMsegmenter_fromdir":
+        # Train an hmm segmenter-classifier from a list of
+        # WAVs and annotations stored in a folder
+        trainHMMsegmenter_fromdir(args.input, args.output, args.mtwin,
+                                  args.mtstep)
+    elif args.task == "segmentClassifyFile":
+        # Apply a classifier (svm or knn or randomforest or gradientboosting
+        # or extratrees) for segmentation-classificaiton to a WAV file
+        segmentclassifyFileWrapper(args.input, args.modelName, args.model)
+    elif args.task == "segmentClassifyFileHMM":
+        # Apply an hmm for segmentation-classificaiton to a WAV file
+        segmentclassifyFileWrapperHMM(args.input, args.hmm)
+    elif args.task == "segmentationEvaluation":
+        # Evaluate segmentation-classification for a list of WAV files
+        # (and ground truth CSVs) stored in a folder
+        segmentationEvaluation(args.input, args.modelName, args.model)
+    elif args.task == "regressionFile":
+        # Apply a regression model to an audio signal stored in a WAV file
+        regressionFileWrapper(args.input, args.model, args.regression)
+    elif args.task == "classifyFolder":
+        # Classify every WAV file in a given path
+        classifyFolderWrapper(args.input, args.model, args.classifier,
+                              args.details)
+    elif args.task == "regressionFolder":
+        # Apply a regression model on every WAV file in a given path
+        regressionFolderWrapper(args.input, args.model, args.regression)
+    elif args.task == "silenceRemoval":
+        # Detect non-silent segments in a WAV file and
+        # output to seperate WAV files
+        silenceRemovalWrapper(args.input, args.smoothing, args.weight)
+    elif args.task == "speakerDiarization":
+        # Perform speaker diarization on a WAV file
+        speakerDiarizationWrapper(args.input, args.num, args.flsd)
+    elif args.task == "speakerDiarizationScriptEval":
+        # Evaluate speaker diarization given a folder that contains
+        # WAV files and .segment (Groundtruth files)
+        aS.speakerDiarizationEvaluateScript(args.input, args.LDAs)
+    elif args.task == "thumbnail":
+        # Audio thumbnailing
+        thumbnailWrapper(args.input, args.size)
diff --git a/audioAnalysisRecordAlsa.py b/audioAnalysisRecordAlsa.py
new file mode 100644
index 0000000..2d5ca5b
--- /dev/null
+++ b/audioAnalysisRecordAlsa.py
@@ -0,0 +1,149 @@
+import sys, os, alsaaudio, time, audioop, numpy, glob,  scipy, subprocess, wave, cPickle, threading, shutil
+import matplotlib.pyplot as plt
+import scipy.io.wavfile as wavfile
+from scipy.fftpack import rfft
+import audioFeatureExtraction as aF	
+import audioTrainTest as aT
+import audioSegmentation as aS
+from scipy.fftpack import fft
+import matplotlib
+matplotlib.use('TkAgg')
+
+Fs = 16000
+
+def recordAudioSegments(RecordPath, BLOCKSIZE):	
+	# This function is used for recording audio segments (until ctr+c is pressed)
+	# ARGUMENTS:
+	# - RecordPath:		the path where the wav segments will be stored
+	# - BLOCKSIZE:		segment recording size (in seconds)
+	# 
+	# NOTE: filenames are based on clock() value
+	
+	print "Press Ctr+C to stop recording"
+	RecordPath += os.sep
+	d = os.path.dirname(RecordPath)
+	if os.path.exists(d) and RecordPath!=".":
+		shutil.rmtree(RecordPath)	
+	os.makedirs(RecordPath)	
+
+	inp = alsaaudio.PCM(alsaaudio.PCM_CAPTURE,alsaaudio.PCM_NONBLOCK)
+	inp.setchannels(1)
+	inp.setrate(Fs)
+	inp.setformat(alsaaudio.PCM_FORMAT_S16_LE)
+	inp.setperiodsize(512)
+	midTermBufferSize = int(Fs*BLOCKSIZE)
+	midTermBuffer = []
+	curWindow = []
+	elapsedTime = "%08.3f" % (time.time())
+	while 1:
+			l,data = inp.read()		   
+		    	if l:
+				for i in range(len(data)/2):
+					curWindow.append(audioop.getsample(data, 2, i))
+		
+				if (len(curWindow)+len(midTermBuffer)>midTermBufferSize):
+					samplesToCopyToMidBuffer = midTermBufferSize - len(midTermBuffer)
+				else:
+					samplesToCopyToMidBuffer = len(curWindow)
+
+				midTermBuffer = midTermBuffer + curWindow[0:samplesToCopyToMidBuffer];
+				del(curWindow[0:samplesToCopyToMidBuffer])
+			
+
+			if len(midTermBuffer) == midTermBufferSize:
+				# allData = allData + midTermBuffer				
+				curWavFileName = RecordPath + os.sep + str(elapsedTime) + ".wav"				
+				midTermBufferArray = numpy.int16(midTermBuffer)
+				wavfile.write(curWavFileName, Fs, midTermBufferArray)
+				print "AUDIO  OUTPUT: Saved " + curWavFileName
+				midTermBuffer = []
+				elapsedTime = "%08.3f" % (time.time())	
+	
+def recordAnalyzeAudio(duration, outputWavFile, midTermBufferSizeSec, modelName, modelType):
+	'''
+	recordAnalyzeAudio(duration, outputWavFile, midTermBufferSizeSec, modelName, modelType)
+
+	This function is used to record and analyze audio segments, in a fix window basis.
+
+	ARGUMENTS: 
+	- duration			total recording duration
+	- outputWavFile			path of the output WAV file
+	- midTermBufferSizeSec		(fix)segment length in seconds
+	- modelName			classification model name
+	- modelType			classification model type
+
+	'''
+
+	if modelType=='svm':
+		[Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadSVModel(modelName)
+	elif modelType=='knn':
+		[Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadKNNModel(modelName)
+	else:
+		Classifier = None
+
+	inp = alsaaudio.PCM(alsaaudio.PCM_CAPTURE, alsaaudio.PCM_NONBLOCK)
+	inp.setchannels(1)
+	inp.setrate(Fs)
+	inp.setformat(alsaaudio.PCM_FORMAT_S16_LE)
+	inp.setperiodsize(512)
+	midTermBufferSize = int(midTermBufferSizeSec * Fs)
+	allData = []
+	midTermBuffer = []
+	curWindow = []
+	count = 0
+
+	while len(allData)<duration*Fs:
+		# Read data from device
+		l,data = inp.read()
+	    	if l:
+			for i in range(l):
+				curWindow.append(audioop.getsample(data, 2, i))		
+			if (len(curWindow)+len(midTermBuffer)>midTermBufferSize):
+				samplesToCopyToMidBuffer = midTermBufferSize - len(midTermBuffer)
+			else:
+				samplesToCopyToMidBuffer = len(curWindow)
+			midTermBuffer = midTermBuffer + curWindow[0:samplesToCopyToMidBuffer];
+			del(curWindow[0:samplesToCopyToMidBuffer])
+		if len(midTermBuffer) == midTermBufferSize:
+			count += 1						
+			if Classifier!=None:
+				[mtFeatures, stFeatures] = aF.mtFeatureExtraction(midTermBuffer, Fs, 2.0*Fs, 2.0*Fs, 0.020*Fs, 0.020*Fs)
+				curFV = (mtFeatures[:,0] - MEAN) / STD;
+				[result, P] = aT.classifierWrapper(Classifier, modelType, curFV)
+				print classNames[int(result)]
+			allData = allData + midTermBuffer
+
+			plt.clf()
+			plt.plot(midTermBuffer)
+			plt.show(block = False)
+			plt.draw()
+
+
+			midTermBuffer = []
+
+	allDataArray = numpy.int16(allData)
+	wavfile.write(outputWavFile, Fs, allDataArray)
+
+def main(argv):
+	if argv[1] == '-recordSegments':		# record input
+		if (len(argv)==4): 			# record segments (until ctrl+c pressed)
+			recordAudioSegments(argv[2], float(argv[3]))
+		else:
+			print "Error.\nSyntax: " + argv[0] + " -recordSegments <recordingPath> <segmentDuration>"
+
+	if argv[1] == '-recordAndClassifySegments':	# record input
+		if (len(argv)==6):			# recording + audio analysis
+			duration = int(argv[2])
+			outputWavFile = argv[3]
+			modelName = argv[4]
+			modelType = argv[5]
+			if modelType not in ["svm", "knn"]:
+				raise Exception("ModelType has to be either svm or knn!")
+			if not os.path.isfile(modelName):
+				raise Exception("Input modelName not found!")
+			recordAnalyzeAudio(duration, outputWavFile, 2.0, modelName, modelType)
+		else:
+			print "Error.\nSyntax: " + argv[0] + " -recordAndClassifySegments <duration> <outputWafFile> <modelName> <modelType>"
+	
+if __name__ == '__main__':
+	main(sys.argv)
diff --git a/audioBasicIO.py b/audioBasicIO.py
new file mode 100644
index 0000000..1e4cff6
--- /dev/null
+++ b/audioBasicIO.py
@@ -0,0 +1,130 @@
+import os, glob, eyed3, ntpath, shutil, numpy
+import scipy.io.wavfile as wavfile
+import pydub
+from pydub import AudioSegment
+
+def convertDirMP3ToWav(dirName, Fs, nC, useMp3TagsAsName = False):
+    '''
+    This function converts the MP3 files stored in a folder to WAV. If required, the output names of the WAV files are based on MP3 tags, otherwise the same names are used.
+    ARGUMENTS:
+     - dirName:     the path of the folder where the MP3s are stored
+     - Fs:          the sampling rate of the generated WAV files
+     - nC:          the number of channesl of the generated WAV files
+     - useMp3TagsAsName:    True if the WAV filename is generated on MP3 tags
+    '''
+
+    types = (dirName+os.sep+'*.mp3',) # the tuple of file types
+    filesToProcess = [] 
+
+    for files in types:
+        filesToProcess.extend(glob.glob(files))     
+
+    for f in filesToProcess:
+        #tag.link(f)
+        audioFile = eyed3.load(f)               
+        if useMp3TagsAsName and audioFile.tag != None:          
+            artist = audioFile.tag.artist
+            title = audioFile.tag.title
+            if artist!=None and title!=None:
+                if len(title)>0 and len(artist)>0:
+                    wavFileName = ntpath.split(f)[0] + os.sep + artist.replace(","," ") + " --- " + title.replace(","," ") + ".wav"
+                else:
+                    wavFileName = f.replace(".mp3",".wav")  
+            else:
+                wavFileName = f.replace(".mp3",".wav")                      
+        else:
+            wavFileName = f.replace(".mp3",".wav")      
+        command = "avconv -i \"" + f + "\" -ar " +str(Fs) + " -ac " + str(nC) + " \"" + wavFileName + "\"";
+        print command
+        os.system(command.decode('unicode_escape').encode('ascii','ignore').replace("\0",""))
+
+def convertFsDirWavToWav(dirName, Fs, nC):
+    '''
+    This function converts the WAV files stored in a folder to WAV using a different sampling freq and number of channels.
+    ARGUMENTS:
+     - dirName:     the path of the folder where the WAVs are stored
+     - Fs:          the sampling rate of the generated WAV files
+     - nC:          the number of channesl of the generated WAV files
+    '''
+
+    types = (dirName+os.sep+'*.wav',) # the tuple of file types
+    filesToProcess = []
+
+    for files in types:
+        filesToProcess.extend(glob.glob(files))     
+
+    newDir = dirName + os.sep + "Fs" + str(Fs) + "_" + "NC"+str(nC)
+    if os.path.exists(newDir) and newDir!=".":
+        shutil.rmtree(newDir)   
+    os.makedirs(newDir) 
+
+    for f in filesToProcess:    
+        _, wavFileName = ntpath.split(f)    
+        command = "avconv -i \"" + f + "\" -ar " +str(Fs) + " -ac " + str(nC) + " \"" + newDir + os.sep + wavFileName + "\"";
+        print command
+        os.system(command)
+
+def readAudioFile(path):
+    '''
+    This function returns a numpy array that stores the audio samples of a specified WAV of AIFF file
+    '''
+    extension = os.path.splitext(path)[1]
+
+    try:
+        #if extension.lower() == '.wav':
+            #[Fs, x] = wavfile.read(path)
+        if extension.lower() == '.aif' or extension.lower() == '.aiff':
+            s = aifc.open(path, 'r')
+            nframes = s.getnframes()
+            strsig = s.readframes(nframes)
+            x = numpy.fromstring(strsig, numpy.short).byteswap()
+            Fs = s.getframerate()
+        elif extension.lower() == '.mp3' or extension.lower() == '.wav' or extension.lower() == '.au':            
+            try:
+                audiofile = AudioSegment.from_file(path)
+            #except pydub.exceptions.CouldntDecodeError:
+            except:
+                print "Error: file not found or other I/O error. (DECODING FAILED)"
+                return (-1,-1)                
+
+            if audiofile.sample_width==2:                
+                data = numpy.fromstring(audiofile._data, numpy.int16)
+            elif audiofile.sample_width==4:
+                data = numpy.fromstring(audiofile._data, numpy.int32)
+            else:
+                return (-1, -1)
+            Fs = audiofile.frame_rate
+            x = []
+            for chn in xrange(audiofile.channels):
+                x.append(data[chn::audiofile.channels])
+            x = numpy.array(x).T
+        else:
+            print "Error in readAudioFile(): Unknown file type!"
+            return (-1,-1)
+    except IOError: 
+        print "Error: file not found or other I/O error."
+        return (-1,-1)
+
+    if x.ndim==2:
+        if x.shape[1]==1:
+            x = x.flatten()
+
+    return (Fs, x)
+
+def stereo2mono(x):
+    '''
+    This function converts the input signal (stored in a numpy array) to MONO (if it is STEREO)
+    '''
+    if isinstance(x, int):
+        return -1
+    if x.ndim==1:
+        return x
+    elif x.ndim==2:
+        if x.shape[1]==1:
+            return x.flatten()
+        else:
+            if x.shape[1]==2:
+                return ( (x[:,1] / 2) + (x[:,0] / 2) )
+            else:
+                return -1
+
diff --git a/audioFeatureExtraction.py b/audioFeatureExtraction.py
new file mode 100644
index 0000000..d714812
--- /dev/null
+++ b/audioFeatureExtraction.py
@@ -0,0 +1,893 @@
+import sys
+import time
+import os
+import glob
+import numpy
+import cPickle
+import aifc
+import math
+from numpy import NaN, Inf, arange, isscalar, array
+from scipy.fftpack import rfft
+from scipy.fftpack import fft
+from scipy.fftpack.realtransforms import dct
+from scipy.signal import fftconvolve
+from matplotlib.mlab import find
+import matplotlib.pyplot as plt
+from scipy import linalg as la
+import audioTrainTest as aT
+import audioBasicIO
+import utilities
+from scipy.signal import lfilter, hamming
+#from scikits.talkbox import lpc
+
+reload(sys)  
+sys.setdefaultencoding('utf8')
+
+eps = 0.00000001
+
+""" Time-domain audio features """
+
+
+def stZCR(frame):
+    """Computes zero crossing rate of frame"""
+    count = len(frame)
+    countZ = numpy.sum(numpy.abs(numpy.diff(numpy.sign(frame)))) / 2
+    return (numpy.float64(countZ) / numpy.float64(count-1.0))
+
+
+def stEnergy(frame):
+    """Computes signal energy of frame"""
+    return numpy.sum(frame ** 2) / numpy.float64(len(frame))
+
+
+def stEnergyEntropy(frame, numOfShortBlocks=10):
+    """Computes entropy of energy"""
+    Eol = numpy.sum(frame ** 2)    # total frame energy
+    L = len(frame)
+    subWinLength = int(numpy.floor(L / numOfShortBlocks))
+    if L != subWinLength * numOfShortBlocks:
+            frame = frame[0:subWinLength * numOfShortBlocks]
+    # subWindows is of size [numOfShortBlocks x L]
+    subWindows = frame.reshape(subWinLength, numOfShortBlocks, order='F').copy()
+
+    # Compute normalized sub-frame energies:
+    s = numpy.sum(subWindows ** 2, axis=0) / (Eol + eps)
+
+    # Compute entropy of the normalized sub-frame energies:
+    Entropy = -numpy.sum(s * numpy.log2(s + eps))
+    return Entropy
+
+
+""" Frequency-domain audio features """
+
+
+def stSpectralCentroidAndSpread(X, fs):
+    """Computes spectral centroid of frame (given abs(FFT))"""
+    ind = (numpy.arange(1, len(X) + 1)) * (fs/(2.0 * len(X)))
+
+    Xt = X.copy()
+    Xt = Xt / Xt.max()
+    NUM = numpy.sum(ind * Xt)
+    DEN = numpy.sum(Xt) + eps
+
+    # Centroid:
+    C = (NUM / DEN)
+
+    # Spread:
+    S = numpy.sqrt(numpy.sum(((ind - C) ** 2) * Xt) / DEN)
+
+    # Normalize:
+    C = C / (fs / 2.0)
+    S = S / (fs / 2.0)
+
+    return (C, S)
+
+
+def stSpectralEntropy(X, numOfShortBlocks=10):
+    """Computes the spectral entropy"""
+    L = len(X)                         # number of frame samples
+    Eol = numpy.sum(X ** 2)            # total spectral energy
+
+    subWinLength = int(numpy.floor(L / numOfShortBlocks))   # length of sub-frame
+    if L != subWinLength * numOfShortBlocks:
+        X = X[0:subWinLength * numOfShortBlocks]
+
+    subWindows = X.reshape(subWinLength, numOfShortBlocks, order='F').copy()  # define sub-frames (using matrix reshape)
+    s = numpy.sum(subWindows ** 2, axis=0) / (Eol + eps)                      # compute spectral sub-energies
+    En = -numpy.sum(s*numpy.log2(s + eps))                                    # compute spectral entropy
+
+    return En
+
+
+def stSpectralFlux(X, Xprev):
+    """
+    Computes the spectral flux feature of the current frame
+    ARGUMENTS:
+        X:        the abs(fft) of the current frame
+        Xpre:        the abs(fft) of the previous frame
+    """
+    # compute the spectral flux as the sum of square distances:
+    sumX = numpy.sum(X + eps)
+    sumPrevX = numpy.sum(Xprev + eps)
+    F = numpy.sum((X / sumX - Xprev/sumPrevX) ** 2)
+
+    return F
+
+
+def stSpectralRollOff(X, c, fs):
+    """Computes spectral roll-off"""
+    totalEnergy = numpy.sum(X ** 2)
+    fftLength = len(X)
+    Thres = c*totalEnergy
+    # Ffind the spectral rolloff as the frequency position where the respective spectral energy is equal to c*totalEnergy
+    CumSum = numpy.cumsum(X ** 2) + eps
+    [a, ] = numpy.nonzero(CumSum > Thres)
+    if len(a) > 0:
+        mC = numpy.float64(a[0]) / (float(fftLength))
+    else:
+        mC = 0.0
+    return (mC)
+
+
+def stHarmonic(frame, fs):
+    """
+    Computes harmonic ratio and pitch
+    """
+    M = numpy.round(0.016 * fs) - 1
+    R = numpy.correlate(frame, frame, mode='full')
+
+    g = R[len(frame)-1]
+    R = R[len(frame):-1]
+
+    # estimate m0 (as the first zero crossing of R)
+    [a, ] = numpy.nonzero(numpy.diff(numpy.sign(R)))
+
+    if len(a) == 0:
+        m0 = len(R)-1
+    else:
+        m0 = a[0]
+    if M > len(R):
+        M = len(R) - 1
+
+    Gamma = numpy.zeros((M), dtype=numpy.float64)
+    CSum = numpy.cumsum(frame ** 2)
+    Gamma[m0:M] = R[m0:M] / (numpy.sqrt((g * CSum[M:m0:-1])) + eps)
+
+    ZCR = stZCR(Gamma)
+
+    if ZCR > 0.15:
+        HR = 0.0
+        f0 = 0.0
+    else:
+        if len(Gamma) == 0:
+            HR = 1.0
+            blag = 0.0
+            Gamma = numpy.zeros((M), dtype=numpy.float64)
+        else:
+            HR = numpy.max(Gamma)
+            blag = numpy.argmax(Gamma)
+
+        # Get fundamental frequency:
+        f0 = fs / (blag + eps)
+        if f0 > 5000:
+            f0 = 0.0
+        if HR < 0.1:
+            f0 = 0.0
+
+    return (HR, f0)
+
+
+def mfccInitFilterBanks(fs, nfft):
+    """
+    Computes the triangular filterbank for MFCC computation (used in the stFeatureExtraction function before the stMFCC function call)
+    This function is taken from the scikits.talkbox library (MIT Licence):
+    https://pypi.python.org/pypi/scikits.talkbox
+    """
+
+    # filter bank params:
+    lowfreq = 133.33
+    linsc = 200/3.
+    logsc = 1.0711703
+    numLinFiltTotal = 13
+    numLogFilt = 27
+
+    if fs < 8000:
+        nlogfil = 5
+
+    # Total number of filters
+    nFiltTotal = numLinFiltTotal + numLogFilt
+
+    # Compute frequency points of the triangle:
+    freqs = numpy.zeros(nFiltTotal+2)
+    freqs[:numLinFiltTotal] = lowfreq + numpy.arange(numLinFiltTotal) * linsc
+    freqs[numLinFiltTotal:] = freqs[numLinFiltTotal-1] * logsc ** numpy.arange(1, numLogFilt + 3)
+    heights = 2./(freqs[2:] - freqs[0:-2])
+
+    # Compute filterbank coeff (in fft domain, in bins)
+    fbank = numpy.zeros((nFiltTotal, nfft))
+    nfreqs = numpy.arange(nfft) / (1. * nfft) * fs
+
+    for i in range(nFiltTotal):
+        lowTrFreq = freqs[i]
+        cenTrFreq = freqs[i+1]
+        highTrFreq = freqs[i+2]
+
+        lid = numpy.arange(numpy.floor(lowTrFreq * nfft / fs) + 1, numpy.floor(cenTrFreq * nfft / fs) + 1, dtype=numpy.int)
+        lslope = heights[i] / (cenTrFreq - lowTrFreq)
+        rid = numpy.arange(numpy.floor(cenTrFreq * nfft / fs) + 1, numpy.floor(highTrFreq * nfft / fs) + 1, dtype=numpy.int)
+        rslope = heights[i] / (highTrFreq - cenTrFreq)
+        fbank[i][lid] = lslope * (nfreqs[lid] - lowTrFreq)
+        fbank[i][rid] = rslope * (highTrFreq - nfreqs[rid])
+
+    return fbank, freqs
+
+
+def stMFCC(X, fbank, nceps):
+    """
+    Computes the MFCCs of a frame, given the fft mag
+
+    ARGUMENTS:
+        X:        fft magnitude abs(FFT)
+        fbank:    filter bank (see mfccInitFilterBanks)
+    RETURN
+        ceps:     MFCCs (13 element vector)
+
+    Note:    MFCC calculation is, in general, taken from the scikits.talkbox library (MIT Licence),
+    #    with a small number of modifications to make it more compact and suitable for the pyAudioAnalysis Lib
+    """
+
+    mspec = numpy.log10(numpy.dot(X, fbank.T)+eps)
+    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:nceps]
+    return ceps
+
+
+def stChromaFeaturesInit(nfft, fs):
+    """
+    This function initializes the chroma matrices used in the calculation of the chroma features
+    """
+    freqs = numpy.array([((f + 1) * fs) / (2 * nfft) for f in range(nfft)])    
+    Cp = 27.50    
+    nChroma = numpy.round(12.0 * numpy.log2(freqs / Cp)).astype(int)
+
+    nFreqsPerChroma = numpy.zeros((nChroma.shape[0], ))
+
+    uChroma = numpy.unique(nChroma)
+    for u in uChroma:
+        idx = numpy.nonzero(nChroma == u)
+        nFreqsPerChroma[idx] = idx[0].shape
+    
+    return nChroma, nFreqsPerChroma
+
+
+def stChromaFeatures(X, fs, nChroma, nFreqsPerChroma):
+    #TODO: 1 complexity
+    #TODO: 2 bug with large windows
+
+    chromaNames = ['A', 'A#', 'B', 'C', 'C#', 'D', 'D#', 'E', 'F', 'F#', 'G', 'G#']
+    spec = X**2    
+    if nChroma.max()<nChroma.shape[0]:        
+        C = numpy.zeros((nChroma.shape[0],))
+        C[nChroma] = spec
+        C /= nFreqsPerChroma[nChroma]
+    else:        
+        I = numpy.nonzero(nChroma>nChroma.shape[0])[0][0]        
+        C = numpy.zeros((nChroma.shape[0],))
+        C[nChroma[0:I-1]] = spec            
+        C /= nFreqsPerChroma
+    finalC = numpy.zeros((12, 1))
+    newD = int(numpy.ceil(C.shape[0] / 12.0) * 12)
+    C2 = numpy.zeros((newD, ))
+    C2[0:C.shape[0]] = C
+    C2 = C2.reshape(C2.shape[0]/12, 12)
+    #for i in range(12):
+    #    finalC[i] = numpy.sum(C[i:C.shape[0]:12])
+    finalC = numpy.matrix(numpy.sum(C2, axis=0)).T
+    finalC /= spec.sum()
+
+#    ax = plt.gca()
+#    plt.hold(False)
+#    plt.plot(finalC)
+#    ax.set_xticks(range(len(chromaNames)))
+#    ax.set_xticklabels(chromaNames)
+#    xaxis = numpy.arange(0, 0.02, 0.01);
+#    ax.set_yticks(range(len(xaxis)))
+#    ax.set_yticklabels(xaxis)
+#    plt.show(block=False)
+#    plt.draw()
+
+    return chromaNames, finalC
+
+
+def stChromagram(signal, Fs, Win, Step, PLOT=False):
+    """
+    Short-term FFT mag for spectogram estimation:
+    Returns:
+        a numpy array (nFFT x numOfShortTermWindows)
+    ARGUMENTS:
+        signal:      the input signal samples
+        Fs:          the sampling freq (in Hz)
+        Win:         the short-term window size (in samples)
+        Step:        the short-term window step (in samples)
+        PLOT:        flag, 1 if results are to be ploted
+    RETURNS:
+    """
+    Win = int(Win)
+    Step = int(Step)
+    signal = numpy.double(signal)
+    signal = signal / (2.0 ** 15)
+    DC = signal.mean()
+    MAX = (numpy.abs(signal)).max()
+    signal = (signal - DC) / (MAX - DC)
+
+    N = len(signal)        # total number of signals
+    curPos = 0
+    countFrames = 0
+    nfft = int(Win / 2)
+    nChroma, nFreqsPerChroma = stChromaFeaturesInit(nfft, Fs)
+    chromaGram = numpy.array([], dtype=numpy.float64)
+
+    while (curPos + Win - 1 < N):
+        countFrames += 1
+        x = signal[curPos:curPos + Win]
+        curPos = curPos + Step
+        X = abs(fft(x))
+        X = X[0:nfft]
+        X = X / len(X)
+        chromaNames, C = stChromaFeatures(X, Fs, nChroma, nFreqsPerChroma)
+        C = C[:, 0]
+        if countFrames == 1:
+            chromaGram = C.T
+        else:
+            chromaGram = numpy.vstack((chromaGram, C.T))
+    FreqAxis = chromaNames
+    TimeAxis = [(t * Step) / Fs for t in range(chromaGram.shape[0])]
+
+    if (PLOT):
+        fig, ax = plt.subplots()
+        chromaGramToPlot = chromaGram.transpose()[::-1, :]
+        Ratio = chromaGramToPlot.shape[1] / (3*chromaGramToPlot.shape[0])        
+        if Ratio < 1:
+            Ratio = 1
+        chromaGramToPlot = numpy.repeat(chromaGramToPlot, Ratio, axis=0)
+        imgplot = plt.imshow(chromaGramToPlot)
+        Fstep = int(nfft / 5.0)
+#        FreqTicks = range(0, int(nfft) + Fstep, Fstep)
+#        FreqTicksLabels = [str(Fs/2-int((f*Fs) / (2*nfft))) for f in FreqTicks]
+        ax.set_yticks(range(Ratio / 2, len(FreqAxis) * Ratio, Ratio))
+        ax.set_yticklabels(FreqAxis[::-1])
+        TStep = countFrames / 3
+        TimeTicks = range(0, countFrames, TStep)
+        TimeTicksLabels = ['%.2f' % (float(t * Step) / Fs) for t in TimeTicks]
+        ax.set_xticks(TimeTicks)
+        ax.set_xticklabels(TimeTicksLabels)
+        ax.set_xlabel('time (secs)')
+        imgplot.set_cmap('jet')
+        plt.colorbar()
+        plt.show()
+
+    return (chromaGram, TimeAxis, FreqAxis)
+
+
+def phormants(x, Fs):
+    N = len(x)
+    w = numpy.hamming(N)
+
+    # Apply window and high pass filter.
+    x1 = x * w   
+    x1 = lfilter([1], [1., 0.63], x1)
+    
+    # Get LPC.    
+    ncoeff = 2 + Fs / 1000
+    A, e, k = lpc(x1, ncoeff)    
+    #A, e, k = lpc(x1, 8)
+
+    # Get roots.
+    rts = numpy.roots(A)
+    rts = [r for r in rts if numpy.imag(r) >= 0]
+
+    # Get angles.
+    angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))
+
+    # Get frequencies.    
+    frqs = sorted(angz * (Fs / (2 * math.pi)))
+
+    return frqs
+def beatExtraction(stFeatures, winSize, PLOT=False):
+    """
+    This function extracts an estimate of the beat rate for a musical signal.
+    ARGUMENTS:
+     - stFeatures:     a numpy array (numOfFeatures x numOfShortTermWindows)
+     - winSize:        window size in seconds
+    RETURNS:
+     - BPM:            estimates of beats per minute
+     - Ratio:          a confidence measure
+    """
+
+    # Features that are related to the beat tracking task:
+    toWatch = [0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]
+
+    maxBeatTime = int(round(2.0 / winSize))
+    HistAll = numpy.zeros((maxBeatTime,))
+    for ii, i in enumerate(toWatch):                                        # for each feature
+        DifThres = 2.0 * (numpy.abs(stFeatures[i, 0:-1] - stFeatures[i, 1::])).mean()    # dif threshold (3 x Mean of Difs)
+        if DifThres<=0:
+            DifThres = 0.0000000000000001        
+        [pos1, _] = utilities.peakdet(stFeatures[i, :], DifThres)           # detect local maxima
+        posDifs = []                                                        # compute histograms of local maxima changes
+        for j in range(len(pos1)-1):
+            posDifs.append(pos1[j+1]-pos1[j])
+        [HistTimes, HistEdges] = numpy.histogram(posDifs, numpy.arange(0.5, maxBeatTime + 1.5))
+        HistCenters = (HistEdges[0:-1] + HistEdges[1::]) / 2.0
+        HistTimes = HistTimes.astype(float) / stFeatures.shape[1]
+        HistAll += HistTimes
+        if PLOT:
+            plt.subplot(9, 2, ii + 1)
+            plt.plot(stFeatures[i, :], 'k')
+            for k in pos1:
+                plt.plot(k, stFeatures[i, k], 'k*')
+            f1 = plt.gca()
+            f1.axes.get_xaxis().set_ticks([])
+            f1.axes.get_yaxis().set_ticks([])
+
+    if PLOT:
+        plt.show(block=False)
+        plt.figure()
+
+    # Get beat as the argmax of the agregated histogram:
+    I = numpy.argmax(HistAll)
+    BPMs = 60 / (HistCenters * winSize)
+    BPM = BPMs[I]
+    # ... and the beat ratio:
+    Ratio = HistAll[I] / HistAll.sum()
+
+    if PLOT:
+        # filter out >500 beats from plotting:
+        HistAll = HistAll[BPMs < 500]
+        BPMs = BPMs[BPMs < 500]
+
+        plt.plot(BPMs, HistAll, 'k')
+        plt.xlabel('Beats per minute')
+        plt.ylabel('Freq Count')
+        plt.show(block=True)
+
+    return BPM, Ratio
+
+
+def stSpectogram(signal, Fs, Win, Step, PLOT=False):
+    """
+    Short-term FFT mag for spectogram estimation:
+    Returns:
+        a numpy array (nFFT x numOfShortTermWindows)
+    ARGUMENTS:
+        signal:      the input signal samples
+        Fs:          the sampling freq (in Hz)
+        Win:         the short-term window size (in samples)
+        Step:        the short-term window step (in samples)
+        PLOT:        flag, 1 if results are to be ploted
+    RETURNS:
+    """
+    Win = int(Win)
+    Step = int(Step)
+    signal = numpy.double(signal)
+    signal = signal / (2.0 ** 15)
+    DC = signal.mean()
+    MAX = (numpy.abs(signal)).max()
+    signal = (signal - DC) / (MAX - DC)
+
+    N = len(signal)        # total number of signals
+    curPos = 0
+    countFrames = 0
+    nfft = int(Win / 2)
+    specgram = numpy.array([], dtype=numpy.float64)
+
+    while (curPos + Win - 1 < N):
+        countFrames += 1
+        x = signal[curPos:curPos+Win]
+        curPos = curPos + Step
+        X = abs(fft(x))
+        X = X[0:nfft]
+        X = X / len(X)
+
+        if countFrames == 1:
+            specgram = X ** 2
+        else:
+            specgram = numpy.vstack((specgram, X))
+
+    FreqAxis = [((f + 1) * Fs) / (2 * nfft) for f in range(specgram.shape[1])]
+    TimeAxis = [(t * Step) / Fs for t in range(specgram.shape[0])]
+
+    if (PLOT):
+        fig, ax = plt.subplots()
+        imgplot = plt.imshow(specgram.transpose()[::-1, :])
+        Fstep = int(nfft / 5.0)
+        FreqTicks = range(0, int(nfft) + Fstep, Fstep)
+        FreqTicksLabels = [str(Fs / 2 - int((f * Fs) / (2 * nfft))) for f in FreqTicks]
+        ax.set_yticks(FreqTicks)
+        ax.set_yticklabels(FreqTicksLabels)
+        TStep = countFrames/3
+        TimeTicks = range(0, countFrames, TStep)
+        TimeTicksLabels = ['%.2f' % (float(t * Step) / Fs) for t in TimeTicks]
+        ax.set_xticks(TimeTicks)
+        ax.set_xticklabels(TimeTicksLabels)
+        ax.set_xlabel('time (secs)')
+        ax.set_ylabel('freq (Hz)')
+        imgplot.set_cmap('jet')
+        plt.colorbar()
+        plt.show()
+
+    return (specgram, TimeAxis, FreqAxis)
+
+
+""" Windowing and feature extraction """
+
+
+def stFeatureExtraction(signal, Fs, Win, Step):
+    """
+    This function implements the shor-term windowing process. For each short-term window a set of features is extracted.
+    This results to a sequence of feature vectors, stored in a numpy matrix.
+
+    ARGUMENTS
+        signal:       the input signal samples
+        Fs:           the sampling freq (in Hz)
+        Win:          the short-term window size (in samples)
+        Step:         the short-term window step (in samples)
+    RETURNS
+        stFeatures:   a numpy array (numOfFeatures x numOfShortTermWindows)
+    """
+
+    Win = int(Win)
+    Step = int(Step)
+
+    # Signal normalization
+    signal = numpy.double(signal)
+
+    signal = signal / (2.0 ** 15)
+    DC = signal.mean()
+    MAX = (numpy.abs(signal)).max()
+    signal = (signal - DC) / (MAX + 0.0000000001)
+
+    N = len(signal)                                # total number of samples
+    curPos = 0
+    countFrames = 0
+    nFFT = Win / 2
+
+    [fbank, freqs] = mfccInitFilterBanks(Fs, nFFT)                # compute the triangular filter banks used in the mfcc calculation
+    nChroma, nFreqsPerChroma = stChromaFeaturesInit(nFFT, Fs)
+
+    numOfTimeSpectralFeatures = 8
+    numOfHarmonicFeatures = 0
+    nceps = 13
+    numOfChromaFeatures = 13
+    totalNumOfFeatures = numOfTimeSpectralFeatures + nceps + numOfHarmonicFeatures + numOfChromaFeatures
+#    totalNumOfFeatures = numOfTimeSpectralFeatures + nceps + numOfHarmonicFeatures
+
+    stFeatures = []
+    while (curPos + Win - 1 < N):                        # for each short-term window until the end of signal
+        countFrames += 1
+        x = signal[curPos:curPos+Win]                    # get current window
+        curPos = curPos + Step                           # update window position
+        X = abs(fft(x))                                  # get fft magnitude
+        X = X[0:nFFT]                                    # normalize fft
+        X = X / len(X)
+        if countFrames == 1:
+            Xprev = X.copy()                             # keep previous fft mag (used in spectral flux)
+        curFV = numpy.zeros((totalNumOfFeatures, 1))
+        curFV[0] = stZCR(x)                              # zero crossing rate
+        curFV[1] = stEnergy(x)                           # short-term energy
+        curFV[2] = stEnergyEntropy(x)                    # short-term entropy of energy
+        [curFV[3], curFV[4]] = stSpectralCentroidAndSpread(X, Fs)    # spectral centroid and spread
+        curFV[5] = stSpectralEntropy(X)                  # spectral entropy
+        curFV[6] = stSpectralFlux(X, Xprev)              # spectral flux
+        curFV[7] = stSpectralRollOff(X, 0.90, Fs)        # spectral rolloff
+        curFV[numOfTimeSpectralFeatures:numOfTimeSpectralFeatures+nceps, 0] = stMFCC(X, fbank, nceps).copy()    # MFCCs
+
+        chromaNames, chromaF = stChromaFeatures(X, Fs, nChroma, nFreqsPerChroma)
+        curFV[numOfTimeSpectralFeatures + nceps: numOfTimeSpectralFeatures + nceps + numOfChromaFeatures - 1] = chromaF
+        curFV[numOfTimeSpectralFeatures + nceps + numOfChromaFeatures - 1] = chromaF.std()
+        stFeatures.append(curFV)
+        # delta features
+        '''
+        if countFrames>1:
+            delta = curFV - prevFV
+            curFVFinal = numpy.concatenate((curFV, delta))            
+        else:
+            curFVFinal = numpy.concatenate((curFV, curFV))
+        prevFV = curFV
+        stFeatures.append(curFVFinal)        
+        '''
+        # end of delta
+        Xprev = X.copy()
+
+    stFeatures = numpy.concatenate(stFeatures, 1)
+    return stFeatures
+
+
+def mtFeatureExtraction(signal, Fs, mtWin, mtStep, stWin, stStep):
+    """
+    Mid-term feature extraction
+    """
+
+    mtWinRatio = int(round(mtWin / stStep))
+    mtStepRatio = int(round(mtStep / stStep))
+
+    mtFeatures = []
+
+    stFeatures = stFeatureExtraction(signal, Fs, stWin, stStep)
+    numOfFeatures = len(stFeatures)
+    numOfStatistics = 2
+
+    mtFeatures = []
+    #for i in range(numOfStatistics * numOfFeatures + 1):
+    for i in range(numOfStatistics * numOfFeatures):
+        mtFeatures.append([])
+
+    for i in range(numOfFeatures):        # for each of the short-term features:
+        curPos = 0
+        N = len(stFeatures[i])
+        while (curPos < N):
+            N1 = curPos
+            N2 = curPos + mtWinRatio
+            if N2 > N:
+                N2 = N
+            curStFeatures = stFeatures[i][N1:N2]
+
+            mtFeatures[i].append(numpy.mean(curStFeatures))
+            mtFeatures[i+numOfFeatures].append(numpy.std(curStFeatures))
+            #mtFeatures[i+2*numOfFeatures].append(numpy.std(curStFeatures) / (numpy.mean(curStFeatures)+0.00000010))
+            curPos += mtStepRatio
+
+    return numpy.array(mtFeatures), stFeatures
+
+
+# TODO
+def stFeatureSpeed(signal, Fs, Win, Step):
+
+    signal = numpy.double(signal)
+    signal = signal / (2.0 ** 15)
+    DC = signal.mean()
+    MAX = (numpy.abs(signal)).max()
+    signal = (signal - DC) / MAX
+    # print (numpy.abs(signal)).max()
+
+    N = len(signal)        # total number of signals
+    curPos = 0
+    countFrames = 0
+
+    lowfreq = 133.33
+    linsc = 200/3.
+    logsc = 1.0711703
+    nlinfil = 13
+    nlogfil = 27
+    nceps = 13
+    nfil = nlinfil + nlogfil
+    nfft = Win / 2
+    if Fs < 8000:
+        nlogfil = 5
+        nfil = nlinfil + nlogfil
+        nfft = Win / 2
+
+    # compute filter banks for mfcc:
+    [fbank, freqs] = mfccInitFilterBanks(Fs, nfft, lowfreq, linsc, logsc, nlinfil, nlogfil)
+
+    numOfTimeSpectralFeatures = 8
+    numOfHarmonicFeatures = 1
+    totalNumOfFeatures = numOfTimeSpectralFeatures + nceps + numOfHarmonicFeatures
+    #stFeatures = numpy.array([], dtype=numpy.float64)
+    stFeatures = []
+
+    while (curPos + Win - 1 < N):
+        countFrames += 1
+        x = signal[curPos:curPos + Win]
+        curPos = curPos + Step
+        X = abs(fft(x))
+        X = X[0:nfft]
+        X = X / len(X)
+        Ex = 0.0
+        El = 0.0
+        X[0:4] = 0
+#        M = numpy.round(0.016 * fs) - 1
+#        R = numpy.correlate(frame, frame, mode='full')
+        stFeatures.append(stHarmonic(x, Fs))
+#        for i in range(len(X)):
+            #if (i < (len(X) / 8)) and (i > (len(X)/40)):
+            #    Ex += X[i]*X[i]
+            #El += X[i]*X[i]
+#        stFeatures.append(Ex / El)
+#        stFeatures.append(numpy.argmax(X))
+#        if curFV[numOfTimeSpectralFeatures+nceps+1]>0:
+#            print curFV[numOfTimeSpectralFeatures+nceps], curFV[numOfTimeSpectralFeatures+nceps+1]
+    return numpy.array(stFeatures)
+
+
+""" Feature Extraction Wrappers
+
+ - The first two feature extraction wrappers are used to extract long-term averaged
+   audio features for a list of WAV files stored in a given category.
+   It is important to note that, one single feature is extracted per WAV file (not the whole sequence of feature vectors)
+
+ """
+
+
+def dirWavFeatureExtraction(dirName, mtWin, mtStep, stWin, stStep, computeBEAT=False):
+    """
+    This function extracts the mid-term features of the WAVE files of a particular folder.
+
+    The resulting feature vector is extracted by long-term averaging the mid-term features.
+    Therefore ONE FEATURE VECTOR is extracted for each WAV file.
+
+    ARGUMENTS:
+        - dirName:        the path of the WAVE directory
+        - mtWin, mtStep:    mid-term window and step (in seconds)
+        - stWin, stStep:    short-term window and step (in seconds)
+    """
+
+    allMtFeatures = numpy.array([])
+    processingTimes = []
+
+    types = ('*.wav', '*.aif',  '*.aiff', '*.mp3','*.au')
+    wavFilesList = []
+    for files in types:
+        wavFilesList.extend(glob.glob(os.path.join(dirName, files)))
+
+    wavFilesList = sorted(wavFilesList)    
+    wavFilesList2 = []
+    for i, wavFile in enumerate(wavFilesList):        
+        print "Analyzing file {0:d} of {1:d}: {2:s}".format(i+1, len(wavFilesList), wavFile.encode('utf-8'))
+        if os.stat(wavFile).st_size == 0:
+            print "   (EMPTY FILE -- SKIPPING)"
+            continue        
+        [Fs, x] = audioBasicIO.readAudioFile(wavFile)            # read file    
+        if isinstance(x, int):
+            continue        
+
+        t1 = time.clock()        
+        x = audioBasicIO.stereo2mono(x)                          # convert stereo to mono                
+        if x.shape[0]<float(Fs)/10:
+            print "  (AUDIO FILE TOO SMALL - SKIPPING)"
+            continue
+        wavFilesList2.append(wavFile)
+        if computeBEAT:                                          # mid-term feature extraction for current file
+            [MidTermFeatures, stFeatures] = mtFeatureExtraction(x, Fs, round(mtWin * Fs), round(mtStep * Fs), round(Fs * stWin), round(Fs * stStep))
+            [beat, beatConf] = beatExtraction(stFeatures, stStep)
+        else:
+            [MidTermFeatures, _] = mtFeatureExtraction(x, Fs, round(mtWin * Fs), round(mtStep * Fs), round(Fs * stWin), round(Fs * stStep))
+
+        MidTermFeatures = numpy.transpose(MidTermFeatures)
+        MidTermFeatures = MidTermFeatures.mean(axis=0)         # long term averaging of mid-term statistics
+        if (not numpy.isnan(MidTermFeatures).any()) and (not numpy.isinf(MidTermFeatures).any()):            
+            if computeBEAT:
+                MidTermFeatures = numpy.append(MidTermFeatures, beat)
+                MidTermFeatures = numpy.append(MidTermFeatures, beatConf)
+            if len(allMtFeatures) == 0:                              # append feature vector
+                allMtFeatures = MidTermFeatures
+            else:
+                allMtFeatures = numpy.vstack((allMtFeatures, MidTermFeatures))
+            t2 = time.clock()
+            duration = float(len(x)) / Fs
+            processingTimes.append((t2 - t1) / duration)
+    if len(processingTimes) > 0:
+        print "Feature extraction complexity ratio: {0:.1f} x realtime".format((1.0 / numpy.mean(numpy.array(processingTimes))))
+    return (allMtFeatures, wavFilesList2)
+
+
+def dirsWavFeatureExtraction(dirNames, mtWin, mtStep, stWin, stStep, computeBEAT=False):
+    '''
+    Same as dirWavFeatureExtraction, but instead of a single dir it takes a list of paths as input and returns a list of feature matrices.
+    EXAMPLE:
+    [features, classNames] =
+           a.dirsWavFeatureExtraction(['audioData/classSegmentsRec/noise','audioData/classSegmentsRec/speech',
+                                       'audioData/classSegmentsRec/brush-teeth','audioData/classSegmentsRec/shower'], 1, 1, 0.02, 0.02);
+
+    It can be used during the training process of a classification model ,
+    in order to get feature matrices from various audio classes (each stored in a seperate path)
+    '''
+
+    # feature extraction for each class:
+    features = []
+    classNames = []
+    fileNames = []
+    for i, d in enumerate(dirNames):
+        [f, fn] = dirWavFeatureExtraction(d, mtWin, mtStep, stWin, stStep, computeBEAT=computeBEAT)
+        if f.shape[0] > 0:       # if at least one audio file has been found in the provided folder:
+            features.append(f)
+            fileNames.append(fn)
+            if d[-1] == "/":
+                classNames.append(d.split(os.sep)[-2])
+            else:
+                classNames.append(d.split(os.sep)[-1])
+    return features, classNames, fileNames
+
+
+def dirWavFeatureExtractionNoAveraging(dirName, mtWin, mtStep, stWin, stStep):
+    """
+    This function extracts the mid-term features of the WAVE files of a particular folder without averaging each file.
+
+    ARGUMENTS:
+        - dirName:          the path of the WAVE directory
+        - mtWin, mtStep:    mid-term window and step (in seconds)
+        - stWin, stStep:    short-term window and step (in seconds)
+    RETURNS:
+        - X:                A feature matrix
+        - Y:                A matrix of file labels
+        - filenames:
+    """
+
+    allMtFeatures = numpy.array([])
+    signalIndices = numpy.array([])
+    processingTimes = []
+
+    types = ('*.wav', '*.aif',  '*.aiff')
+    wavFilesList = []
+    for files in types:
+        wavFilesList.extend(glob.glob(os.path.join(dirName, files)))
+
+    wavFilesList = sorted(wavFilesList)
+
+    for i, wavFile in enumerate(wavFilesList):
+        [Fs, x] = audioBasicIO.readAudioFile(wavFile)            # read file
+        if isinstance(x, int):
+            continue        
+        
+        x = audioBasicIO.stereo2mono(x)                          # convert stereo to mono
+        [MidTermFeatures, _] = mtFeatureExtraction(x, Fs, round(mtWin * Fs), round(mtStep * Fs), round(Fs * stWin), round(Fs * stStep))  # mid-term feature
+
+        MidTermFeatures = numpy.transpose(MidTermFeatures)
+#        MidTermFeatures = MidTermFeatures.mean(axis=0)        # long term averaging of mid-term statistics
+        if len(allMtFeatures) == 0:                # append feature vector
+            allMtFeatures = MidTermFeatures
+            signalIndices = numpy.zeros((MidTermFeatures.shape[0], ))
+        else:
+            allMtFeatures = numpy.vstack((allMtFeatures, MidTermFeatures))
+            signalIndices = numpy.append(signalIndices, i * numpy.ones((MidTermFeatures.shape[0], )))
+
+    return (allMtFeatures, signalIndices, wavFilesList)
+
+
+# The following two feature extraction wrappers extract features for given audio files, however
+# NO LONG-TERM AVERAGING is performed. Therefore, the output for each audio file is NOT A SINGLE FEATURE VECTOR
+# but a whole feature matrix.
+#
+# Also, another difference between the following two wrappers and the previous is that they NO LONG-TERM AVERAGING IS PERFORMED.
+# In other words, the WAV files in these functions are not used as uniform samples that need to be averaged but as sequences
+
+def mtFeatureExtractionToFile(fileName, midTermSize, midTermStep, shortTermSize, shortTermStep, outPutFile,
+                              storeStFeatures=False, storeToCSV=False, PLOT=False):
+    """
+    This function is used as a wrapper to:
+    a) read the content of a WAV file
+    b) perform mid-term feature extraction on that signal
+    c) write the mid-term feature sequences to a numpy file
+    """
+    [Fs, x] = audioBasicIO.readAudioFile(fileName)            # read the wav file
+    x = audioBasicIO.stereo2mono(x)                           # convert to MONO if required
+    if storeStFeatures:
+        [mtF, stF] = mtFeatureExtraction(x, Fs, round(Fs * midTermSize), round(Fs * midTermStep), round(Fs * shortTermSize), round(Fs * shortTermStep))
+    else:
+        [mtF, _] = mtFeatureExtraction(x, Fs, round(Fs*midTermSize), round(Fs * midTermStep), round(Fs * shortTermSize), round(Fs * shortTermStep))
+
+    numpy.save(outPutFile, mtF)                              # save mt features to numpy file
+    if PLOT:
+        print "Mid-term numpy file: " + outPutFile + ".npy saved"
+    if storeToCSV:
+        numpy.savetxt(outPutFile+".csv", mtF.T, delimiter=",")
+        if PLOT:
+            print "Mid-term CSV file: " + outPutFile + ".csv saved"
+
+    if storeStFeatures:
+        numpy.save(outPutFile+"_st", stF)                    # save st features to numpy file
+        if PLOT:
+            print "Short-term numpy file: " + outPutFile + "_st.npy saved"
+        if storeToCSV:
+            numpy.savetxt(outPutFile+"_st.csv", stF.T, delimiter=",")    # store st features to CSV file
+            if PLOT:
+                print "Short-term CSV file: " + outPutFile + "_st.csv saved"
+
+
+def mtFeatureExtractionToFileDir(dirName, midTermSize, midTermStep, shortTermSize, shortTermStep, storeStFeatures=False, storeToCSV=False, PLOT=False):
+    types = (dirName + os.sep + '*.wav', )
+    filesToProcess = []
+    for files in types:
+        filesToProcess.extend(glob.glob(files))
+    for f in filesToProcess:
+        outPath = f
+        mtFeatureExtractionToFile(f, midTermSize, midTermStep, shortTermSize, shortTermStep, outPath, storeStFeatures, storeToCSV, PLOT)
diff --git a/audioSegmentation.py b/audioSegmentation.py
new file mode 100644
index 0000000..d2d58eb
--- /dev/null
+++ b/audioSegmentation.py
@@ -0,0 +1,1030 @@
+import numpy
+import sklearn.cluster
+import time
+import scipy
+import os
+import audioFeatureExtraction as aF
+import audioTrainTest as aT
+import audioBasicIO
+import matplotlib.pyplot as plt
+from scipy.spatial import distance
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import sklearn.discriminant_analysis
+import csv
+import os.path
+import sklearn
+import sklearn.cluster
+import hmmlearn.hmm
+import cPickle
+import glob
+
+""" General utility functions """
+
+
+def smoothMovingAvg(inputSignal, windowLen=11):
+    windowLen = int(windowLen)
+    if inputSignal.ndim != 1:
+        raise ValueError("")
+    if inputSignal.size < windowLen:
+        raise ValueError("Input vector needs to be bigger than window size.")
+    if windowLen < 3:
+        return inputSignal
+    s = numpy.r_[2*inputSignal[0] - inputSignal[windowLen-1::-1], inputSignal, 2*inputSignal[-1]-inputSignal[-1:-windowLen:-1]]
+    w = numpy.ones(windowLen, 'd')
+    y = numpy.convolve(w/w.sum(), s, mode='same')
+    return y[windowLen:-windowLen+1]
+
+
+def selfSimilarityMatrix(featureVectors):
+    '''
+    This function computes the self-similarity matrix for a sequence of feature vectors.
+    ARGUMENTS:
+     - featureVectors:     a numpy matrix (nDims x nVectors) whose i-th column corresponds to the i-th feature vector
+
+    RETURNS:
+     - S:             the self-similarity matrix (nVectors x nVectors)
+    '''
+
+    [nDims, nVectors] = featureVectors.shape
+    [featureVectors2, MEAN, STD] = aT.normalizeFeatures([featureVectors.T])
+    featureVectors2 = featureVectors2[0].T
+    S = 1.0 - distance.squareform(distance.pdist(featureVectors2.T, 'cosine'))
+    return S
+
+
+def flags2segs(Flags, window):
+    '''
+    ARGUMENTS:
+     - Flags:     a sequence of class flags (per time window)
+     - window:    window duration (in seconds)
+
+    RETURNS:
+     - segs:    a sequence of segment's limits: segs[i,0] is start and segs[i,1] are start and end point of segment i
+     - classes:    a sequence of class flags: class[i] is the class ID of the i-th segment
+    '''
+
+    preFlag = 0
+    curFlag = 0
+    numOfSegments = 0
+
+    curVal = Flags[curFlag]
+    segsList = []
+    classes = []
+    while (curFlag < len(Flags) - 1):
+        stop = 0
+        preFlag = curFlag
+        preVal = curVal
+        while (stop == 0):
+            curFlag = curFlag + 1
+            tempVal = Flags[curFlag]
+            if ((tempVal != curVal) | (curFlag == len(Flags) - 1)):  # stop
+                numOfSegments = numOfSegments + 1
+                stop = 1
+                curSegment = curVal
+                curVal = Flags[curFlag]
+                segsList.append((curFlag * window))
+                classes.append(preVal)
+    segs = numpy.zeros((len(segsList), 2))
+
+    for i in range(len(segsList)):
+        if i > 0:
+            segs[i, 0] = segsList[i-1]
+        segs[i, 1] = segsList[i]
+    return (segs, classes)
+
+
+def segs2flags(segStart, segEnd, segLabel, winSize):
+    '''
+    This function converts segment endpoints and respective segment labels to fix-sized class labels.
+    ARGUMENTS:
+     - segStart:    segment start points (in seconds)
+     - segEnd:    segment endpoints (in seconds)
+     - segLabel:    segment labels
+      - winSize:    fix-sized window (in seconds)
+    RETURNS:
+     - flags:    numpy array of class indices
+     - classNames:    list of classnames (strings)
+    '''
+    flags = []
+    classNames = list(set(segLabel))
+    curPos = winSize / 2.0
+    while curPos < segEnd[-1]:
+        for i in range(len(segStart)):
+            if curPos > segStart[i] and curPos <= segEnd[i]:
+                break
+        flags.append(classNames.index(segLabel[i]))
+        curPos += winSize
+    return numpy.array(flags), classNames
+
+def computePreRec(CM, classNames):
+    '''
+    This function computes the Precision, Recall and F1 measures, given a confusion matrix
+    '''
+    numOfClasses = CM.shape[0]
+    if len(classNames) != numOfClasses:
+        print "Error in computePreRec! Confusion matrix and classNames list must be of the same size!"
+        return
+    Precision = []
+    Recall = []
+    F1 = []    
+    for i, c in enumerate(classNames):
+        Precision.append(CM[i,i] / numpy.sum(CM[:,i]))
+        Recall.append(CM[i,i] / numpy.sum(CM[i,:]))
+        F1.append( 2 * Precision[-1] * Recall[-1] / (Precision[-1] + Recall[-1]))
+    return Recall, Precision, F1
+
+
+def readSegmentGT(gtFile):
+    '''
+    This function reads a segmentation ground truth file, following a simple CSV format with the following columns:
+    <segment start>,<segment end>,<class label>
+
+    ARGUMENTS:
+     - gtFile:       the path of the CSV segment file
+    RETURNS:
+     - segStart:     a numpy array of segments' start positions
+     - segEnd:       a numpy array of segments' ending positions
+     - segLabel:     a list of respective class labels (strings)
+    '''
+    f = open(gtFile, "rb")
+    reader = csv.reader(f, delimiter=',')
+    segStart = []
+    segEnd = []
+    segLabel = []
+    for row in reader:
+        if len(row) == 3:
+            segStart.append(float(row[0]))
+            segEnd.append(float(row[1]))
+            #if row[2]!="other":
+            #    segLabel.append((row[2]))
+            #else:
+            #    segLabel.append("silence")
+            segLabel.append((row[2]))
+    return numpy.array(segStart), numpy.array(segEnd), segLabel
+
+
+def plotSegmentationResults(flagsInd, flagsIndGT, classNames, mtStep, ONLY_EVALUATE=False):
+    '''
+    This function plots statistics on the classification-segmentation results produced either by the fix-sized supervised method or the HMM method.
+    It also computes the overall accuracy achieved by the respective method if ground-truth is available.
+    '''    
+    flags = [classNames[int(f)] for f in flagsInd]
+    (segs, classes) = flags2segs(flags, mtStep)    
+    minLength = min(flagsInd.shape[0], flagsIndGT.shape[0])    
+    if minLength > 0:
+        accuracy = numpy.sum(flagsInd[0:minLength] == flagsIndGT[0:minLength]) / float(minLength)
+    else:
+        accuracy = -1
+
+    if not ONLY_EVALUATE:
+        Duration = segs[-1, 1]
+        SPercentages = numpy.zeros((len(classNames), 1))
+        Percentages = numpy.zeros((len(classNames), 1))
+        AvDurations = numpy.zeros((len(classNames), 1))
+
+        for iSeg in range(segs.shape[0]):
+            SPercentages[classNames.index(classes[iSeg])] += (segs[iSeg, 1]-segs[iSeg, 0])
+
+        for i in range(SPercentages.shape[0]):
+            Percentages[i] = 100.0 * SPercentages[i] / Duration
+            S = sum(1 for c in classes if c == classNames[i])
+            if S > 0:
+                AvDurations[i] = SPercentages[i] / S
+            else:
+                AvDurations[i] = 0.0
+
+        for i in range(Percentages.shape[0]):
+            print classNames[i], Percentages[i], AvDurations[i]
+
+        font = {'size': 10}
+        plt.rc('font', **font)
+
+        fig = plt.figure()
+        ax1 = fig.add_subplot(211)
+        ax1.set_yticks(numpy.array(range(len(classNames))))
+        ax1.axis((0, Duration, -1, len(classNames)))
+        ax1.set_yticklabels(classNames)
+        ax1.plot(numpy.array(range(len(flagsInd))) * mtStep + mtStep / 2.0, flagsInd)
+        if flagsIndGT.shape[0] > 0:
+            ax1.plot(numpy.array(range(len(flagsIndGT))) * mtStep + mtStep / 2.0, flagsIndGT + 0.05, '--r')
+        plt.xlabel("time (seconds)")
+        if accuracy >= 0:
+            plt.title('Accuracy = {0:.1f}%'.format(100.0 * accuracy))
+
+        ax2 = fig.add_subplot(223)
+        plt.title("Classes percentage durations")
+        ax2.axis((0, len(classNames) + 1, 0, 100))
+        ax2.set_xticks(numpy.array(range(len(classNames) + 1)))
+        ax2.set_xticklabels([" "] + classNames)
+        ax2.bar(numpy.array(range(len(classNames))) + 0.5, Percentages)
+
+        ax3 = fig.add_subplot(224)
+        plt.title("Segment average duration per class")
+        ax3.axis((0, len(classNames)+1, 0, AvDurations.max()))
+        ax3.set_xticks(numpy.array(range(len(classNames) + 1)))
+        ax3.set_xticklabels([" "] + classNames)
+        ax3.bar(numpy.array(range(len(classNames))) + 0.5, AvDurations)
+        fig.tight_layout()
+        plt.show()
+    return accuracy
+
+
+def evaluateSpeakerDiarization(flags, flagsGT):
+
+    minLength = min(flags.shape[0], flagsGT.shape[0])
+    flags = flags[0:minLength]
+    flagsGT = flagsGT[0:minLength]
+
+    uFlags = numpy.unique(flags)
+    uFlagsGT = numpy.unique(flagsGT)
+
+    # compute contigency table:
+    cMatrix = numpy.zeros((uFlags.shape[0], uFlagsGT.shape[0]))
+    for i in range(minLength):
+        cMatrix[int(numpy.nonzero(uFlags == flags[i])[0]), int(numpy.nonzero(uFlagsGT == flagsGT[i])[0])] += 1.0
+
+    Nc, Ns = cMatrix.shape
+    N_s = numpy.sum(cMatrix, axis=0)
+    N_c = numpy.sum(cMatrix, axis=1)
+    N = numpy.sum(cMatrix)
+
+    purityCluster = numpy.zeros((Nc, ))
+    puritySpeaker = numpy.zeros((Ns, ))
+    # compute cluster purity:
+    for i in range(Nc):
+        purityCluster[i] = numpy.max((cMatrix[i, :])) / (N_c[i])
+
+    for j in range(Ns):
+        puritySpeaker[j] = numpy.max((cMatrix[:, j])) / (N_s[j])
+
+    purityClusterMean = numpy.sum(purityCluster * N_c) / N
+    puritySpeakerMean = numpy.sum(puritySpeaker * N_s) / N
+
+    return purityClusterMean, puritySpeakerMean
+
+
+def trainHMM_computeStatistics(features, labels):
+    '''
+    This function computes the statistics used to train an HMM joint segmentation-classification model
+    using a sequence of sequential features and respective labels
+
+    ARGUMENTS:
+     - features:    a numpy matrix of feature vectors (numOfDimensions x numOfWindows)
+     - labels:    a numpy array of class indices (numOfWindows x 1)
+    RETURNS:
+     - startprob:    matrix of prior class probabilities (numOfClasses x 1)
+     - transmat:    transition matrix (numOfClasses x numOfClasses)
+     - means:    means matrix (numOfDimensions x 1)
+     - cov:        deviation matrix (numOfDimensions x 1)
+    '''
+    uLabels = numpy.unique(labels)
+    nComps = len(uLabels)
+
+    nFeatures = features.shape[0]
+
+    if features.shape[1] < labels.shape[0]:
+        print "trainHMM warning: number of short-term feature vectors must be greater or equal to the labels length!"
+        labels = labels[0:features.shape[1]]
+
+    # compute prior probabilities:
+    startprob = numpy.zeros((nComps,))
+    for i, u in enumerate(uLabels):
+        startprob[i] = numpy.count_nonzero(labels == u)
+    startprob = startprob / startprob.sum()                # normalize prior probabilities
+
+    # compute transition matrix:
+    transmat = numpy.zeros((nComps, nComps))
+    for i in range(labels.shape[0]-1):
+        transmat[int(labels[i]), int(labels[i + 1])] += 1
+    for i in range(nComps):                     # normalize rows of transition matrix:
+        transmat[i, :] /= transmat[i, :].sum()
+
+    means = numpy.zeros((nComps, nFeatures))
+    for i in range(nComps):
+        means[i, :] = numpy.matrix(features[:, numpy.nonzero(labels == uLabels[i])[0]].mean(axis=1))
+
+    cov = numpy.zeros((nComps, nFeatures))
+    for i in range(nComps):
+        #cov[i,:,:] = numpy.cov(features[:,numpy.nonzero(labels==uLabels[i])[0]])  # use this lines if HMM using full gaussian distributions are to be used!
+        cov[i, :] = numpy.std(features[:, numpy.nonzero(labels == uLabels[i])[0]], axis=1)
+
+    return startprob, transmat, means, cov
+
+
+def trainHMM_fromFile(wavFile, gtFile, hmmModelName, mtWin, mtStep):
+    '''
+    This function trains a HMM model for segmentation-classification using a single annotated audio file
+    ARGUMENTS:
+     - wavFile:        the path of the audio filename
+     - gtFile:         the path of the ground truth filename
+                       (a csv file of the form <segment start in seconds>,<segment end in seconds>,<segment label> in each row
+     - hmmModelName:   the name of the HMM model to be stored
+     - mtWin:          mid-term window size
+     - mtStep:         mid-term window step
+    RETURNS:
+     - hmm:            an object to the resulting HMM
+     - classNames:     a list of classNames
+
+    After training, hmm, classNames, along with the mtWin and mtStep values are stored in the hmmModelName file
+    '''
+
+    [segStart, segEnd, segLabels] = readSegmentGT(gtFile)                        # read ground truth data
+    flags, classNames = segs2flags(segStart, segEnd, segLabels, mtStep)          # convert to fix-sized sequence of flags
+
+    [Fs, x] = audioBasicIO.readAudioFile(wavFile)                                # read audio data
+    #F = aF.stFeatureExtraction(x, Fs, 0.050*Fs, 0.050*Fs);
+    [F, _] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs * 0.050), round(Fs * 0.050))    # feature extraction
+    startprob, transmat, means, cov = trainHMM_computeStatistics(F, flags)                    # compute HMM statistics (priors, transition matrix, etc)
+    
+    hmm = hmmlearn.hmm.GaussianHMM(startprob.shape[0], "diag")            # hmm training
+
+    hmm.startprob_ = startprob
+    hmm.transmat_ = transmat    
+    hmm.means_ = means
+    hmm.covars_ = cov
+    
+    fo = open(hmmModelName, "wb")                                                             # output to file
+    cPickle.dump(hmm, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    cPickle.dump(classNames, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    cPickle.dump(mtWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    cPickle.dump(mtStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    fo.close()
+
+    return hmm, classNames
+
+
+def trainHMM_fromDir(dirPath, hmmModelName, mtWin, mtStep):
+    '''
+    This function trains a HMM model for segmentation-classification using a where WAV files and .segment (ground-truth files) are stored
+    ARGUMENTS:
+     - dirPath:        the path of the data diretory
+     - hmmModelName:    the name of the HMM model to be stored
+     - mtWin:        mid-term window size
+     - mtStep:        mid-term window step
+    RETURNS:
+     - hmm:            an object to the resulting HMM
+     - classNames:        a list of classNames
+
+    After training, hmm, classNames, along with the mtWin and mtStep values are stored in the hmmModelName file
+    '''
+
+    flagsAll = numpy.array([])
+    classesAll = []
+    for i, f in enumerate(glob.glob(dirPath + os.sep + '*.wav')):               # for each WAV file
+        wavFile = f
+        gtFile = f.replace('.wav', '.segments')                                 # open for annotated file
+        if not os.path.isfile(gtFile):                                          # if current WAV file does not have annotation -> skip
+            continue
+        [segStart, segEnd, segLabels] = readSegmentGT(gtFile)                   # read GT data
+        flags, classNames = segs2flags(segStart, segEnd, segLabels, mtStep)     # convert to flags
+        for c in classNames:                                                    # update classnames:
+            if c not in classesAll:
+                classesAll.append(c)
+        [Fs, x] = audioBasicIO.readAudioFile(wavFile)                           # read audio data
+        [F, _] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs * 0.050), round(Fs * 0.050))     # feature extraction
+
+        lenF = F.shape[1]
+        lenL = len(flags)
+        MIN = min(lenF, lenL)
+        F = F[:, 0:MIN]
+        flags = flags[0:MIN]
+
+        flagsNew = []
+        for j, fl in enumerate(flags):      # append features and labels
+            flagsNew.append(classesAll.index(classNames[flags[j]]))
+
+        flagsAll = numpy.append(flagsAll, numpy.array(flagsNew))
+
+        if i == 0:
+            Fall = F
+        else:
+            Fall = numpy.concatenate((Fall, F), axis=1)
+    startprob, transmat, means, cov = trainHMM_computeStatistics(Fall, flagsAll)        # compute HMM statistics
+    hmm = hmmlearn.hmm.GaussianHMM(startprob.shape[0], "diag")      # train HMM
+    hmm.startprob_ = startprob
+    hmm.transmat_ = transmat        
+    hmm.means_ = means
+    hmm.covars_ = cov
+
+    fo = open(hmmModelName, "wb")   # save HMM model
+    cPickle.dump(hmm, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    cPickle.dump(classesAll, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    cPickle.dump(mtWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    cPickle.dump(mtStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+    fo.close()
+
+    return hmm, classesAll
+
+
+def hmmSegmentation(wavFileName, hmmModelName, PLOT=False, gtFileName=""):
+    [Fs, x] = audioBasicIO.readAudioFile(wavFileName)          # read audio data
+
+    try:
+        fo = open(hmmModelName, "rb")
+    except IOError:
+        print "didn't find file"
+        return
+
+    try:
+        hmm = cPickle.load(fo)
+        classesAll = cPickle.load(fo)
+        mtWin = cPickle.load(fo)
+        mtStep = cPickle.load(fo)
+    except:
+        fo.close()
+    fo.close()
+
+    #Features = audioFeatureExtraction.stFeatureExtraction(x, Fs, 0.050*Fs, 0.050*Fs);    # feature extraction
+    [Features, _] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs * 0.050), round(Fs * 0.050))
+    flagsInd = hmm.predict(Features.T)                            # apply model
+    #for i in range(len(flagsInd)):
+    #    if classesAll[flagsInd[i]]=="silence":
+    #        flagsInd[i]=classesAll.index("speech")
+                   
+                                                                             # plot results
+    if os.path.isfile(gtFileName):
+        [segStart, segEnd, segLabels] = readSegmentGT(gtFileName)
+        flagsGT, classNamesGT = segs2flags(segStart, segEnd, segLabels, mtStep)
+        flagsGTNew = []
+        for j, fl in enumerate(flagsGT):                        # "align" labels with GT
+            if classNamesGT[flagsGT[j]] in classesAll:
+                flagsGTNew.append(classesAll.index(classNamesGT[flagsGT[j]]))
+            else:
+                flagsGTNew.append(-1)
+        CM = numpy.zeros((len(classNamesGT), len(classNamesGT)))
+        flagsIndGT = numpy.array(flagsGTNew)
+        for i in range(min(flagsInd.shape[0], flagsIndGT.shape[0])):
+            CM[int(flagsIndGT[i]),int(flagsInd[i])] += 1                
+    else:
+        flagsIndGT = numpy.array([])    
+    acc = plotSegmentationResults(flagsInd, flagsIndGT, classesAll, mtStep, not PLOT)
+    if acc >= 0:
+        print "Overall Accuracy: {0:.2f}".format(acc)
+        return (flagsInd, classNamesGT, acc, CM)
+    else:
+        return (flagsInd, classesAll, -1, -1)
+
+
+
+def mtFileClassification(inputFile, modelName, modelType, plotResults=False, gtFile=""):
+    '''
+    This function performs mid-term classification of an audio stream.
+    Towards this end, supervised knowledge is used, i.e. a pre-trained classifier.
+    ARGUMENTS:
+        - inputFile:        path of the input WAV file
+        - modelName:        name of the classification model
+        - modelType:        svm or knn depending on the classifier type
+        - plotResults:      True if results are to be plotted using matplotlib along with a set of statistics
+
+    RETURNS:
+          - segs:           a sequence of segment's endpoints: segs[i] is the endpoint of the i-th segment (in seconds)
+          - classes:        a sequence of class flags: class[i] is the class ID of the i-th segment
+    '''
+
+    if not os.path.isfile(modelName):
+        print "mtFileClassificationError: input modelType not found!"
+        return (-1, -1, -1, -1)
+    # Load classifier:
+    if (modelType == 'svm') or (modelType == 'svm_rbf'):
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadSVModel(modelName)
+    elif modelType == 'knn':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadKNNModel(modelName)
+    elif modelType == 'randomforest':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadRandomForestModel(modelName)
+    elif modelType == 'gradientboosting':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadGradientBoostingModel(modelName)
+    elif modelType == 'extratrees':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = aT.loadExtraTreesModel(modelName)
+
+
+    if computeBEAT:
+        print "Model " + modelName + " contains long-term music features (beat etc) and cannot be used in segmentation"
+        return (-1, -1, -1, -1)
+    [Fs, x] = audioBasicIO.readAudioFile(inputFile)        # load input file
+    if Fs == -1:                                           # could not read file
+        return (-1, -1, -1, -1)
+    x = audioBasicIO.stereo2mono(x)                        # convert stereo (if) to mono
+    Duration = len(x) / Fs
+    # mid-term feature extraction:
+    [MidTermFeatures, _] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs * stWin), round(Fs * stStep))
+    flags = []
+    Ps = []
+    flagsInd = []
+    for i in range(MidTermFeatures.shape[1]):              # for each feature vector (i.e. for each fix-sized segment):
+        curFV = (MidTermFeatures[:, i] - MEAN) / STD       # normalize current feature vector
+        [Result, P] = aT.classifierWrapper(Classifier, modelType, curFV)    # classify vector
+        flagsInd.append(Result)
+        flags.append(classNames[int(Result)])              # update class label matrix
+        Ps.append(numpy.max(P))                            # update probability matrix
+    flagsInd = numpy.array(flagsInd)
+
+    # 1-window smoothing
+    for i in range(1, len(flagsInd) - 1):
+        if flagsInd[i-1] == flagsInd[i + 1]:
+            flagsInd[i] = flagsInd[i + 1]
+    (segs, classes) = flags2segs(flags, mtStep)            # convert fix-sized flags to segments and classes
+    segs[-1] = len(x) / float(Fs)
+
+    # Load grount-truth:        
+    if os.path.isfile(gtFile):
+        [segStartGT, segEndGT, segLabelsGT] = readSegmentGT(gtFile)
+        flagsGT, classNamesGT = segs2flags(segStartGT, segEndGT, segLabelsGT, mtStep)
+        flagsIndGT = []
+        for j, fl in enumerate(flagsGT):                    # "align" labels with GT
+            if classNamesGT[flagsGT[j]] in classNames:
+                flagsIndGT.append(classNames.index(classNamesGT[flagsGT[j]]))
+            else:
+                flagsIndGT.append(-1)
+        flagsIndGT = numpy.array(flagsIndGT)        
+        CM = numpy.zeros((len(classNamesGT), len(classNamesGT)))
+        for i in range(min(flagsInd.shape[0], flagsIndGT.shape[0])):
+            CM[int(flagsIndGT[i]),int(flagsInd[i])] += 1        
+    else:
+        CM = []
+        flagsIndGT = numpy.array([])
+    acc = plotSegmentationResults(flagsInd, flagsIndGT, classNames, mtStep, not plotResults)
+    if acc >= 0:
+        print "Overall Accuracy: {0:.3f}".format(acc)    
+        return (flagsInd, classNamesGT, acc, CM)
+    else:
+        return (flagsInd, classNames, acc, CM)
+
+
+def evaluateSegmentationClassificationDir(dirName, modelName, methodName):
+    flagsAll = numpy.array([])
+    classesAll = []
+    accuracys = []
+    
+    for i, f in enumerate(glob.glob(dirName + os.sep + '*.wav')):            # for each WAV file
+        wavFile = f
+        print wavFile
+        gtFile = f.replace('.wav', '.segments')                             # open for annotated file
+
+        if methodName.lower() in ["svm", "svm_rbf", "knn","randomforest","gradientboosting","extratrees"]:
+            flagsInd, classNames, acc, CMt = mtFileClassification(wavFile, modelName, methodName, False, gtFile)
+        else:
+            flagsInd, classNames, acc, CMt = hmmSegmentation(wavFile, modelName, False, gtFile)
+        if acc > -1:
+            if i==0:
+                CM = numpy.copy(CMt)
+            else:                
+                CM = CM + CMt
+            accuracys.append(acc)
+            print CMt, classNames
+            print CM
+            [Rec, Pre, F1] = computePreRec(CMt, classNames)
+
+    CM = CM / numpy.sum(CM)
+    [Rec, Pre, F1] = computePreRec(CM, classNames)
+
+    print " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
+    print "Average Accuracy: {0:.1f}".format(100.0*numpy.array(accuracys).mean())
+    print "Average Recall: {0:.1f}".format(100.0*numpy.array(Rec).mean())
+    print "Average Precision: {0:.1f}".format(100.0*numpy.array(Pre).mean())
+    print "Average F1: {0:.1f}".format(100.0*numpy.array(F1).mean())    
+    print "Median Accuracy: {0:.1f}".format(100.0*numpy.median(numpy.array(accuracys)))
+    print "Min Accuracy: {0:.1f}".format(100.0*numpy.array(accuracys).min())
+    print "Max Accuracy: {0:.1f}".format(100.0*numpy.array(accuracys).max())
+
+
+def silenceRemoval(x, Fs, stWin, stStep, smoothWindow=0.5, Weight=0.5, plot=False):
+    '''
+    Event Detection (silence removal)
+    ARGUMENTS:
+         - x:                the input audio signal
+         - Fs:               sampling freq
+         - stWin, stStep:    window size and step in seconds
+         - smoothWindow:     (optinal) smooth window (in seconds)
+         - Weight:           (optinal) weight factor (0 < Weight < 1) the higher, the more strict
+         - plot:             (optinal) True if results are to be plotted
+    RETURNS:
+         - segmentLimits:    list of segment limits in seconds (e.g [[0.1, 0.9], [1.4, 3.0]] means that
+                    the resulting segments are (0.1 - 0.9) seconds and (1.4, 3.0) seconds
+    '''
+
+    if Weight >= 1:
+        Weight = 0.99
+    if Weight <= 0:
+        Weight = 0.01
+
+    # Step 1: feature extraction
+    x = audioBasicIO.stereo2mono(x)                        # convert to mono
+    ShortTermFeatures = aF.stFeatureExtraction(x, Fs, stWin * Fs, stStep * Fs)        # extract short-term features
+
+    # Step 2: train binary SVM classifier of low vs high energy frames
+    EnergySt = ShortTermFeatures[1, :]                  # keep only the energy short-term sequence (2nd feature)
+    E = numpy.sort(EnergySt)                            # sort the energy feature values:
+    L1 = int(len(E) / 10)                               # number of 10% of the total short-term windows
+    T1 = numpy.mean(E[0:L1]) + 0.000000000000001                 # compute "lower" 10% energy threshold
+    T2 = numpy.mean(E[-L1:-1]) + 0.000000000000001                # compute "higher" 10% energy threshold
+    Class1 = ShortTermFeatures[:, numpy.where(EnergySt <= T1)[0]]         # get all features that correspond to low energy
+    Class2 = ShortTermFeatures[:, numpy.where(EnergySt >= T2)[0]]         # get all features that correspond to high energy
+    featuresSS = [Class1.T, Class2.T]                                    # form the binary classification task and ...
+
+    [featuresNormSS, MEANSS, STDSS] = aT.normalizeFeatures(featuresSS)   # normalize and ...
+    SVM = aT.trainSVM(featuresNormSS, 1.0)                               # train the respective SVM probabilistic model (ONSET vs SILENCE)
+
+    # Step 3: compute onset probability based on the trained SVM
+    ProbOnset = []
+    for i in range(ShortTermFeatures.shape[1]):                    # for each frame
+        curFV = (ShortTermFeatures[:, i] - MEANSS) / STDSS         # normalize feature vector
+        ProbOnset.append(SVM.predict_proba(curFV.reshape(1,-1))[0][1])           # get SVM probability (that it belongs to the ONSET class)
+    ProbOnset = numpy.array(ProbOnset)
+    ProbOnset = smoothMovingAvg(ProbOnset, smoothWindow / stStep)  # smooth probability
+
+    # Step 4A: detect onset frame indices:
+    ProbOnsetSorted = numpy.sort(ProbOnset)                        # find probability Threshold as a weighted average of top 10% and lower 10% of the values
+    Nt = ProbOnsetSorted.shape[0] / 10
+    T = (numpy.mean((1 - Weight) * ProbOnsetSorted[0:Nt]) + Weight * numpy.mean(ProbOnsetSorted[-Nt::]))
+
+    MaxIdx = numpy.where(ProbOnset > T)[0]                         # get the indices of the frames that satisfy the thresholding
+    i = 0
+    timeClusters = []
+    segmentLimits = []
+
+    # Step 4B: group frame indices to onset segments
+    while i < len(MaxIdx):                                         # for each of the detected onset indices
+        curCluster = [MaxIdx[i]]
+        if i == len(MaxIdx)-1:
+            break
+        while MaxIdx[i+1] - curCluster[-1] <= 2:
+            curCluster.append(MaxIdx[i+1])
+            i += 1
+            if i == len(MaxIdx)-1:
+                break
+        i += 1
+        timeClusters.append(curCluster)
+        segmentLimits.append([curCluster[0] * stStep, curCluster[-1] * stStep])
+
+    # Step 5: Post process: remove very small segments:
+    minDuration = 0.2
+    segmentLimits2 = []
+    for s in segmentLimits:
+        if s[1] - s[0] > minDuration:
+            segmentLimits2.append(s)
+    segmentLimits = segmentLimits2
+
+    if plot:
+        timeX = numpy.arange(0, x.shape[0] / float(Fs), 1.0 / Fs)
+
+        plt.subplot(2, 1, 1)
+        plt.plot(timeX, x)
+        for s in segmentLimits:
+            plt.axvline(x=s[0])
+            plt.axvline(x=s[1])
+        plt.subplot(2, 1, 2)
+        plt.plot(numpy.arange(0, ProbOnset.shape[0] * stStep, stStep), ProbOnset)
+        plt.title('Signal')
+        for s in segmentLimits:
+            plt.axvline(x=s[0])
+            plt.axvline(x=s[1])
+        plt.title('SVM Probability')
+        plt.show()
+
+    return segmentLimits
+
+
+def speakerDiarization(fileName, numOfSpeakers, mtSize=2.0, mtStep=0.2, stWin=0.05, LDAdim=35, PLOT=False):
+    '''
+    ARGUMENTS:
+        - fileName:        the name of the WAV file to be analyzed
+        - numOfSpeakers    the number of speakers (clusters) in the recording (<=0 for unknown)
+        - mtSize (opt)     mid-term window size
+        - mtStep (opt)     mid-term window step
+        - stWin  (opt)     short-term window size
+        - LDAdim (opt)     LDA dimension (0 for no LDA)
+        - PLOT     (opt)   0 for not plotting the results 1 for plottingy
+    '''
+    [Fs, x] = audioBasicIO.readAudioFile(fileName)
+    x = audioBasicIO.stereo2mono(x)
+    Duration = len(x) / Fs
+
+    [Classifier1, MEAN1, STD1, classNames1, mtWin1, mtStep1, stWin1, stStep1, computeBEAT1] = aT.loadKNNModel(os.path.join("data","knnSpeakerAll"))
+    [Classifier2, MEAN2, STD2, classNames2, mtWin2, mtStep2, stWin2, stStep2, computeBEAT2] = aT.loadKNNModel(os.path.join("data","knnSpeakerFemaleMale"))
+
+    [MidTermFeatures, ShortTermFeatures] = aF.mtFeatureExtraction(x, Fs, mtSize * Fs, mtStep * Fs, round(Fs * stWin), round(Fs*stWin * 0.5))
+
+    MidTermFeatures2 = numpy.zeros((MidTermFeatures.shape[0] + len(classNames1) + len(classNames2), MidTermFeatures.shape[1]))
+
+    for i in range(MidTermFeatures.shape[1]):
+        curF1 = (MidTermFeatures[:, i] - MEAN1) / STD1
+        curF2 = (MidTermFeatures[:, i] - MEAN2) / STD2
+        [Result, P1] = aT.classifierWrapper(Classifier1, "knn", curF1)
+        [Result, P2] = aT.classifierWrapper(Classifier2, "knn", curF2)
+        MidTermFeatures2[0:MidTermFeatures.shape[0], i] = MidTermFeatures[:, i]
+        MidTermFeatures2[MidTermFeatures.shape[0]:MidTermFeatures.shape[0]+len(classNames1), i] = P1 + 0.0001
+        MidTermFeatures2[MidTermFeatures.shape[0] + len(classNames1)::, i] = P2 + 0.0001
+
+    MidTermFeatures = MidTermFeatures2    # TODO
+    # SELECT FEATURES:
+    #iFeaturesSelect = [8,9,10,11,12,13,14,15,16,17,18,19,20];                                                                                         # SET 0A
+    #iFeaturesSelect = [8,9,10,11,12,13,14,15,16,17,18,19,20, 99,100];                                                                                 # SET 0B
+    #iFeaturesSelect = [8,9,10,11,12,13,14,15,16,17,18,19,20, 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,
+    #   97,98, 99,100];     # SET 0C
+
+    iFeaturesSelect = [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]                           # SET 1A
+    #iFeaturesSelect = [8,9,10,11,12,13,14,15,16,17,18,19,20,41,42,43,44,45,46,47,48,49,50,51,52,53, 99,100];                                          # SET 1B
+    #iFeaturesSelect = [8,9,10,11,12,13,14,15,16,17,18,19,20,41,42,43,44,45,46,47,48,49,50,51,52,53, 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98, 99,100];     # SET 1C
+
+    #iFeaturesSelect = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53];             # SET 2A
+    #iFeaturesSelect = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53, 99,100];     # SET 2B
+    #iFeaturesSelect = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53, 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98, 99,100];     # SET 2C
+
+    #iFeaturesSelect = range(100);                                                                                                    # SET 3
+    #MidTermFeatures += numpy.random.rand(MidTermFeatures.shape[0], MidTermFeatures.shape[1]) * 0.000000010
+
+    MidTermFeatures = MidTermFeatures[iFeaturesSelect, :]
+
+    (MidTermFeaturesNorm, MEAN, STD) = aT.normalizeFeatures([MidTermFeatures.T])
+    MidTermFeaturesNorm = MidTermFeaturesNorm[0].T
+    numOfWindows = MidTermFeatures.shape[1]
+
+    # remove outliers:
+    DistancesAll = numpy.sum(distance.squareform(distance.pdist(MidTermFeaturesNorm.T)), axis=0)
+    MDistancesAll = numpy.mean(DistancesAll)
+    iNonOutLiers = numpy.nonzero(DistancesAll < 1.2 * MDistancesAll)[0]
+
+    # TODO: Combine energy threshold for outlier removal:
+    #EnergyMin = numpy.min(MidTermFeatures[1,:])
+    #EnergyMean = numpy.mean(MidTermFeatures[1,:])
+    #Thres = (1.5*EnergyMin + 0.5*EnergyMean) / 2.0
+    #iNonOutLiers = numpy.nonzero(MidTermFeatures[1,:] > Thres)[0]
+    #print iNonOutLiers
+
+    perOutLier = (100.0 * (numOfWindows - iNonOutLiers.shape[0])) / numOfWindows
+    MidTermFeaturesNormOr = MidTermFeaturesNorm
+    MidTermFeaturesNorm = MidTermFeaturesNorm[:, iNonOutLiers]
+
+    # LDA dimensionality reduction:
+    if LDAdim > 0:
+        #[mtFeaturesToReduce, _] = aF.mtFeatureExtraction(x, Fs, mtSize * Fs, stWin * Fs, round(Fs*stWin), round(Fs*stWin));
+        # extract mid-term features with minimum step:
+        mtWinRatio = int(round(mtSize / stWin))
+        mtStepRatio = int(round(stWin / stWin))
+        mtFeaturesToReduce = []
+        numOfFeatures = len(ShortTermFeatures)
+        numOfStatistics = 2
+        #for i in range(numOfStatistics * numOfFeatures + 1):
+        for i in range(numOfStatistics * numOfFeatures):
+            mtFeaturesToReduce.append([])
+
+        for i in range(numOfFeatures):        # for each of the short-term features:
+            curPos = 0
+            N = len(ShortTermFeatures[i])
+            while (curPos < N):
+                N1 = curPos
+                N2 = curPos + mtWinRatio
+                if N2 > N:
+                    N2 = N
+                curStFeatures = ShortTermFeatures[i][N1:N2]
+                mtFeaturesToReduce[i].append(numpy.mean(curStFeatures))
+                mtFeaturesToReduce[i+numOfFeatures].append(numpy.std(curStFeatures))
+                curPos += mtStepRatio
+        mtFeaturesToReduce = numpy.array(mtFeaturesToReduce)
+        mtFeaturesToReduce2 = numpy.zeros((mtFeaturesToReduce.shape[0] + len(classNames1) + len(classNames2), mtFeaturesToReduce.shape[1]))
+        for i in range(mtFeaturesToReduce.shape[1]):
+            curF1 = (mtFeaturesToReduce[:, i] - MEAN1) / STD1
+            curF2 = (mtFeaturesToReduce[:, i] - MEAN2) / STD2
+            [Result, P1] = aT.classifierWrapper(Classifier1, "knn", curF1)
+            [Result, P2] = aT.classifierWrapper(Classifier2, "knn", curF2)
+            mtFeaturesToReduce2[0:mtFeaturesToReduce.shape[0], i] = mtFeaturesToReduce[:, i]
+            mtFeaturesToReduce2[mtFeaturesToReduce.shape[0]:mtFeaturesToReduce.shape[0] + len(classNames1), i] = P1 + 0.0001
+            mtFeaturesToReduce2[mtFeaturesToReduce.shape[0]+len(classNames1)::, i] = P2 + 0.0001
+        mtFeaturesToReduce = mtFeaturesToReduce2
+        mtFeaturesToReduce = mtFeaturesToReduce[iFeaturesSelect, :]
+        #mtFeaturesToReduce += numpy.random.rand(mtFeaturesToReduce.shape[0], mtFeaturesToReduce.shape[1]) * 0.0000010
+        (mtFeaturesToReduce, MEAN, STD) = aT.normalizeFeatures([mtFeaturesToReduce.T])
+        mtFeaturesToReduce = mtFeaturesToReduce[0].T
+        #DistancesAll = numpy.sum(distance.squareform(distance.pdist(mtFeaturesToReduce.T)), axis=0)
+        #MDistancesAll = numpy.mean(DistancesAll)
+        #iNonOutLiers2 = numpy.nonzero(DistancesAll < 3.0*MDistancesAll)[0]
+        #mtFeaturesToReduce = mtFeaturesToReduce[:, iNonOutLiers2]
+        Labels = numpy.zeros((mtFeaturesToReduce.shape[1], ));
+        LDAstep = 1.0
+        LDAstepRatio = LDAstep / stWin
+        #print LDAstep, LDAstepRatio
+        for i in range(Labels.shape[0]):
+            Labels[i] = int(i*stWin/LDAstepRatio);        
+        clf = sklearn.discriminant_analysis.LinearDiscriminantAnalysis(n_components=LDAdim)
+        clf.fit(mtFeaturesToReduce.T, Labels)
+        MidTermFeaturesNorm = (clf.transform(MidTermFeaturesNorm.T)).T
+
+    if numOfSpeakers <= 0:
+        sRange = range(2, 10)
+    else:
+        sRange = [numOfSpeakers]
+    clsAll = []
+    silAll = []
+    centersAll = []
+    
+    for iSpeakers in sRange:        
+        k_means = sklearn.cluster.KMeans(n_clusters = iSpeakers)
+        k_means.fit(MidTermFeaturesNorm.T)
+        cls = k_means.labels_        
+        means = k_means.cluster_centers_
+
+        # Y = distance.squareform(distance.pdist(MidTermFeaturesNorm.T))
+        clsAll.append(cls)
+        centersAll.append(means)
+        silA = []; silB = []
+        for c in range(iSpeakers):                                # for each speaker (i.e. for each extracted cluster)
+            clusterPerCent = numpy.nonzero(cls==c)[0].shape[0] / float(len(cls))
+            if clusterPerCent < 0.020:
+                silA.append(0.0)
+                silB.append(0.0)
+            else:
+                MidTermFeaturesNormTemp = MidTermFeaturesNorm[:,cls==c]            # get subset of feature vectors
+                Yt = distance.pdist(MidTermFeaturesNormTemp.T)                # compute average distance between samples that belong to the cluster (a values)
+                silA.append(numpy.mean(Yt)*clusterPerCent)
+                silBs = []
+                for c2 in range(iSpeakers):                        # compute distances from samples of other clusters
+                    if c2!=c:
+                        clusterPerCent2 = numpy.nonzero(cls==c2)[0].shape[0] / float(len(cls))
+                        MidTermFeaturesNormTemp2 = MidTermFeaturesNorm[:,cls==c2]
+                        Yt = distance.cdist(MidTermFeaturesNormTemp.T, MidTermFeaturesNormTemp2.T)
+                        silBs.append(numpy.mean(Yt)*(clusterPerCent+clusterPerCent2)/2.0)
+                silBs = numpy.array(silBs)                            
+                silB.append(min(silBs))                            # ... and keep the minimum value (i.e. the distance from the "nearest" cluster)
+        silA = numpy.array(silA); 
+        silB = numpy.array(silB); 
+        sil = []
+        for c in range(iSpeakers):                                # for each cluster (speaker)
+            sil.append( ( silB[c] - silA[c]) / (max(silB[c],  silA[c])+0.00001)  )        # compute silhouette
+
+        silAll.append(numpy.mean(sil))                                # keep the AVERAGE SILLOUETTE
+
+    #silAll = silAll * (1.0/(numpy.power(numpy.array(sRange),0.5)))
+    imax = numpy.argmax(silAll)                                    # position of the maximum sillouette value
+    nSpeakersFinal = sRange[imax]                                    # optimal number of clusters
+
+    # generate the final set of cluster labels
+    # (important: need to retrieve the outlier windows: this is achieved by giving them the value of their nearest non-outlier window)
+    cls = numpy.zeros((numOfWindows,))
+    for i in range(numOfWindows):
+        j = numpy.argmin(numpy.abs(i-iNonOutLiers))        
+        cls[i] = clsAll[imax][j]
+        
+    # Post-process method 1: hmm smoothing
+    for i in range(1):
+        startprob, transmat, means, cov = trainHMM_computeStatistics(MidTermFeaturesNormOr, cls)
+        hmm = hmmlearn.hmm.GaussianHMM(startprob.shape[0], "diag")            # hmm training        
+        hmm.startprob_ = startprob
+        hmm.transmat_ = transmat            
+        hmm.means_ = means; hmm.covars_ = cov
+        cls = hmm.predict(MidTermFeaturesNormOr.T)                    
+    
+    # Post-process method 2: median filtering:
+    cls = scipy.signal.medfilt(cls, 13)
+    cls = scipy.signal.medfilt(cls, 11)
+
+    sil = silAll[imax]                                        # final sillouette
+    classNames = ["speaker{0:d}".format(c) for c in range(nSpeakersFinal)];
+
+
+    # load ground-truth if available
+    gtFile = fileName.replace('.wav', '.segments');                            # open for annotated file
+    if os.path.isfile(gtFile):                                    # if groundturh exists
+        [segStart, segEnd, segLabels] = readSegmentGT(gtFile)                    # read GT data
+        flagsGT, classNamesGT = segs2flags(segStart, segEnd, segLabels, mtStep)            # convert to flags
+
+    if PLOT:
+        fig = plt.figure()    
+        if numOfSpeakers>0:
+            ax1 = fig.add_subplot(111)
+        else:
+            ax1 = fig.add_subplot(211)
+        ax1.set_yticks(numpy.array(range(len(classNames))))
+        ax1.axis((0, Duration, -1, len(classNames)))
+        ax1.set_yticklabels(classNames)
+        ax1.plot(numpy.array(range(len(cls)))*mtStep+mtStep/2.0, cls)
+
+    if os.path.isfile(gtFile):
+        if PLOT:
+            ax1.plot(numpy.array(range(len(flagsGT)))*mtStep+mtStep/2.0, flagsGT, 'r')
+        purityClusterMean, puritySpeakerMean = evaluateSpeakerDiarization(cls, flagsGT)
+        print "{0:.1f}\t{1:.1f}".format(100*purityClusterMean, 100*puritySpeakerMean)
+        if PLOT:
+            plt.title("Cluster purity: {0:.1f}% - Speaker purity: {1:.1f}%".format(100*purityClusterMean, 100*puritySpeakerMean) )
+    if PLOT:
+        plt.xlabel("time (seconds)")
+        #print sRange, silAll    
+        if numOfSpeakers<=0:
+            plt.subplot(212)
+            plt.plot(sRange, silAll)
+            plt.xlabel("number of clusters");
+            plt.ylabel("average clustering's sillouette");
+        plt.show()
+    return cls
+    
+def speakerDiarizationEvaluateScript(folderName, LDAs):
+    '''
+        This function prints the cluster purity and speaker purity for each WAV file stored in a provided directory (.SEGMENT files are needed as ground-truth)
+        ARGUMENTS:
+            - folderName:     the full path of the folder where the WAV and SEGMENT (ground-truth) files are stored
+            - LDAs:            a list of LDA dimensions (0 for no LDA)
+    '''
+    types = ('*.wav',  )
+    wavFilesList = []
+    for files in types:
+        wavFilesList.extend(glob.glob(os.path.join(folderName, files)))    
+    
+    wavFilesList = sorted(wavFilesList)
+
+    # get number of unique speakers per file (from ground-truth)    
+    N = []
+    for wavFile in wavFilesList:        
+        gtFile = wavFile.replace('.wav', '.segments');
+        if os.path.isfile(gtFile):
+            [segStart, segEnd, segLabels] = readSegmentGT(gtFile)                            # read GT data
+            N.append(len(list(set(segLabels))))
+        else:
+            N.append(-1)
+    
+    for l in LDAs:
+        print "LDA = {0:d}".format(l)
+        for i, wavFile in enumerate(wavFilesList):
+            speakerDiarization(wavFile, N[i], 2.0, 0.2, 0.05, l, PLOT = False)            
+        print
+        
+def musicThumbnailing(x, Fs, shortTermSize=1.0, shortTermStep=0.5, thumbnailSize=10.0, Limit1 = 0, Limit2 = 1):
+    '''
+    This function detects instances of the most representative part of a music recording, also called "music thumbnails".
+    A technique similar to the one proposed in [1], however a wider set of audio features is used instead of chroma features.
+    In particular the following steps are followed:
+     - Extract short-term audio features. Typical short-term window size: 1 second
+     - Compute the self-silimarity matrix, i.e. all pairwise similarities between feature vectors
+      - Apply a diagonal mask is as a moving average filter on the values of the self-similarty matrix. 
+       The size of the mask is equal to the desirable thumbnail length.
+      - Find the position of the maximum value of the new (filtered) self-similarity matrix.
+       The audio segments that correspond to the diagonial around that position are the selected thumbnails
+    
+
+    ARGUMENTS:
+     - x:            input signal
+     - Fs:            sampling frequency
+     - shortTermSize:     window size (in seconds)
+     - shortTermStep:    window step (in seconds)
+     - thumbnailSize:    desider thumbnail size (in seconds)
+    
+    RETURNS:
+     - A1:            beginning of 1st thumbnail (in seconds)
+     - A2:            ending of 1st thumbnail (in seconds)
+     - B1:            beginning of 2nd thumbnail (in seconds)
+     - B2:            ending of 2nd thumbnail (in seconds)
+
+    USAGE EXAMPLE:
+       import audioFeatureExtraction as aF
+     [Fs, x] = basicIO.readAudioFile(inputFile)
+     [A1, A2, B1, B2] = musicThumbnailing(x, Fs)
+
+    [1] Bartsch, M. A., & Wakefield, G. H. (2005). Audio thumbnailing of popular music using chroma-based representations. 
+    Multimedia, IEEE Transactions on, 7(1), 96-104.
+    '''
+    x = audioBasicIO.stereo2mono(x);
+    # feature extraction:
+    stFeatures = aF.stFeatureExtraction(x, Fs, Fs*shortTermSize, Fs*shortTermStep)
+
+    # self-similarity matrix
+    S = selfSimilarityMatrix(stFeatures)
+
+    # moving filter:
+    M = int(round(thumbnailSize / shortTermStep))
+    B = numpy.eye(M,M)
+    S = scipy.signal.convolve2d(S, B, 'valid')
+
+
+    # post-processing (remove main diagonal elements)
+    MIN = numpy.min(S)
+    for i in range(S.shape[0]):
+        for j in range(S.shape[1]):
+            if abs(i-j) < 5.0 / shortTermStep or i > j:
+                S[i,j] = MIN;
+
+    # find max position:
+    S[0:int(Limit1*S.shape[0]), :] = MIN
+    S[:, 0:int(Limit1*S.shape[0])] = MIN
+    S[int(Limit2*S.shape[0])::, :] = MIN
+    S[:, int(Limit2*S.shape[0])::] = MIN
+
+    maxVal = numpy.max(S)        
+    [I, J] = numpy.unravel_index(S.argmax(), S.shape)
+    #plt.imshow(S)
+    #plt.show()
+    # expand:
+    i1 = I; i2 = I
+    j1 = J; j2 = J
+
+    while i2-i1<M: 
+        if i1 <=0 or j1<=0 or i2>=S.shape[0]-2 or j2>=S.shape[1]-2:       
+            break
+        if S[i1-1, j1-1] > S[i2+1,j2+1]:            
+            i1 -= 1
+            j1 -= 1            
+        else:            
+            i2 += 1
+            j2 += 1            
+
+    return (shortTermStep*i1, shortTermStep*i2, shortTermStep*j1, shortTermStep*j2, S)
+
+
diff --git a/audioTrainTest.py b/audioTrainTest.py
new file mode 100644
index 0000000..cd91673
--- /dev/null
+++ b/audioTrainTest.py
@@ -0,0 +1,1141 @@
+import sys
+import numpy
+import os
+import glob
+import cPickle
+import signal
+import csv
+import ntpath
+import audioFeatureExtraction as aF
+import audioBasicIO
+from scipy import linalg as la
+from scipy.spatial import distance
+import sklearn.svm
+import sklearn.decomposition
+import sklearn.ensemble
+
+def signal_handler(signal, frame):
+    print 'You pressed Ctrl+C! - EXIT'
+    os.system("stty -cbreak echo")
+    sys.exit(0)
+
+signal.signal(signal.SIGINT, signal_handler)
+
+shortTermWindow = 0.050
+shortTermStep = 0.050
+eps = 0.00000001
+
+
+class kNN:
+    def __init__(self, X, Y, k):
+        self.X = X
+        self.Y = Y
+        self.k = k
+
+    def classify(self, testSample):
+        nClasses = numpy.unique(self.Y).shape[0]
+        YDist = (distance.cdist(self.X, testSample.reshape(1, testSample.shape[0]), 'euclidean')).T
+        iSort = numpy.argsort(YDist)
+        P = numpy.zeros((nClasses,))
+        for i in range(nClasses):
+            P[i] = numpy.nonzero(self.Y[iSort[0][0:self.k]] == i)[0].shape[0] / float(self.k)
+        return (numpy.argmax(P), P)
+
+
+def classifierWrapper(classifier, classifierType, testSample):
+    '''
+    This function is used as a wrapper to pattern classification.
+    ARGUMENTS:
+        - classifier:        a classifier object of type sklearn.svm.SVC or kNN (defined in this library) or sklearn.ensemble.RandomForestClassifier or sklearn.ensemble.GradientBoostingClassifier  or sklearn.ensemble.ExtraTreesClassifier
+        - classifierType:    "svm" or "knn" or "randomforests" or "gradientboosting" or "extratrees"
+        - testSample:        a feature vector (numpy array)
+    RETURNS:
+        - R:            class ID
+        - P:            probability estimate
+
+    EXAMPLE (for some audio signal stored in array x):
+        import audioFeatureExtraction as aF
+        import audioTrainTest as aT
+        # load the classifier (here SVM, for kNN use loadKNNModel instead):
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep] = aT.loadSVModel(modelName)
+        # mid-term feature extraction:
+        [MidTermFeatures, _] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs*stWin), round(Fs*stStep));
+        # feature normalization:
+        curFV = (MidTermFeatures[:, i] - MEAN) / STD;
+        # classification
+        [Result, P] = classifierWrapper(Classifier, modelType, curFV)
+    '''
+    R = -1
+    P = -1
+    if classifierType == "knn":
+        [R, P] = classifier.classify(testSample)
+    elif classifierType == "svm" or classifierType == "randomforest" or classifierType == "gradientboosting" or "extratrees":
+        R = classifier.predict(testSample.reshape(1,-1))[0]
+        P = classifier.predict_proba(testSample.reshape(1,-1))[0]
+    return [R, P]
+
+
+def regressionWrapper(model, modelType, testSample):
+    '''
+    This function is used as a wrapper to pattern classification.
+    ARGUMENTS:
+        - model:        regression model
+        - modelType:        "svm" or "knn" (TODO)
+        - testSample:        a feature vector (numpy array)
+    RETURNS:
+        - R:            regression result (estimated value)
+
+    EXAMPLE (for some audio signal stored in array x):
+        TODO
+    '''
+    if modelType == "svm" or modelType == "randomforest" or modelType == "svm_rbf":
+        return (model.predict(testSample.reshape(1,-1))[0])
+
+    #    elif classifierType == "knn":
+    #    TODO
+
+    return None
+
+
+def randSplitFeatures(features, partTrain):
+    '''
+    def randSplitFeatures(features):
+
+    This function splits a feature set for training and testing.
+
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features.
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - partTrain:        percentage
+    RETURNS:
+        - featuresTrains:    a list of training data for each class
+        - featuresTest:        a list of testing data for each class
+    '''
+
+    featuresTrain = []
+    featuresTest = []
+    for i, f in enumerate(features):
+        [numOfSamples, numOfDims] = f.shape
+        randperm = numpy.random.permutation(range(numOfSamples))
+        nTrainSamples = int(round(partTrain * numOfSamples))
+        featuresTrain.append(f[randperm[0:nTrainSamples]])
+        featuresTest.append(f[randperm[nTrainSamples::]])
+    return (featuresTrain, featuresTest)
+
+
+def trainKNN(features, K):
+    '''
+    Train a kNN  classifier.
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features.
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - K:                parameter K
+    RETURNS:
+        - kNN:              the trained kNN variable
+
+    '''
+    [Xt, Yt] = listOfFeatures2Matrix(features)
+    knn = kNN(Xt, Yt, K)
+    return knn
+
+
+def trainSVM(features, Cparam):
+    '''
+    Train a multi-class probabilitistic SVM classifier.
+    Note:     This function is simply a wrapper to the sklearn functionality for SVM training
+              See function trainSVM_feature() to use a wrapper on both the feature extraction and the SVM training (and parameter tuning) processes.
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - Cparam:           SVM parameter C (cost of constraints violation)
+    RETURNS:
+        - svm:              the trained SVM variable
+
+    NOTE:
+        This function trains a linear-kernel SVM for a given C value. For a different kernel, other types of parameters should be provided.
+    '''
+
+    [X, Y] = listOfFeatures2Matrix(features)
+    svm = sklearn.svm.SVC(C = Cparam, kernel = 'linear',  probability = True)        
+    svm.fit(X,Y)
+
+    return svm
+
+def trainSVM_RBF(features, Cparam):
+    '''
+    Train a multi-class probabilitistic SVM classifier.
+    Note:     This function is simply a wrapper to the sklearn functionality for SVM training
+              See function trainSVM_feature() to use a wrapper on both the feature extraction and the SVM training (and parameter tuning) processes.
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - Cparam:           SVM parameter C (cost of constraints violation)
+    RETURNS:
+        - svm:              the trained SVM variable
+
+    NOTE:
+        This function trains a linear-kernel SVM for a given C value. For a different kernel, other types of parameters should be provided.
+    '''
+
+    [X, Y] = listOfFeatures2Matrix(features)
+    svm = sklearn.svm.SVC(C = Cparam, kernel = 'rbf',  probability = True)        
+    svm.fit(X,Y)
+
+    return svm
+
+
+def trainRandomForest(features, n_estimators):
+    '''
+    Train a multi-class decision tree classifier.
+    Note:     This function is simply a wrapper to the sklearn functionality for SVM training
+              See function trainSVM_feature() to use a wrapper on both the feature extraction and the SVM training (and parameter tuning) processes.
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - n_estimators:     number of trees in the forest
+    RETURNS:
+        - svm:              the trained SVM variable
+
+    NOTE:
+        This function trains a linear-kernel SVM for a given C value. For a different kernel, other types of parameters should be provided.
+    '''
+
+    [X, Y] = listOfFeatures2Matrix(features)
+    rf = sklearn.ensemble.RandomForestClassifier(n_estimators = n_estimators)
+    rf.fit(X,Y)
+
+    return rf
+
+def trainGradientBoosting(features, n_estimators):
+    '''
+    Train a gradient boosting classifier
+    Note:     This function is simply a wrapper to the sklearn functionality for SVM training
+              See function trainSVM_feature() to use a wrapper on both the feature extraction and the SVM training (and parameter tuning) processes.
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - n_estimators:     number of trees in the forest
+    RETURNS:
+        - svm:              the trained SVM variable
+
+    NOTE:
+        This function trains a linear-kernel SVM for a given C value. For a different kernel, other types of parameters should be provided.
+    '''
+
+    [X, Y] = listOfFeatures2Matrix(features)
+    rf = sklearn.ensemble.GradientBoostingClassifier(n_estimators = n_estimators)
+    rf.fit(X,Y)
+
+    return rf
+
+def trainExtraTrees(features, n_estimators):
+    '''
+    Train a gradient boosting classifier
+    Note:     This function is simply a wrapper to the sklearn functionality for extra tree classifiers
+              See function trainSVM_feature() to use a wrapper on both the feature extraction and the SVM training (and parameter tuning) processes.
+    ARGUMENTS:
+        - features:         a list ([numOfClasses x 1]) whose elements containt numpy matrices of features
+                            each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        - n_estimators:     number of trees in the forest
+    RETURNS:
+        - svm:              the trained SVM variable
+
+    NOTE:
+        This function trains a linear-kernel SVM for a given C value. For a different kernel, other types of parameters should be provided.
+    '''
+
+    [X, Y] = listOfFeatures2Matrix(features)
+    et = sklearn.ensemble.ExtraTreesClassifier(n_estimators = n_estimators)
+    et.fit(X,Y)
+
+    return et
+
+
+def trainSVMregression(Features, Y, Cparam):    
+    svm = sklearn.svm.SVR(C = Cparam, kernel = 'linear')    
+    svm.fit(Features,Y)    
+    trainError = numpy.mean(numpy.abs(svm.predict(Features) - Y))
+    return svm, trainError
+
+
+def trainSVMregression_rbf(Features, Y, Cparam):    
+    svm = sklearn.svm.SVR(C = Cparam, kernel = 'rbf')    
+    svm.fit(Features,Y)    
+    trainError = numpy.mean(numpy.abs(svm.predict(Features) - Y))
+    return svm, trainError
+
+
+def trainRandomForestRegression(Features, Y, n_estimators):    
+    rf = sklearn.ensemble.RandomForestRegressor(n_estimators = n_estimators)
+    rf.fit(Features,Y)
+    trainError = numpy.mean(numpy.abs(rf.predict(Features) - Y))
+    return rf, trainError
+
+def featureAndTrain(listOfDirs, mtWin, mtStep, stWin, stStep, classifierType, modelName, computeBEAT=False, perTrain=0.90):
+    '''
+    This function is used as a wrapper to segment-based audio feature extraction and classifier training.
+    ARGUMENTS:
+        listOfDirs:        list of paths of directories. Each directory contains a signle audio class whose samples are stored in seperate WAV files.
+        mtWin, mtStep:        mid-term window length and step
+        stWin, stStep:        short-term window and step
+        classifierType:        "svm" or "knn" or "randomforest" or "gradientboosting" or "extratrees"
+        modelName:        name of the model to be saved
+    RETURNS:
+        None. Resulting classifier along with the respective model parameters are saved on files.
+    '''
+
+    # STEP A: Feature Extraction:
+    [features, classNames, _] = aF.dirsWavFeatureExtraction(listOfDirs, mtWin, mtStep, stWin, stStep, computeBEAT=computeBEAT)
+
+    if len(features) == 0:
+        print "trainSVM_feature ERROR: No data found in any input folder!"
+        return
+
+    numOfFeatures = features[0].shape[1]
+    featureNames = ["features" + str(d + 1) for d in range(numOfFeatures)]
+
+    writeTrainDataToARFF(modelName, features, classNames, featureNames)
+
+    for i, f in enumerate(features):
+        if len(f) == 0:
+            print "trainSVM_feature ERROR: " + listOfDirs[i] + " folder is empty or non-existing!"
+            return
+
+    # STEP B: Classifier Evaluation and Parameter Selection:
+    if classifierType == "svm" or classifierType == "svm_rbf":
+        classifierParams = numpy.array([0.001, 0.01,  0.5, 1.0, 5.0, 10.0, 20.0])
+    elif classifierType == "randomforest":
+        classifierParams = numpy.array([10, 25, 50, 100,200,500])
+    elif classifierType == "knn":
+        classifierParams = numpy.array([1, 3, 5, 7, 9, 11, 13, 15])        
+    elif classifierType == "gradientboosting":
+        classifierParams = numpy.array([10, 25, 50, 100,200,500])        
+    elif classifierType == "extratrees":
+        classifierParams = numpy.array([10, 25, 50, 100,200,500])        
+
+    # get optimal classifeir parameter:
+    features2 = []
+    for f in features:        
+        fTemp = []
+        for i in range(f.shape[0]):
+            temp = f[i,:]
+            if (not numpy.isnan(temp).any()) and (not numpy.isinf(temp).any()) :
+                fTemp.append(temp.tolist())
+            else:
+                print "NaN Found! Feature vector not used for training"
+        features2.append(numpy.array(fTemp))
+    features = features2
+
+    bestParam = evaluateClassifier(features, classNames, 100, classifierType, classifierParams, 0, perTrain)
+
+    print "Selected params: {0:.5f}".format(bestParam)
+
+    C = len(classNames)
+    [featuresNorm, MEAN, STD] = normalizeFeatures(features)        # normalize features
+    MEAN = MEAN.tolist()
+    STD = STD.tolist()
+    featuresNew = featuresNorm
+
+    # STEP C: Save the classifier to file
+    if classifierType == "svm":
+        Classifier = trainSVM(featuresNew, bestParam)        
+    elif classifierType == "svm_rbf":
+        Classifier = trainSVM_RBF(featuresNew, bestParam)
+    elif classifierType == "randomforest":
+        Classifier = trainRandomForest(featuresNew, bestParam)
+    elif classifierType == "gradientboosting":
+        Classifier = trainGradientBoosting(featuresNew, bestParam)
+    elif classifierType == "extratrees":
+        Classifier = trainExtraTrees(featuresNew, bestParam)
+
+
+    if classifierType == "knn":
+        [X, Y] = listOfFeatures2Matrix(featuresNew)
+        X = X.tolist()
+        Y = Y.tolist()
+        fo = open(modelName, "wb")
+        cPickle.dump(X, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(Y,  fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(MEAN, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(STD,  fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(classNames,  fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(bestParam,  fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(mtWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(mtStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(stWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(stStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(computeBEAT, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        fo.close()
+    elif classifierType == "svm" or classifierType == "svm_rbf" or classifierType == "randomforest" or classifierType == "gradientboosting" or classifierType == "extratrees":
+        with open(modelName, 'wb') as fid:                                            # save to file
+            cPickle.dump(Classifier, fid)            
+        fo = open(modelName + "MEANS", "wb")
+        cPickle.dump(MEAN, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(STD, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(classNames, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(mtWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(mtStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(stWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(stStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        cPickle.dump(computeBEAT, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+        fo.close()        
+
+
+def featureAndTrainRegression(dirName, mtWin, mtStep, stWin, stStep, modelType, modelName, computeBEAT=False):
+    '''
+    This function is used as a wrapper to segment-based audio feature extraction and classifier training.
+    ARGUMENTS:
+        dirName:        path of directory containing the WAV files and Regression CSVs
+        mtWin, mtStep:        mid-term window length and step
+        stWin, stStep:        short-term window and step
+        modelType:        "svm" or "knn" or "randomforest"
+        modelName:        name of the model to be saved
+    RETURNS:
+        None. Resulting regression model along with the respective model parameters are saved on files.
+    '''
+    # STEP A: Feature Extraction:
+    [features, _, fileNames] = aF.dirsWavFeatureExtraction([dirName], mtWin, mtStep, stWin, stStep, computeBEAT=computeBEAT)
+    features = features[0]
+    fileNames = [ntpath.basename(f) for f in fileNames[0]]
+    featuresFinal = []
+
+    # Read CSVs:
+    CSVs = glob.glob(dirName + os.sep + "*.csv")
+    regressionLabels = []
+    regressionNames = []
+    featuresFinal = []
+    for c in CSVs:                                                            # for each CSV
+        #curRegressionLabels = numpy.zeros((len(fileNames, )))                 # read filenames, map to "fileNames" and append respective values in the regressionLabels
+        curRegressionLabels = []
+        featuresTemp = []
+        with open(c, 'rb') as csvfile:                                        # open the csv file that contains the current target value's annotations
+            CSVreader = csv.reader(csvfile, delimiter=',', quotechar='|')
+            for row in CSVreader:
+                if len(row) == 2:                                             # if the current row contains two fields (filename, target value)
+                    if row[0] in fileNames:                                   # ... and if the current filename exists in the list of filenames
+                        index = fileNames.index(row[0])
+                        #curRegressionLabels[index] = float(row[1])
+                        curRegressionLabels.append(float(row[1]))
+                        featuresTemp.append(features[index,:])
+
+        featuresFinal.append(numpy.array(featuresTemp))
+        regressionLabels.append(numpy.array(curRegressionLabels))                          # curRegressionLabels is the list of values for the current regression problem
+        regressionNames.append(ntpath.basename(c).replace(".csv", ""))        # regression task name   
+        if len(features) == 0:
+            print "ERROR: No data found in any input folder!"
+            return
+
+    numOfFeatures = featuresFinal[0].shape[1]
+
+    # TODO: ARRF WRITE????
+    # STEP B: Classifier Evaluation and Parameter Selection:
+    if modelType == "svm" or modelType == "svm_rbf":
+        modelParams = numpy.array([0.001, 0.005, 0.01, 0.05, 0.1, 0.25, 0.5, 1.0, 5.0, 10.0])        
+    elif modelType == "randomforest":
+        modelParams = numpy.array([5, 10, 25, 50, 100])
+
+#    elif modelType == "knn":
+#        modelParams = numpy.array([1, 3, 5, 7, 9, 11, 13, 15]);
+    errors = []
+    errorsBase = []
+    bestParams = []
+
+    for iRegression, r in enumerate(regressionNames):
+        # get optimal classifeir parameter:
+        print "Regression task " + r
+        bestParam, error, berror = evaluateRegression(featuresFinal[iRegression], regressionLabels[iRegression], 100, modelType, modelParams)
+        errors.append(error)
+        errorsBase.append(berror)
+        bestParams.append(bestParam)
+        print "Selected params: {0:.5f}".format(bestParam)
+
+        [featuresNorm, MEAN, STD] = normalizeFeatures([featuresFinal[iRegression]])        # normalize features
+
+        # STEP C: Save the model to file
+        if modelType == "svm":
+            Classifier, _ = trainSVMregression(featuresNorm[0], regressionLabels[iRegression], bestParam)
+        if modelType == "svm_rbf":
+            Classifier, _ = trainSVMregression_rbf(featuresNorm[0], regressionLabels[iRegression], bestParam)
+        if modelType == "randomforest":
+            Classifier, _ = trainRandomForestRegression(featuresNorm[0], regressionLabels[iRegression], bestParam)
+
+        if modelType == "svm" or modelType == "svm_rbf" or modelType == "randomforest":
+            with open(modelName + "_" + r, 'wb') as fid:                                            # save to file
+                cPickle.dump(Classifier, fid)            
+            fo = open(modelName + "_" + r + "MEANS", "wb")
+            cPickle.dump(MEAN, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            cPickle.dump(STD,  fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            cPickle.dump(mtWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            cPickle.dump(mtStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            cPickle.dump(stWin, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            cPickle.dump(stStep, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            cPickle.dump(computeBEAT, fo, protocol=cPickle.HIGHEST_PROTOCOL)
+            fo.close()
+    return errors, errorsBase, bestParams
+
+
+def loadKNNModel(kNNModelName, isRegression=False):
+    try:
+        fo = open(kNNModelName, "rb")
+    except IOError:
+        print "didn't find file"
+        return
+    try:
+        X = cPickle.load(fo)
+        Y = cPickle.load(fo)
+        MEAN = cPickle.load(fo)
+        STD = cPickle.load(fo)
+        if not isRegression:
+            classNames = cPickle.load(fo)
+        K = cPickle.load(fo)
+        mtWin = cPickle.load(fo)
+        mtStep = cPickle.load(fo)
+        stWin = cPickle.load(fo)
+        stStep = cPickle.load(fo)
+        computeBEAT = cPickle.load(fo)
+    except:
+        fo.close()
+    fo.close()
+
+    X = numpy.array(X)
+    Y = numpy.array(Y)
+    MEAN = numpy.array(MEAN)
+    STD = numpy.array(STD)
+
+    Classifier = kNN(X, Y, K)  # Note: a direct call to the kNN constructor is used here
+
+    if isRegression:
+        return(Classifier, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT)
+    else:
+        return(Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT)
+
+
+def loadSVModel(SVMmodelName, isRegression=False):
+    '''
+    This function loads an SVM model either for classification or training.
+    ARGMUMENTS:
+        - SVMmodelName:     the path of the model to be loaded
+        - isRegression:        a flag indigating whereas this model is regression or not
+    '''
+    try:
+        fo = open(SVMmodelName+"MEANS", "rb")
+    except IOError:
+            print "Load SVM Model: Didn't find file"
+            return
+    try:
+        MEAN = cPickle.load(fo)
+        STD = cPickle.load(fo)
+        if not isRegression:
+            classNames = cPickle.load(fo)
+        mtWin = cPickle.load(fo)
+        mtStep = cPickle.load(fo)
+        stWin = cPickle.load(fo)
+        stStep = cPickle.load(fo)
+        computeBEAT = cPickle.load(fo)
+
+    except:
+        fo.close()
+    fo.close()
+
+    MEAN = numpy.array(MEAN)
+    STD = numpy.array(STD)
+
+    COEFF = []
+    with open(SVMmodelName, 'rb') as fid:
+        SVM = cPickle.load(fid)    
+
+    if isRegression:
+        return(SVM, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT)
+    else:
+        return(SVM, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT)
+
+
+def loadRandomForestModel(RFmodelName, isRegression=False):
+    '''
+    This function loads an SVM model either for classification or training.
+    ARGMUMENTS:
+        - SVMmodelName:     the path of the model to be loaded
+        - isRegression:     a flag indigating whereas this model is regression or not
+    '''
+    try:
+        fo = open(RFmodelName+"MEANS", "rb")
+    except IOError:
+            print "Load Random Forest Model: Didn't find file"
+            return
+    try:
+        MEAN = cPickle.load(fo)
+        STD = cPickle.load(fo)
+        if not isRegression:
+            classNames = cPickle.load(fo)
+        mtWin = cPickle.load(fo)
+        mtStep = cPickle.load(fo)
+        stWin = cPickle.load(fo)
+        stStep = cPickle.load(fo)
+        computeBEAT = cPickle.load(fo)
+
+    except:
+        fo.close()
+    fo.close()
+
+    MEAN = numpy.array(MEAN)
+    STD = numpy.array(STD)
+
+    COEFF = []
+    with open(RFmodelName, 'rb') as fid:
+        RF = cPickle.load(fid)    
+
+    if isRegression:
+        return(RF, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT)
+    else:
+        return(RF, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT)
+
+def loadGradientBoostingModel(GBModelName, isRegression=False):
+    '''
+    This function loads gradient boosting either for classification or training.
+    ARGMUMENTS:
+        - SVMmodelName:     the path of the model to be loaded
+        - isRegression:     a flag indigating whereas this model is regression or not
+    '''
+    try:
+        fo = open(GBModelName+"MEANS", "rb")
+    except IOError:
+            print "Load Random Forest Model: Didn't find file"
+            return
+    try:
+        MEAN = cPickle.load(fo)
+        STD = cPickle.load(fo)
+        if not isRegression:
+            classNames = cPickle.load(fo)
+        mtWin = cPickle.load(fo)
+        mtStep = cPickle.load(fo)
+        stWin = cPickle.load(fo)
+        stStep = cPickle.load(fo)
+        computeBEAT = cPickle.load(fo)
+
+    except:
+        fo.close()
+    fo.close()
+
+    MEAN = numpy.array(MEAN)
+    STD = numpy.array(STD)
+
+    COEFF = []
+    with open(GBModelName, 'rb') as fid:
+        GB = cPickle.load(fid)    
+
+    if isRegression:
+        return(GB, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT)
+    else:
+        return(GB, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT)
+
+def loadExtraTreesModel(ETmodelName, isRegression=False):
+    '''
+    This function loads extra trees either for classification or training.
+    ARGMUMENTS:
+        - SVMmodelName:     the path of the model to be loaded
+        - isRegression:     a flag indigating whereas this model is regression or not
+    '''
+    try:
+        fo = open(ETmodelName+"MEANS", "rb")
+    except IOError:
+            print "Load Random Forest Model: Didn't find file"
+            return
+    try:
+        MEAN = cPickle.load(fo)
+        STD = cPickle.load(fo)
+        if not isRegression:
+            classNames = cPickle.load(fo)
+        mtWin = cPickle.load(fo)
+        mtStep = cPickle.load(fo)
+        stWin = cPickle.load(fo)
+        stStep = cPickle.load(fo)
+        computeBEAT = cPickle.load(fo)
+
+    except:
+        fo.close()
+    fo.close()
+
+    MEAN = numpy.array(MEAN)
+    STD = numpy.array(STD)
+
+    COEFF = []
+    with open(ETmodelName, 'rb') as fid:
+        GB = cPickle.load(fid)    
+
+    if isRegression:
+        return(GB, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT)
+    else:
+        return(GB, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT)
+
+
+def evaluateClassifier(features, ClassNames, nExp, ClassifierName, Params, parameterMode, perTrain=0.90):
+    '''
+    ARGUMENTS:
+        features:     a list ([numOfClasses x 1]) whose elements containt numpy matrices of features.
+                each matrix features[i] of class i is [numOfSamples x numOfDimensions]
+        ClassNames:    list of class names (strings)
+        nExp:        number of cross-validation experiments
+        ClassifierName: svm or knn or randomforest
+        Params:        list of classifier parameters (for parameter tuning during cross-validation)
+        parameterMode:    0: choose parameters that lead to maximum overall classification ACCURACY
+                1: choose parameters that lead to maximum overall F1 MEASURE
+    RETURNS:
+         bestParam:    the value of the input parameter that optimizes the selected performance measure
+    '''
+
+    # feature normalization:
+    (featuresNorm, MEAN, STD) = normalizeFeatures(features)
+    #featuresNorm = features;
+    nClasses = len(features)
+    CAll = []
+    acAll = []
+    F1All = []
+    PrecisionClassesAll = []
+    RecallClassesAll = []
+    ClassesAll = []
+    F1ClassesAll = []
+    CMsAll = []
+
+    # compute total number of samples:
+    nSamplesTotal = 0
+    for f in features:
+        nSamplesTotal += f.shape[0]
+    if nSamplesTotal > 1000 and nExp > 50:
+        nExp = 50
+        print "Number of training experiments changed to 50 due to high number of samples"
+    if nSamplesTotal > 2000 and nExp > 10:
+        nExp = 10
+        print "Number of training experiments changed to 10 due to high number of samples"
+
+    for Ci, C in enumerate(Params):                # for each param value
+                CM = numpy.zeros((nClasses, nClasses))
+                for e in range(nExp):              # for each cross-validation iteration:
+                    print "Param = {0:.5f} - Classifier Evaluation Experiment {1:d} of {2:d}".format(C, e+1, nExp)
+                    # split features:
+                    featuresTrain, featuresTest = randSplitFeatures(featuresNorm, perTrain)
+                    # train multi-class svms:
+                    if ClassifierName == "svm":
+                        Classifier = trainSVM(featuresTrain, C)
+                    elif ClassifierName == "svm_rbf":
+                        Classifier = trainSVM_RBF(featuresTrain, C)
+                    elif ClassifierName == "knn":
+                        Classifier = trainKNN(featuresTrain, C)
+                    elif ClassifierName == "randomforest":
+                        Classifier = trainRandomForest(featuresTrain, C)
+                    elif ClassifierName == "gradientboosting":
+                        Classifier = trainGradientBoosting(featuresTrain, C)
+                    elif ClassifierName == "extratrees":
+                        Classifier = trainExtraTrees(featuresTrain, C)
+
+                    CMt = numpy.zeros((nClasses, nClasses))
+                    for c1 in range(nClasses):
+                        #Results = Classifier.pred(featuresTest[c1])
+                        nTestSamples = len(featuresTest[c1])
+                        Results = numpy.zeros((nTestSamples, 1))
+                        for ss in range(nTestSamples):
+                            [Results[ss], _] = classifierWrapper(Classifier, ClassifierName, featuresTest[c1][ss])
+                        for c2 in range(nClasses):
+                            CMt[c1][c2] = float(len(numpy.nonzero(Results == c2)[0]))
+                    CM = CM + CMt
+                CM = CM + 0.0000000010
+                Rec = numpy.zeros((CM.shape[0], ))
+                Pre = numpy.zeros((CM.shape[0], ))
+
+                for ci in range(CM.shape[0]):
+                    Rec[ci] = CM[ci, ci] / numpy.sum(CM[ci, :])
+                    Pre[ci] = CM[ci, ci] / numpy.sum(CM[:, ci])
+                PrecisionClassesAll.append(Pre)
+                RecallClassesAll.append(Rec)
+                F1 = 2 * Rec * Pre / (Rec + Pre)
+                F1ClassesAll.append(F1)
+                acAll.append(numpy.sum(numpy.diagonal(CM)) / numpy.sum(CM))
+
+                CMsAll.append(CM)
+                F1All.append(numpy.mean(F1))
+                # print "{0:6.4f}{1:6.4f}{2:6.1f}{3:6.1f}".format(nu, g, 100.0*acAll[-1], 100.0*F1All[-1])
+
+    print ("\t\t"),
+    for i, c in enumerate(ClassNames):
+        if i == len(ClassNames)-1:
+            print "{0:s}\t\t".format(c),
+        else:
+            print "{0:s}\t\t\t".format(c),
+    print ("OVERALL")
+    print ("\tC"),
+    for c in ClassNames:
+        print "\tPRE\tREC\tF1",
+    print "\t{0:s}\t{1:s}".format("ACC", "F1")
+    bestAcInd = numpy.argmax(acAll)
+    bestF1Ind = numpy.argmax(F1All)
+    for i in range(len(PrecisionClassesAll)):
+        print "\t{0:.3f}".format(Params[i]),
+        for c in range(len(PrecisionClassesAll[i])):
+            print "\t{0:.1f}\t{1:.1f}\t{2:.1f}".format(100.0 * PrecisionClassesAll[i][c], 100.0 * RecallClassesAll[i][c], 100.0 * F1ClassesAll[i][c]),
+        print "\t{0:.1f}\t{1:.1f}".format(100.0 * acAll[i], 100.0 * F1All[i]),
+        if i == bestF1Ind:
+            print "\t best F1",
+        if i == bestAcInd:
+            print "\t best Acc",
+        print
+
+    if parameterMode == 0:    # keep parameters that maximize overall classification accuracy:
+        print "Confusion Matrix:"
+        printConfusionMatrix(CMsAll[bestAcInd], ClassNames)
+        return Params[bestAcInd]
+    elif parameterMode == 1:  # keep parameters that maximize overall F1 measure:
+        print "Confusion Matrix:"
+        printConfusionMatrix(CMsAll[bestF1Ind], ClassNames)
+        return Params[bestF1Ind]
+
+
+def evaluateRegression(features, labels, nExp, MethodName, Params):
+    '''
+    ARGUMENTS:
+        features:     numpy matrices of features [numOfSamples x numOfDimensions]
+        labels:       list of sample labels
+        nExp:         number of cross-validation experiments
+        MethodName:   "svm" or "randomforest"
+        Params:       list of classifier params to be evaluated
+    RETURNS:
+         bestParam:   the value of the input parameter that optimizes the selected performance measure
+    '''
+
+    # feature normalization:
+    (featuresNorm, MEAN, STD) = normalizeFeatures([features])
+    featuresNorm = featuresNorm[0]
+    nSamples = labels.shape[0]
+    partTrain = 0.9
+    ErrorsAll = []
+    ErrorsTrainAll = []
+    ErrorsBaselineAll = []
+    for Ci, C in enumerate(Params):                # for each param value
+                Errors = []
+                ErrorsTrain = []
+                ErrorsBaseline = []
+                for e in range(nExp):             # for each cross-validation iteration:
+                    # split features:
+                    randperm = numpy.random.permutation(range(nSamples))
+                    nTrain = int(round(partTrain * nSamples))
+                    featuresTrain = [featuresNorm[randperm[i]] for i in range(nTrain)]
+                    featuresTest = [featuresNorm[randperm[i+nTrain]] for i in range(nSamples - nTrain)]
+                    labelsTrain = [labels[randperm[i]] for i in range(nTrain)]
+                    labelsTest = [labels[randperm[i + nTrain]] for i in range(nSamples - nTrain)]
+
+                    # train multi-class svms:                    
+                    featuresTrain = numpy.matrix(featuresTrain)                                 
+                    if MethodName == "svm":                                        
+                        [Classifier, trainError] = trainSVMregression(featuresTrain, labelsTrain, C)     
+                    elif MethodName == "svm_rbf":                      
+                        [Classifier, trainError] = trainSVMregression_rbf(featuresTrain, labelsTrain, C)                                             
+                    elif MethodName == "randomforest":
+                        [Classifier, trainError] = trainRandomForestRegression(featuresTrain, labelsTrain, C)
+                    ErrorTest = []
+                    ErrorTestBaseline = []
+                    for itest, fTest in enumerate(featuresTest):
+                        R = regressionWrapper(Classifier, MethodName, fTest)
+                        Rbaseline = numpy.mean(labelsTrain)
+                        ErrorTest.append((R - labelsTest[itest]) * (R - labelsTest[itest]))
+                        ErrorTestBaseline.append((Rbaseline - labelsTest[itest]) * (Rbaseline - labelsTest[itest]))
+                    Error = numpy.array(ErrorTest).mean()
+                    ErrorBaseline = numpy.array(ErrorTestBaseline).mean()
+                    Errors.append(Error)
+                    ErrorsTrain.append(trainError)
+                    ErrorsBaseline.append(ErrorBaseline)
+                ErrorsAll.append(numpy.array(Errors).mean())
+                ErrorsTrainAll.append(numpy.array(ErrorsTrain).mean())
+                ErrorsBaselineAll.append(numpy.array(ErrorsBaseline).mean())
+
+    bestInd = numpy.argmin(ErrorsAll)
+
+    print "{0:s}\t\t{1:s}\t\t{2:s}\t\t{3:s}".format("Param", "MSE", "T-MSE", "R-MSE")
+    for i in range(len(ErrorsAll)):
+        print "{0:.4f}\t\t{1:.2f}\t\t{2:.2f}\t\t{3:.2f}".format(Params[i], ErrorsAll[i], ErrorsTrainAll[i], ErrorsBaselineAll[i]),
+        if i == bestInd:
+            print "\t\t best",
+        print
+    return Params[bestInd], ErrorsAll[bestInd], ErrorsBaselineAll[bestInd]
+
+
+def printConfusionMatrix(CM, ClassNames):
+    '''
+    This function prints a confusion matrix for a particular classification task.
+    ARGUMENTS:
+        CM:            a 2-D numpy array of the confusion matrix
+                       (CM[i,j] is the number of times a sample from class i was classified in class j)
+        ClassNames:    a list that contains the names of the classes
+    '''
+
+    if CM.shape[0] != len(ClassNames):
+        print "printConfusionMatrix: Wrong argument sizes\n"
+        return
+
+    for c in ClassNames:
+        if len(c) > 4:
+            c = c[0:3]
+        print "\t{0:s}".format(c),
+    print
+
+    for i, c in enumerate(ClassNames):
+        if len(c) > 4:
+            c = c[0:3]
+        print "{0:s}".format(c),
+        for j in range(len(ClassNames)):
+            print "\t{0:.2f}".format(100.0 * CM[i][j] / numpy.sum(CM)),
+        print
+
+
+def normalizeFeatures(features):
+    '''
+    This function normalizes a feature set to 0-mean and 1-std.
+    Used in most classifier trainning cases.
+
+    ARGUMENTS:
+        - features:    list of feature matrices (each one of them is a numpy matrix)
+    RETURNS:
+        - featuresNorm:    list of NORMALIZED feature matrices
+        - MEAN:        mean vector
+        - STD:        std vector
+    '''
+    X = numpy.array([])
+
+    for count, f in enumerate(features):
+        if f.shape[0] > 0:
+            if count == 0:
+                X = f
+            else:
+                X = numpy.vstack((X, f))
+            count += 1
+
+    MEAN = numpy.mean(X, axis=0) + 0.00000000000001;
+    STD = numpy.std(X, axis=0) + 0.00000000000001;
+
+    featuresNorm = []
+    for f in features:
+        ft = f.copy()
+        for nSamples in range(f.shape[0]):
+            ft[nSamples, :] = (ft[nSamples, :] - MEAN) / STD
+        featuresNorm.append(ft)
+    return (featuresNorm, MEAN, STD)
+
+
+def listOfFeatures2Matrix(features):
+    '''
+    listOfFeatures2Matrix(features)
+
+    This function takes a list of feature matrices as argument and returns a single concatenated feature matrix and the respective class labels.
+
+    ARGUMENTS:
+        - features:        a list of feature matrices
+
+    RETURNS:
+        - X:            a concatenated matrix of features
+        - Y:            a vector of class indeces
+    '''
+
+    X = numpy.array([])
+    Y = numpy.array([])
+    for i, f in enumerate(features):
+        if i == 0:
+            X = f
+            Y = i * numpy.ones((len(f), 1))
+        else:
+            X = numpy.vstack((X, f))
+            Y = numpy.append(Y, i * numpy.ones((len(f), 1)))
+    return (X, Y)
+
+
+def pcaDimRed(features, nDims):
+    [X, Y] = listOfFeatures2Matrix(features)
+    pca = sklearn.decomposition.PCA(n_components = nDims)
+    pca.fit(X)
+    coeff = pca.components_
+    coeff = coeff[:, 0:nDims]
+
+    featuresNew = []
+    for f in features:
+        ft = f.copy()
+#        ft = pca.transform(ft, k=nDims)
+        ft = numpy.dot(f, coeff)
+        featuresNew.append(ft)
+
+    return (featuresNew, coeff)
+
+
+def fileClassification(inputFile, modelName, modelType):
+    # Load classifier:
+
+    if not os.path.isfile(modelName):
+        print "fileClassification: input modelName not found!"
+        return (-1, -1, -1)
+
+    if not os.path.isfile(inputFile):
+        print "fileClassification: wav file not found!"
+        return (-1, -1, -1)
+
+    if (modelType) == 'svm' or (modelType == 'svm_rbf'):
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = loadSVModel(modelName)
+    elif modelType == 'knn':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = loadKNNModel(modelName)
+    elif modelType == 'randomforest':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = loadRandomForestModel(modelName)
+    elif modelType == 'gradientboosting':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = loadGradientBoostingModel(modelName)
+    elif modelType == 'extratrees':
+        [Classifier, MEAN, STD, classNames, mtWin, mtStep, stWin, stStep, computeBEAT] = loadExtraTreesModel(modelName)
+
+    [Fs, x] = audioBasicIO.readAudioFile(inputFile)        # read audio file and convert to mono
+    x = audioBasicIO.stereo2mono(x)
+
+    if isinstance(x, int):                                 # audio file IO problem
+        return (-1, -1, -1)
+    if x.shape[0] / float(Fs) <= mtWin:
+        return (-1, -1, -1)
+
+    # feature extraction:
+    [MidTermFeatures, s] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs * stWin), round(Fs * stStep))
+    MidTermFeatures = MidTermFeatures.mean(axis=1)        # long term averaging of mid-term statistics
+    if computeBEAT:
+        [beat, beatConf] = aF.beatExtraction(s, stStep)
+        MidTermFeatures = numpy.append(MidTermFeatures, beat)
+        MidTermFeatures = numpy.append(MidTermFeatures, beatConf)
+    curFV = (MidTermFeatures - MEAN) / STD                # normalization
+
+    [Result, P] = classifierWrapper(Classifier, modelType, curFV)    # classification        
+    return Result, P, classNames
+
+
+def fileRegression(inputFile, modelName, modelType):
+    # Load classifier:
+
+    if not os.path.isfile(inputFile):
+        print "fileClassification: wav file not found!"
+        return (-1, -1, -1)
+
+    regressionModels = glob.glob(modelName + "_*")
+    regressionModels2 = []
+    for r in regressionModels:
+        if r[-5::] != "MEANS":
+            regressionModels2.append(r)
+    regressionModels = regressionModels2
+    regressionNames = []
+    for r in regressionModels:
+        regressionNames.append(r[r.rfind("_")+1::])
+
+    # FEATURE EXTRACTION
+    # LOAD ONLY THE FIRST MODEL (for mtWin, etc)
+    if modelType == 'svm' or modelType == "svm_rbf":        
+        [_, _, _, mtWin, mtStep, stWin, stStep, computeBEAT] = loadSVModel(regressionModels[0], True)
+    elif modelType == 'randomforest':
+        [_, _, _, mtWin, mtStep, stWin, stStep, computeBEAT] = loadRandomForestModel(regressionModels[0], True)
+
+    [Fs, x] = audioBasicIO.readAudioFile(inputFile)        # read audio file and convert to mono
+    x = audioBasicIO.stereo2mono(x)
+    # feature extraction:
+    [MidTermFeatures, s] = aF.mtFeatureExtraction(x, Fs, mtWin * Fs, mtStep * Fs, round(Fs * stWin), round(Fs * stStep))
+    MidTermFeatures = MidTermFeatures.mean(axis=1)        # long term averaging of mid-term statistics
+    if computeBEAT:
+        [beat, beatConf] = aF.beatExtraction(s, stStep)
+        MidTermFeatures = numpy.append(MidTermFeatures, beat)
+        MidTermFeatures = numpy.append(MidTermFeatures, beatConf)
+
+    # REGRESSION
+    R = []
+    for ir, r in enumerate(regressionModels):
+        if not os.path.isfile(r):
+            print "fileClassification: input modelName not found!"
+            return (-1, -1, -1)
+        if modelType == 'svm' or modelType == "svm_rbf":
+            [Model, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT] = loadSVModel(r, True)
+        elif modelType == 'randomforest':
+            [Model, MEAN, STD, mtWin, mtStep, stWin, stStep, computeBEAT] = loadRandomForestModel(r, True)
+        curFV = (MidTermFeatures - MEAN) / STD                  # normalization
+        R.append(regressionWrapper(Model, modelType, curFV))    # classification
+    return R, regressionNames
+
+
+def lda(data, labels, redDim):
+    # Centre data
+    data -= data.mean(axis=0)
+    nData = numpy.shape(data)[0]
+    nDim = numpy.shape(data)[1]
+    print nData, nDim
+    Sw = numpy.zeros((nDim, nDim))
+    Sb = numpy.zeros((nDim, nDim))
+
+    C = numpy.cov((data.T))
+
+    # Loop over classes
+    classes = numpy.unique(labels)
+    for i in range(len(classes)):
+        # Find relevant datapoints
+        indices = (numpy.where(labels == classes[i]))
+        d = numpy.squeeze(data[indices, :])
+        classcov = numpy.cov((d.T))
+        Sw += float(numpy.shape(indices)[0])/nData * classcov
+
+    Sb = C - Sw
+    # Now solve for W
+    # Compute eigenvalues, eigenvectors and sort into order
+    #evals,evecs = linalg.eig(dot(linalg.pinv(Sw),sqrt(Sb)))
+    evals, evecs = la.eig(Sw, Sb)
+    indices = numpy.argsort(evals)
+    indices = indices[::-1]
+    evecs = evecs[:, indices]
+    evals = evals[indices]
+    w = evecs[:, :redDim]
+    #print evals, w
+
+    newData = numpy.dot(data, w)
+    #for i in range(newData.shape[0]):
+    #    plt.text(newData[i,0],newData[i,1],str(labels[i]))
+
+    #plt.xlim([newData[:,0].min(), newData[:,0].max()])
+    #plt.ylim([newData[:,1].min(), newData[:,1].max()])
+    #plt.show()
+    return newData, w
+
+
+def writeTrainDataToARFF(modelName, features, classNames, featureNames):
+    f = open(modelName + ".arff", 'w')
+    f.write('@RELATION ' + modelName + '\n')
+    for fn in featureNames:
+        f.write('@ATTRIBUTE ' + fn + ' NUMERIC\n')
+    f.write('@ATTRIBUTE class {')
+    for c in range(len(classNames)-1):
+        f.write(classNames[c] + ',')
+    f.write(classNames[-1] + '}\n\n')
+    f.write('@DATA\n')
+    for c, fe in enumerate(features):
+        for i in range(fe.shape[0]):
+            for j in range(fe.shape[1]):
+                f.write("{0:f},".format(fe[i, j]))
+            f.write(classNames[c]+"\n")
+    f.close()
+
+
+def trainSpeakerModelsScript():
+    '''
+    This script is used to train the speaker-related models (NOTE: data paths are hard-coded and NOT included in the library, the models are, however included)
+         import audioTrainTest as aT
+        aT.trainSpeakerModelsScript()
+
+    '''
+    mtWin = 2.0
+    mtStep = 2.0
+    stWin = 0.020
+    stStep = 0.020
+
+    dirName = "DIARIZATION_ALL/all"
+    listOfDirs = [os.path.join(dirName, name) for name in os.listdir(dirName) if os.path.isdir(os.path.join(dirName, name))]
+    featureAndTrain(listOfDirs, mtWin, mtStep, stWin, stStep, "knn", "data/knnSpeakerAll", computeBEAT=False, perTrain=0.50)
+
+    dirName = "DIARIZATION_ALL/female_male"
+    listOfDirs = [os.path.join(dirName, name) for name in os.listdir(dirName) if os.path.isdir(os.path.join(dirName, name))]
+    featureAndTrain(listOfDirs, mtWin, mtStep, stWin, stStep, "knn", "data/knnSpeakerFemaleMale", computeBEAT=False, perTrain=0.50)
+
+
+def main(argv):
+    return 0
+
+if __name__ == '__main__':
+    main(sys.argv)
diff --git a/audioVisualization.py b/audioVisualization.py
new file mode 100644
index 0000000..0f9b477
--- /dev/null
+++ b/audioVisualization.py
@@ -0,0 +1,217 @@
+import os, sys, shutil, struct, simplejson
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.spatial import distance
+from pylab import *
+import ntpath
+import audioFeatureExtraction as aF    
+import audioTrainTest as aT
+import sklearn
+import sklearn.discriminant_analysis
+
+def generateColorMap():
+    '''
+    This function generates a 256 jet colormap of HTML-like hex string colors (e.g. FF88AA)
+    '''
+    Map = cm.jet(np.arange(256))
+    stringColors = []
+    for i in range(Map.shape[0]):
+        rgb = (int(255*Map[i][0]), int(255*Map[i][1]), int(255*Map[i][2]))
+        stringColors.append(struct.pack('BBB',*rgb).encode('hex'))
+    return stringColors
+
+def levenshtein(str1, s2):
+    '''
+    Distance between two strings
+    '''
+    N1 = len(str1)
+    N2 = len(s2)
+
+    stringRange = [range(N1 + 1)] * (N2 + 1)
+    for i in range(N2 + 1):
+        stringRange[i] = range(i,i + N1 + 1)
+    for i in range(0,N2):
+        for j in range(0,N1):
+            if str1[j] == s2[i]:
+                stringRange[i+1][j+1] = min(stringRange[i+1][j] + 1, stringRange[i][j+1] + 1, stringRange[i][j])
+            else:
+                stringRange[i+1][j+1] = min(stringRange[i+1][j] + 1, stringRange[i][j+1] + 1, stringRange[i][j] + 1)
+    return stringRange[N2][N1]
+
+def textListToColors(names):
+    '''
+    Generates a list of colors based on a list of names (strings). Similar strings correspond to similar colors.
+    '''
+    # STEP A: compute strings distance between all combnations of strings
+    Dnames = np.zeros( (len(names), len(names)) )
+    for i in range(len(names)):
+        for j in range(len(names)):
+            Dnames[i,j] = 1 - 2.0 * levenshtein(names[i], names[j]) / float(len(names[i]+names[j]))
+
+    # STEP B: pca dimanesionality reduction to a single-dimension (from the distance space)
+    pca = sklearn.decomposition.PCA(n_components = 1)
+    pca.fit(Dnames)    
+    
+    # STEP C: mapping of 1-dimensional values to colors in a jet-colormap
+    textToColor = pca.transform(Dnames)
+    textToColor = 255 * (textToColor - textToColor.min()) / (textToColor.max() - textToColor.min())
+    textmaps = generateColorMap();
+    colors = [textmaps[int(c)] for c in textToColor]
+    return colors
+
+def textListToColorsSimple(names):
+    '''
+    Generates a list of colors based on a list of names (strings). Similar strings correspond to similar colors. 
+    '''
+    uNames = list(set(names))
+    uNames.sort()
+    textToColor = [ uNames.index(n) for n in names ]
+    textToColor = np.array(textToColor)
+    textToColor = 255 * (textToColor - textToColor.min()) / (textToColor.max() - textToColor.min())
+    textmaps = generateColorMap();
+    colors = [textmaps[int(c)] for c in textToColor]
+    
+    # colors = [c for (n, c) in sorted(zip(names, colors))]
+
+    return colors
+
+def chordialDiagram(fileStr, SM, Threshold, names, namesCategories):
+    '''
+    Generates a d3js chordial diagram that illustrates similarites
+    '''
+    #colors = textListToColors(namesCategories)
+    colors = textListToColorsSimple(namesCategories)
+    SM2 = SM.copy()
+    SM2 = (SM2 + SM2.T) / 2.0
+    for i in range(SM2.shape[0]):
+        M = Threshold
+#        a = np.sort(SM2[i,:])[::-1]
+#        M = np.mean(a[0:int(SM2.shape[1]/3+1)])
+        SM2[i,SM2[i,:]<M] = 0;
+
+    dirChordial = fileStr + "_Chordial"
+    if not os.path.isdir(dirChordial):
+        os.mkdir(dirChordial)
+    jsonPath         = dirChordial + os.sep + "matrix.json"
+    namesPath        = dirChordial + os.sep + "Names.csv"
+ 
+    jsonSMMatrix = simplejson.dumps(SM2.tolist())
+    f = open(jsonPath,'w'); f.write(jsonSMMatrix);  f.close()
+    f = open(namesPath,'w'); f.write("name,color\n"); 
+    for i, n in enumerate(names):
+        f.write("{0:s},{1:s}\n".format(n,"#"+colors[i]))
+    f.close()
+
+    shutil.copyfile("data/similarities.html", dirChordial+os.sep+"similarities.html")
+    shutil.copyfile("data/style.css", dirChordial+os.sep+"style.css")
+
+def visualizeFeaturesFolder(folder, dimReductionMethod, priorKnowledge = "none"):
+    '''
+    This function generates a chordial visualization for the recordings of the provided path.
+    ARGUMENTS:
+        - folder:        path of the folder that contains the WAV files to be processed
+        - dimReductionMethod:    method used to reduce the dimension of the initial feature space before computing the similarity.
+        - priorKnowledge:    if this is set equal to "artist"
+    '''
+    if dimReductionMethod=="pca":
+        allMtFeatures, wavFilesList = aF.dirWavFeatureExtraction(folder, 30.0, 30.0, 0.050, 0.050, computeBEAT = True)
+        if allMtFeatures.shape[0]==0:
+            print "Error: No data found! Check input folder"
+            return
+        
+        namesCategoryToVisualize = [ntpath.basename(w).replace('.wav','').split(" --- ")[0] for w in wavFilesList]; 
+        namesToVisualize       = [ntpath.basename(w).replace('.wav','') for w in wavFilesList]; 
+
+        (F, MEAN, STD) = aT.normalizeFeatures([allMtFeatures])
+        F = np.concatenate(F)
+        
+        # check that the new PCA dimension is at most equal to the number of samples
+        K1 = 2
+        K2 = 10
+        if K1 > F.shape[0]:
+            K1 = F.shape[0]
+        if K2 > F.shape[0]:
+            K2 = F.shape[0]
+        pca1 = sklearn.decomposition.PCA(n_components = K1)
+        pca1.fit(F)        
+        pca2 = sklearn.decomposition.PCA(n_components = K2)
+        pca2.fit(F)        
+
+        finalDims = pca1.transform(F)
+        finalDims2 = pca2.transform(F)
+    else:    
+        allMtFeatures, Ys, wavFilesList = aF.dirWavFeatureExtractionNoAveraging(folder, 20.0, 5.0, 0.040, 0.040) # long-term statistics cannot be applied in this context (LDA needs mid-term features)
+        if allMtFeatures.shape[0]==0:
+            print "Error: No data found! Check input folder"
+            return
+        
+        namesCategoryToVisualize = [ntpath.basename(w).replace('.wav','').split(" --- ")[0] for w in wavFilesList]; 
+        namesToVisualize       = [ntpath.basename(w).replace('.wav','') for w in wavFilesList]; 
+
+        ldaLabels = Ys
+        if priorKnowledge=="artist":
+            uNamesCategoryToVisualize = list(set(namesCategoryToVisualize))
+            YsNew = np.zeros( Ys.shape )
+            for i, uname in enumerate(uNamesCategoryToVisualize):        # for each unique artist name:
+                indicesUCategories = [j for j, x in enumerate(namesCategoryToVisualize) if x == uname]
+                for j in indicesUCategories:
+                    indices = np.nonzero(Ys==j)
+                    YsNew[indices] = i
+            ldaLabels = YsNew
+
+        (F, MEAN, STD) = aT.normalizeFeatures([allMtFeatures])
+        F = np.array(F[0])
+
+        clf = sklearn.discriminant_analysis.LinearDiscriminantAnalysis(n_components=10)
+        clf.fit(F, ldaLabels)    
+        reducedDims =  clf.transform(F)
+
+        pca = sklearn.decomposition.PCA(n_components = 2)
+        pca.fit(reducedDims)
+        reducedDims = pca.transform(reducedDims)
+
+        # TODO: CHECK THIS ... SHOULD LDA USED IN SEMI-SUPERVISED ONLY????
+
+        uLabels = np.sort(np.unique((Ys)))        # uLabels must have as many labels as the number of wavFilesList elements
+        reducedDimsAvg = np.zeros( (uLabels.shape[0], reducedDims.shape[1] ) )
+        finalDims = np.zeros( (uLabels.shape[0], 2) ) 
+        for i, u in enumerate(uLabels):
+            indices = [j for j, x in enumerate(Ys) if x == u]
+            f = reducedDims[indices, :]
+            finalDims[i, :] = f.mean(axis=0)
+        finalDims2 = reducedDims
+
+    for i in range(finalDims.shape[0]):            
+        plt.text(finalDims[i,0], finalDims[i,1], ntpath.basename(wavFilesList[i].replace('.wav','')), horizontalalignment='center', verticalalignment='center', fontsize=10)
+        plt.plot(finalDims[i,0], finalDims[i,1], '*r')
+    plt.xlim([1.2*finalDims[:,0].min(), 1.2*finalDims[:,0].max()])
+    plt.ylim([1.2*finalDims[:,1].min(), 1.2*finalDims[:,1].max()])            
+    plt.show()
+
+    SM = 1.0 - distance.squareform(distance.pdist(finalDims2, 'cosine'))
+    for i in range(SM.shape[0]):
+        SM[i,i] = 0.0;
+
+
+    chordialDiagram("visualization", SM, 0.50, namesToVisualize, namesCategoryToVisualize)
+
+    SM = 1.0 - distance.squareform(distance.pdist(F, 'cosine'))
+    for i in range(SM.shape[0]):
+        SM[i,i] = 0.0;
+    chordialDiagram("visualizationInitial", SM, 0.50, namesToVisualize, namesCategoryToVisualize)
+
+    # plot super-categories (i.e. artistname
+    uNamesCategoryToVisualize = sort(list(set(namesCategoryToVisualize)))
+    finalDimsGroup = np.zeros( (len(uNamesCategoryToVisualize), finalDims2.shape[1] ) )
+    for i, uname in enumerate(uNamesCategoryToVisualize):
+        indices = [j for j, x in enumerate(namesCategoryToVisualize) if x == uname]
+        f = finalDims2[indices, :]
+        finalDimsGroup[i, :] = f.mean(axis=0)
+
+    SMgroup = 1.0 - distance.squareform(distance.pdist(finalDimsGroup, 'cosine'))
+    for i in range(SMgroup.shape[0]):
+        SMgroup[i,i] = 0.0;
+    chordialDiagram("visualizationGroup", SMgroup, 0.50, uNamesCategoryToVisualize, uNamesCategoryToVisualize)
+
+
+
diff --git a/convertToWav.py b/convertToWav.py
new file mode 100644
index 0000000..04f84d4
--- /dev/null
+++ b/convertToWav.py
@@ -0,0 +1,25 @@
+# USAGE:
+# convertToWav <folder path> <sampling rate> <number of channels>
+#
+
+import glob, sys, os
+
+def getVideoFilesFromFolder(dirPath):
+	types = (dirPath+os.sep+'*.avi', dirPath+os.sep+'*.mkv', dirPath+os.sep+'*.mp4', dirPath+os.sep+'*.mp3', dirPath+os.sep+'*.flac', dirPath+os.sep+'*.ogg') # the tuple of file types
+	files_grabbed = []
+	for files in types:
+		files_grabbed.extend(glob.glob(files))
+	return files_grabbed
+
+def main(argv):
+	if (len(argv)==4):
+		files = getVideoFilesFromFolder(argv[1])
+		samplingRate = int(argv[2])
+		channels = int(argv[3])
+	
+		for f in files:
+			ffmpegString = 'avconv -i ' + '\"' + f + '\"' + ' -ar ' + str(samplingRate) + ' -ac ' + str(channels) + ' ' + '\"' + os.path.splitext(f)[0] + '\"' + '.wav' 
+			os.system(ffmpegString)
+
+if __name__ == '__main__':
+	main(sys.argv)
diff --git a/feature_extraction.py b/feature_extraction.py
new file mode 100644
index 0000000..13951f4
--- /dev/null
+++ b/feature_extraction.py
@@ -0,0 +1,23 @@
+import audioFeatureExtraction
+import audioBasicIO
+import numpy as np
+import matplotlib.pyplot as plt
+
+def tutorial():
+    [Fs, x] = audioBasicIO.readAudioFile("/home/parallels/git_repos/pyAudioAnalysis/data/scottish.wav")
+    F = audioFeatureExtraction.stFeatureExtraction(x, Fs, 0.050 * Fs, 0.025 * Fs)
+    plt.subplot(2, 1, 1)
+    plt.plot(F[0, :])
+    plt.xlabel('Frame no')
+    plt.ylabel('ZCR')
+    plt.subplot(2, 1, 2)
+    plt.plot(F[1, :])
+    plt.xlabel('Frame no')
+    plt.ylabel('Energy')
+    plt.savefig('plot1.png', format='png')
+
+
+
+if __name__ == "__main__":
+    np.load("/home/parallels/git_repos/pyAudioAnalysis/data/speech_music_sample.wav.npy")
+    p
\ No newline at end of file
diff --git a/plot1.png b/plot1.png
new file mode 100644
index 0000000..38d510a
Binary files /dev/null and b/plot1.png differ
diff --git a/utilities.py b/utilities.py
new file mode 100644
index 0000000..7fa8611
--- /dev/null
+++ b/utilities.py
@@ -0,0 +1,103 @@
+import sys, os, numpy
+
+def isfloat(x):
+	"""
+	Check if argument is float
+	"""
+	try:
+		a = float(x)
+	except ValueError:
+		return False
+	else:
+		return True
+
+def isint(x):
+	"""
+	Check if argument is int
+	"""
+	try:
+		a = float(x)
+		b = int(a)
+	except ValueError:
+		return False
+	else:
+		return a == b
+
+def isNum(x):
+	"""
+	Check if string argument is numerical
+	"""
+	return isfloat(x) or isint(x)
+
+
+def peakdet(v, delta, x = None):
+    """
+    Converted from MATLAB script at http://billauer.co.il/peakdet.html
+    
+    Returns two arrays
+    
+    function [maxtab, mintab]=peakdet(v, delta, x)
+    %PEAKDET Detect peaks in a vector
+    %        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
+    %        maxima and minima ("peaks") in the vector V.
+    %        MAXTAB and MINTAB consists of two columns. Column 1
+    %        contains indices in V, and column 2 the found values.
+    %      
+    %        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
+    %        in MAXTAB and MINTAB are replaced with the corresponding
+    %        X-values.
+    %
+    %        A point is considered a maximum peak if it has the maximal
+    %        value, and was preceded (to the left) by a value lower by
+    %        DELTA.
+    
+    % Eli Billauer, 3.4.05
+    % This function is released to the public domain; Any use is allowed.
+    
+    """
+    maxtab = []
+    mintab = []
+       
+    if x is None:
+        x = numpy.arange(len(v))
+    
+    v = numpy.asarray(v)
+    
+    if len(v) != len(x):
+        sys.exit('Input vectors v and x must have same length')
+    
+    if not numpy.isscalar(delta):
+        sys.exit('Input argument delta must be a scalar')
+    
+    if delta <= 0:
+        sys.exit('Input argument delta must be positive')
+    
+    mn, mx = numpy.Inf, -numpy.Inf
+    mnpos, mxpos = numpy.NaN, numpy.NaN
+    
+    lookformax = True
+    
+    for i in numpy.arange(len(v)):
+        this = v[i]
+        if this > mx:
+            mx = this
+            mxpos = x[i]
+        if this < mn:
+            mn = this
+            mnpos = x[i]
+        
+        if lookformax:
+            if this < mx-delta:
+                maxtab.append(mxpos)
+                mn = this
+                mnpos = x[i]
+                lookformax = False
+        else:
+            if this > mn+delta:
+                mintab.append(mnpos)
+                mx = this
+                mxpos = x[i]
+                lookformax = True
+ 
+    return numpy.array(maxtab), numpy.array(mintab)
+