diff --git a/.gitignore b/.gitignore
index ee3d72b9ce85ed0673faf083f70c44b1fef0a0b7..b84d884c191ce7ccec3b5d428e3383921b66e0ae 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ Stats/.Rhistory
 Stats/.Rhistory
 Stats/.Rhistory
 Stats/.Rhistory
+Whistles/__pycache__/WhistleUtils.cpython-39.pyc
diff --git a/Whistles/1-Identification_of_whistles.py b/Whistles/1-Identification_of_whistles.py
index 906ed0f605c6cff15c5025b70b7ea0dd93974904..0db36477438f939d62e9118ce863575c5afb6470 100644
--- a/Whistles/1-Identification_of_whistles.py
+++ b/Whistles/1-Identification_of_whistles.py
@@ -22,7 +22,7 @@ print("Importation of packages complete!")
 #%% Parameters
 # Paths
 print("\rSetting up parameters...", end="\r")
-audio_f = "./../Audio_data"  # Path to recordings  
+audio_f = "./../Audio_data"  # Path to recordings
 csv_f = "./../CSV_data"      # Path to data in csv
 save_f = "./Trajectories"	 # Path were results are stored
 
@@ -99,4 +99,33 @@ for file in range(len(audio_paths)):
 	f.close()
 
 	print(f"\r\tFile {file+1} on {len(audio_paths)}: found {len(values)} whistles", end='\r')
-print("\nDetection of whistles finished!")
\ No newline at end of file
+print("\nDetection of whistles finished!")
+
+#%% Uncomment this section to plot the last audio  in loop
+# from WhistleUtils import plot_spectrums
+# final_traj[final_traj != 0] = 1 # binarisation
+
+# # generate bright colors to differenciate trajectories
+# prism = cm.get_cmap('prism', 256)
+# newcolors = prism(np.linspace(0, 1, np.unique(harmonized_traj).shape[0]))
+# pink = np.array([0/256, 0/256, 0/256, 1])
+# newcolors[0, :] = pink
+# newcmp = ListedColormap(newcolors)
+
+# # By default shows the whole audio
+# start = 0
+# stop = spectrum.shape[1]
+
+# # Create figure
+# fig, axs = plot_spectrums([amplitude_to_db(spectrum), 
+#                     max_loc_per_bin_check1, 
+#                     final_traj, 
+#                     harmonized_traj], 
+#                ['gray_r', 'gray', 'gray', newcmp], 
+#                titles=['Spectrogram (dB scale)', 'Local maxima selection', 'Extraction of continuous trajectories',
+#                'Exclusion of harmonics'], 
+#                ylabels=["Frequency"]*4,
+#                bins=375, title="")
+# # Update view
+# axs[3].set_xlim(start,stop)
+# plt.show(block=False)
\ No newline at end of file
diff --git a/Whistles/WhistleUtils.py b/Whistles/WhistleUtils.py
index 56b8271a34271d9fe9b223ae8c4ef017c5d893ba..487497485a99817462773eb8071c909e0677bfe1 100755
--- a/Whistles/WhistleUtils.py
+++ b/Whistles/WhistleUtils.py
@@ -17,6 +17,8 @@ from scipy.stats.mstats import gmean
 from scipy.spatial import distance
 from scipy import interpolate
 
+import matplotlib.pyplot as plt
+
 #%% Functions process
 def import_csv(csv_name, folder, useless=True, slash="/"):
     """
@@ -137,7 +139,7 @@ def get_local_maxima(spectrum, spectrum2, hardness, threshold=10e-5):
 
     for each_bin in range(spectrum.shape[1]):
         for freq in range(1, spectrum.shape[0]-1):
-            if (spectrum1[freq, each_bin] > spectrum2[freq-1, each_bin]) and \
+            if (spectrum[freq, each_bin] > spectrum2[freq-1, each_bin]) and \
             (spectrum2[freq, each_bin] > spectrum2[freq+1, each_bin]):
                 local_max1[freq, each_bin] = 1
                 if (spectrum[freq, each_bin] > (threshold)):
@@ -385,8 +387,8 @@ def get_category(samples, names, data_frame, category='acoustic'):
         
     Returns
     -------
-    TYPE
-        DESCRIPTION.
+    cat_list : LIST OF STRINGS
+        A list of categories corresponding to inputs.
 
     """
     use_samples = np.array([file[:23] for file in samples])
@@ -425,4 +427,60 @@ def get_category(samples, names, data_frame, category='acoustic'):
                       == 'AV')[0]
         cat_list[NC] = 'NC'
     
-    return cat_list
\ No newline at end of file
+    return cat_list
+
+
+def plot_spectrums(list_of_spectrums, cmaps=[], direction=-1, title="Whistles", 
+    bins=1, titles=[], ylabels=[]):
+    """
+    Parameters
+    ----------
+    list_of_spectrums : LIST OF NUMPY ARRAYS
+        List of images that will be plotted, should be a spectrogram
+        or a binary image (with 0s and 1s). 
+        All arrays should be of the same shape.
+    cmaps : LIST
+        List that must contain the same number of elements as list_of_spectrums.
+        Each element is either the name of a cmap (e.g. "viridis") or a custom cmap.
+    direction : INT, optionnal
+        Must be 1 or -1. If -1 it inverses the order for frequency display.
+        Default is -1.
+    title : STRING, optionnal
+        Title that will be displayed on top of the plot.
+    bins : FLOAT, optionnal
+        Number of bins (windows in spectrogram) per second
+    titles : LIST OF STRINGS
+        List that must contain the same number of elements as list_of_spectrums.
+        Titles associated with each spectrogram in plot.
+    ylabels : LIST OF STRINGS
+        List that must contain the same number of elements as list_of_spectrums.
+        Name of y-axis for each spectrogram.
+
+    Returns
+    -------
+    fig, axs : MATPLOTLIB OBJECT
+        figure and axis object that contain the plot of all given spectrogram.
+
+    """
+    n = len(list_of_spectrums)
+    fig, axs = plt.subplots(nrows=n, sharex=True, sharey=True, figsize=(15,15))
+
+    if len(cmaps) == 0:
+        cmaps = ['gray_r']*n
+
+    for i in range(n):
+        axs[i].imshow(list_of_spectrums[i][::direction], aspect='auto', 
+            interpolation='nearest', cmap=cmaps[i])
+        axs[i].set_yticklabels([])
+        if len(titles) > 0:
+            axs[i].set_title(titles[i], fontsize=30)
+        if len(ylabels) >0:
+            axs[i].set_ylabel(ylabels[i], fontsize=25)
+        axs[i].tick_params(axis='both', which='both', labelsize=20)
+    if bins !=1 :
+        axs[n-1].set_xlabel(f"Time in bins (1 bin = 1/{bins} sec)", fontsize=25)
+    else:
+        axs[n-1].set_xlabel(f"Time in bins)", fontsize=25)
+    fig.suptitle(title)
+    fig.tight_layout(pad=1)
+    return fig, axs
\ No newline at end of file