Source code for libfmp.c6.c6s1_onset_detection

"""
Module: libfmp.c6.c6s1_onset_detection
Author: Meinard Müller, Angel Villar-Corrales
License: The MIT license, https://opensource.org/licenses/MIT

This file is part of the FMP Notebooks (https://www.audiolabs-erlangen.de/FMP)
"""

import numpy as np
from numba import jit
from scipy import signal
from scipy.interpolate import interp1d
from scipy import ndimage
import librosa
import libfmp.b


[docs]def read_annotation_pos(fn_ann, label='', header=True, print_table=False): """Read and convert file containing either list of pairs (number,label) or list of (number) Notebook: C6/C6S1_OnsetDetection.ipynb Args: fn_ann (str): Name of file label (str): Name of label (Default value = '') header (bool): Assumes header (True) or not (False) (Default value = True) print_table (bool): Prints table if True (Default value = False) Returns: ann (list): List of annotations label_keys (dict): Dictionaries specifying color and line style used for labels """ df = libfmp.b.read_csv(fn_ann, header=header) if print_table: print(df) num_col = df.values[0].shape[0] if num_col == 1: df = df.assign(label=[label] * len(df.index)) ann = df.values.tolist() label_keys = {'beat': {'linewidth': 2, 'linestyle': ':', 'color': 'r'}, 'onset': {'linewidth': 1, 'linestyle': ':', 'color': 'r'}} return ann, label_keys
[docs]def compute_novelty_energy(x, Fs=1, N=2048, H=128, gamma=10.0, norm=True): """Compute energy-based novelty function Notebook: C6/C6S1_NoveltyEnergy.ipynb Args: x (np.ndarray): Signal Fs (scalar): Sampling rate (Default value = 1) N (int): Window size (Default value = 2048) H (int): Hop size (Default value = 128) gamma (float): Parameter for logarithmic compression (Default value = 10.0) norm (bool): Apply max norm (if norm==True) (Default value = True) Returns: novelty_energy (np.ndarray): Energy-based novelty function Fs_feature (scalar): Feature rate """ # x_power = x**2 w = signal.hann(N) Fs_feature = Fs / H energy_local = np.convolve(x**2, w**2, 'same') energy_local = energy_local[::H] if gamma is not None: energy_local = np.log(1 + gamma * energy_local) energy_local_diff = np.diff(energy_local) energy_local_diff = np.concatenate((energy_local_diff, np.array([0]))) novelty_energy = np.copy(energy_local_diff) novelty_energy[energy_local_diff < 0] = 0 if norm: max_value = max(novelty_energy) if max_value > 0: novelty_energy = novelty_energy / max_value return novelty_energy, Fs_feature
[docs]@jit(nopython=True) def compute_local_average(x, M): """Compute local average of signal Notebook: C6/C6S1_NoveltySpectral.ipynb Args: x (np.ndarray): Signal M (int): Determines size (2M+1) in samples of centric window used for local average Returns: local_average (np.ndarray): Local average signal """ L = len(x) local_average = np.zeros(L) for m in range(L): a = max(m - M, 0) b = min(m + M + 1, L) local_average[m] = (1 / (2 * M + 1)) * np.sum(x[a:b]) return local_average
[docs]def compute_novelty_spectrum(x, Fs=1, N=1024, H=256, gamma=100.0, M=10, norm=True): """Compute spectral-based novelty function Notebook: C6/C6S1_NoveltySpectral.ipynb Args: x (np.ndarray): Signal Fs (scalar): Sampling rate (Default value = 1) N (int): Window size (Default value = 1024) H (int): Hop size (Default value = 256) gamma (float): Parameter for logarithmic compression (Default value = 100.0) M (int): Size (frames) of local average (Default value = 10) norm (bool): Apply max norm (if norm==True) (Default value = True) Returns: novelty_spectrum (np.ndarray): Energy-based novelty function Fs_feature (scalar): Feature rate """ X = librosa.stft(x, n_fft=N, hop_length=H, win_length=N, window='hann') Fs_feature = Fs / H Y = np.log(1 + gamma * np.abs(X)) Y_diff = np.diff(Y) Y_diff[Y_diff < 0] = 0 novelty_spectrum = np.sum(Y_diff, axis=0) novelty_spectrum = np.concatenate((novelty_spectrum, np.array([0.0]))) if M > 0: local_average = compute_local_average(novelty_spectrum, M) novelty_spectrum = novelty_spectrum - local_average novelty_spectrum[novelty_spectrum < 0] = 0.0 if norm: max_value = max(novelty_spectrum) if max_value > 0: novelty_spectrum = novelty_spectrum / max_value return novelty_spectrum, Fs_feature
[docs]def principal_argument(v): """Principal argument function | Notebook: C6/C6S1_NoveltyPhase.ipynb, see also | Notebook: C8/C8S2_InstantFreqEstimation.ipynb Args: v (float or np.ndarray): Value (or vector of values) Returns: w (float or np.ndarray): Principle value of v """ w = np.mod(v + 0.5, 1) - 0.5 return w
[docs]def compute_novelty_phase(x, Fs=1, N=1024, H=64, M=40, norm=True): """Compute phase-based novelty function Notebook: C6/C6S1_NoveltyPhase.ipynb Args: x (np.ndarray): Signal Fs (scalar): Sampling rate (Default value = 1) N (int): Window size (Default value = 1024) H (int): Hop size (Default value = 64) M (int): Determines size (2M+1) in samples of centric window used for local average (Default value = 40) norm (bool): Apply max norm (if norm==True) (Default value = True) Returns: novelty_phase (np.ndarray): Energy-based novelty function Fs_feature (scalar): Feature rate """ X = librosa.stft(x, n_fft=N, hop_length=H, win_length=N, window='hann') Fs_feature = Fs / H phase = np.angle(X) / (2*np.pi) phase_diff = principal_argument(np.diff(phase, axis=1)) phase_diff2 = principal_argument(np.diff(phase_diff, axis=1)) novelty_phase = np.sum(np.abs(phase_diff2), axis=0) novelty_phase = np.concatenate((novelty_phase, np.array([0, 0]))) if M > 0: local_average = compute_local_average(novelty_phase, M) novelty_phase = novelty_phase - local_average novelty_phase[novelty_phase < 0] = 0 if norm: max_value = np.max(novelty_phase) if max_value > 0: novelty_phase = novelty_phase / max_value return novelty_phase, Fs_feature
[docs]def compute_novelty_complex(x, Fs=1, N=1024, H=64, gamma=10.0, M=40, norm=True): """Compute complex-domain novelty function Notebook: C6/C6S1_NoveltyComplex.ipynb Args: x (np.ndarray): Signal Fs (scalar): Sampling rate (Default value = 1) N (int): Window size (Default value = 1024) H (int): Hop size (Default value = 64) gamma (float): Parameter for logarithmic compression (Default value = 10.0) M (int): Determines size (2M+1) in samples of centric window used for local average (Default value = 40) norm (bool): Apply max norm (if norm==True) (Default value = True) Returns: novelty_complex (np.ndarray): Energy-based novelty function Fs_feature (scalar): Feature rate """ X = librosa.stft(x, n_fft=N, hop_length=H, win_length=N, window='hann') Fs_feature = Fs / H mag = np.abs(X) if gamma > 0: mag = np.log(1 + gamma * mag) phase = np.angle(X) / (2*np.pi) phase_diff = np.diff(phase, axis=1) phase_diff = np.concatenate((phase_diff, np.zeros((phase.shape[0], 1))), axis=1) X_hat = mag * np.exp(2*np.pi*1j*(phase+phase_diff)) X_prime = np.abs(X_hat - X) X_plus = np.copy(X_prime) for n in range(1, X.shape[1]): idx = np.where(mag[:, n] < mag[:, n-1]) X_plus[idx, n] = 0 novelty_complex = np.sum(X_plus, axis=0) if M > 0: local_average = compute_local_average(novelty_complex, M) novelty_complex = novelty_complex - local_average novelty_complex[novelty_complex < 0] = 0 if norm: max_value = np.max(novelty_complex) if max_value > 0: novelty_complex = novelty_complex / max_value return novelty_complex, Fs_feature
[docs]def resample_signal(x_in, Fs_in, Fs_out=100, norm=True, time_max_sec=None, sigma=None): """Resample and smooth signal Notebook: C6/C6S1_NoveltyComparison.ipynb Args: x_in (np.ndarray): Input signal Fs_in (scalar): Sampling rate of input signal Fs_out (scalar): Sampling rate of output signal (Default value = 100) norm (bool): Apply max norm (if norm==True) (Default value = True) time_max_sec (float): Duration of output signal (given in seconds) (Default value = None) sigma (float): Standard deviation for smoothing Gaussian kernel (Default value = None) Returns: x_out (np.ndarray): Output signal Fs_out (scalar): Feature rate of output signal """ if sigma is not None: x_in = ndimage.gaussian_filter(x_in, sigma=sigma) T_coef_in = np.arange(x_in.shape[0]) / Fs_in time_in_max_sec = T_coef_in[-1] if time_max_sec is None: time_max_sec = time_in_max_sec N_out = int(np.ceil(time_max_sec*Fs_out)) T_coef_out = np.arange(N_out) / Fs_out if T_coef_out[-1] > time_in_max_sec: x_in = np.append(x_in, [0]) T_coef_in = np.append(T_coef_in, [T_coef_out[-1]]) x_out = interp1d(T_coef_in, x_in, kind='linear')(T_coef_out) if norm: x_max = max(x_out) if x_max > 0: x_out = x_out / max(x_out) return x_out, Fs_out
[docs]def average_nov_dic(nov_dic, time_max_sec, Fs_out=100, norm=True, sigma=None): """Average respamples set of novelty functions Notebook: C6/C6S1_NoveltyComparison.ipynb Args: nov_dic (dict): Dictionary of novelty functions time_max_sec (float): Duration of output signals (given in seconds) Fs_out (scalar): Sampling rate of output signal (Default value = 100) norm (bool): Apply max norm (if norm==True) (Default value = True) sigma (float): Standard deviation for smoothing Gaussian kernel (Default value = None) Returns: nov_matrix (np.ndarray): Matrix containing resampled output signal (last one is average) Fs_out (scalar): Sampling rate of output signals """ nov_num = len(nov_dic) N_out = int(np.ceil(time_max_sec*Fs_out)) nov_matrix = np.zeros([nov_num + 1, N_out]) for k in range(nov_num): nov = nov_dic[k][0] Fs_nov = nov_dic[k][1] nov_out, Fs_out = resample_signal(nov, Fs_in=Fs_nov, Fs_out=Fs_out, time_max_sec=time_max_sec, sigma=sigma) nov_matrix[k, :] = nov_out nov_average = np.sum(nov_matrix, axis=0)/nov_num if norm: max_value = np.max(nov_average) if max_value > 0: nov_average = nov_average / max_value nov_matrix[nov_num, :] = nov_average return nov_matrix, Fs_out