# --- # jupyter: # jupytext: # text_representation: # extension: .py # format_name: percent # format_version: '1.3' # jupytext_version: 1.16.1 # kernelspec: # display_name: Python 3 (ipykernel) # language: python # name: python3 # --- # %% [markdown] # ## pipeline for all the files # %% import os import numpy as np import librosa from sklearn.preprocessing import StandardScaler import joblib import numpy as np from sklearn.cluster import KMeans from sklearn.decomposition import PCA import matplotlib.pyplot as plt import librosa from IPython.display import Audio, display from sklearn.model_selection import cross_val_score from sklearn.ensemble import RandomForestClassifier # %% ## cockpit for directories # Directory containing the audio files # audio_dir = "data/soundscape_data" audio_dir = ( "../data/SoundMeters_Ingles_Primary-20240519T132658Z-002/SoundMeters_Ingles_Primary" ) # Directory to save features features_dir = "../data/features" os.makedirs(features_dir, exist_ok=True) # Directory to save clusters information clusters_dir = "../data/clusters" os.makedirs(clusters_dir, exist_ok=True) # %% # Parameters for windowing window_size = 10 # window size in seconds hop_size = 10 # hop size in seconds # Define frequency bands (in Hz) bands = { "Sub-bass": (20, 60), "Bass": (60, 250), "Low Midrange": (250, 500), "Midrange": (500, 2000), "Upper Midrange": (2000, 4000), "Presence": (4000, 6000), "Brilliance": (6000, 20000), } # Iterate over each audio file in the directory for filename in os.listdir(audio_dir): if filename.endswith(".wav"): file_path = os.path.join(audio_dir, filename) y, sr = librosa.load(file_path, sr=None) # Convert window and hop size to samples window_samples = int(window_size * sr) hop_samples = int(hop_size * sr) # Total number of windows in the current file num_windows = (len(y) - window_samples) // hop_samples + 1 all_features = [] for i in range(num_windows): start_sample = i * hop_samples end_sample = start_sample + window_samples y_window = y[start_sample:end_sample] # Compute STFT S = librosa.stft(y_window) S_db = librosa.amplitude_to_db(np.abs(S)) # Compute features for each band features = [] for band, (low_freq, high_freq) in bands.items(): low_bin = int(np.floor(low_freq * (S.shape[0] / sr))) high_bin = int(np.ceil(high_freq * (S.shape[0] / sr))) band_energy = np.mean(S_db[low_bin:high_bin, :], axis=0) features.append(band_energy) # Flatten the feature array and add to all_features features_flat = np.concatenate(features) all_features.append(features_flat) # Convert to numpy array all_features = np.array(all_features) # Standardize features scaler = StandardScaler() all_features = scaler.fit_transform(all_features) # Save features to disk feature_file = os.path.join( features_dir, f"{os.path.splitext(filename)[0]}_features.npy" ) joblib.dump((all_features, scaler), feature_file) # %% # Number of clusters n_clusters = 5 # Load all features all_features = [] for feature_file in os.listdir(features_dir): if feature_file.endswith("_features.npy"): features, _ = joblib.load(os.path.join(features_dir, feature_file)) all_features.append(features) # Combine all features into a single array all_features = np.vstack(all_features) # Perform PCA for 2D visualization pca = PCA(n_components=2) features_pca = pca.fit_transform(all_features) # Perform k-means clustering kmeans = KMeans(n_clusters=n_clusters) # Example: 5 clusters clusters = kmeans.fit_predict(all_features) # Save clustering results clustering_results = {"clusters": clusters, "kmeans": kmeans, "pca": pca} joblib.dump(clustering_results, os.path.join(clusters_dir, "clustering_results.pkl")) # %% # Plot the PCA-reduced features with cluster labels using a legend plt.figure(figsize=(10, 6)) # Define a colormap colors = plt.cm.tab10(np.arange(kmeans.n_clusters)) for cluster_label in np.unique(clusters): cluster_points = features_pca[clusters == cluster_label] plt.scatter( cluster_points[:, 0], cluster_points[:, 1], s=50, color=colors[cluster_label], label=f"Cluster {cluster_label}", ) plt.title("PCA of Clustered Frequency Band Features") plt.xlabel("Principal Component 1") plt.ylabel("Principal Component 2") plt.legend() plt.show() # %% # Plot the clusters plt.figure(figsize=(10, 6)) for i in range(n_clusters): plt.plot(all_features[clusters == i].mean(axis=0), label=f"Cluster {i}") plt.legend() plt.title("Clustered Frequency Band Features") plt.xlabel("Feature Index (Frequency Bands)") plt.ylabel("Mean Feature Value (Energy in dB)") plt.show() # %% # Function to plot the spectrogram def plot_spectrogram(y, sr, title): S = librosa.stft(y) S_db = librosa.amplitude_to_db(np.abs(S), ref=np.max) plt.figure(figsize=(10, 6)) librosa.display.specshow(S_db, sr=sr, x_axis="time", y_axis="log") plt.colorbar(format="%+2.0f dB") plt.title(title) plt.show() # %% # Parameters for windowing window_size = 10 # window size in seconds hop_size = 10 # hop size in seconds # Load clustering results clustering_results = joblib.load(os.path.join(clusters_dir, "clustering_results.pkl")) clusters = clustering_results["clusters"] # Load all features all_features = [] audio_segments = [] for feature_file in os.listdir(features_dir): if feature_file.endswith("_features.npy"): features, scaler = joblib.load(os.path.join(features_dir, feature_file)) filename = feature_file.replace("_features.npy", ".wav") file_path = os.path.join(audio_dir, filename) y, sr = librosa.load(file_path, sr=None) # Convert window and hop size to samples window_samples = int(window_size * sr) hop_samples = int(hop_size * sr) num_windows = (len(y) - window_samples) // hop_samples + 1 for i in range(num_windows): start_sample = i * hop_samples end_sample = start_sample + window_samples y_window = y[start_sample:end_sample] audio_segments.append(y_window) all_features.append(features) # Flatten the list of all features all_features = np.vstack(all_features) # %% for cluster_label in np.unique(clusters): try: # Find the first data point in the cluster representative_index = np.where(clusters == cluster_label)[0][0] # Use the original audio segment at the representative index y_representative = audio_segments[representative_index] # Check if y_representative is not empty if y_representative.size == 0: raise ValueError("The audio segment is empty") print(f"Cluster {cluster_label} representative audio:") display(Audio(data=y_representative, rate=sr)) # Plot the spectrogram plot_spectrogram( y_representative, sr, f"Spectrogram for Cluster {cluster_label}" ) except Exception as e: print(f"Could not play audio for cluster {cluster_label}: {e}") # %% scaler = StandardScaler() all_features_scaled = scaler.fit_transform(all_features) # Fit PCA pca = PCA().fit(all_features_scaled) # Method 1: Variance Explained explained_variance = pca.explained_variance_ratio_ cumulative_explained_variance = np.cumsum(explained_variance) # Plot the cumulative explained variance plt.figure(figsize=(10, 6)) plt.plot(cumulative_explained_variance, marker="o") plt.xlabel("Number of Principal Components") plt.ylabel("Cumulative Explained Variance") plt.title("Explained Variance vs. Number of Principal Components") plt.grid(True) plt.show() # %% # Method 2: Scree Plot plt.figure(figsize=(10, 6)) plt.plot(np.arange(1, len(explained_variance) + 1), explained_variance, marker="o") plt.xlabel("Principal Component") plt.ylabel("Explained Variance") plt.title("Scree Plot") plt.grid(True) plt.show() # %% # Method 3: Kaiser Criterion eigenvalues = pca.explained_variance_ kaiser_criterion = np.sum(eigenvalues > 1) # IMO this doesnt make sense at the moment, we need to extract more features print(f"Number of components selected by Kaiser Criterion: {kaiser_criterion}") # %% # %%