Librosa issue #891

librosa-891.py

import librosa
import numpy as np

y, sr = librosa.load(librosa.util.example_audio_file())
y = y[:sr] * (2**31)
S = librosa.feature.melspectrogram(y, sr=sr).astype('float32')
G = librosa.pcen(S, sr=sr, power=1, bias=0).astype('float32')

powers = np.logspace(-12, 0, 13)
biases = [1.0, 2.0, 5.0, 10.0]

fig, axs = plt.subplots(1, len(biases), figsize=(16, 4))

for bias, ax in zip(biases, axs.ravel()):
    
    in_rel_diffs_v1, in_rel_diffs_v2 = [], []
    out_rel_diffs = []
    
    for power in powers:
        # Pow-based
        PCEN_v1 = (G + bias)**power - bias**power
        G_tilde_v1 = (PCEN_v1 + bias**power)**(1/power) - bias
        in_rel_diffs_v1.append(np.max(np.abs(G - G_tilde_v1) / G))

        # Expm1-based exponentiation 
        PCEN_v2 = bias**power * np.expm1(power * np.log1p(G / bias))
        G_tilde_v2 = bias * np.expm1(np.log1p(PCEN_v2/(bias**power)) / power)
        in_rel_diffs_v2.append(np.max(np.abs(G - G_tilde_v2) / G))
        
        out_rel_diffs.append(np.max(np.abs(PCEN_v2 - PCEN_v1) / PCEN_v2))

    ax.plot(np.log10(powers), np.log10(in_rel_diffs_v1), '-o', label='pow-based')
    ax.plot(np.log10(powers), np.log10(in_rel_diffs_v2), '-o', label='expm1-based')
    ax.set_xlabel('log10(power)')
    ax.set_ylabel('log10(relative error)')
    ax.set_title("Bias = {:5.2f}".format(bias))
    ax.set_xticks(np.arange(-12, 2, 2))
    ax.set_yticks(np.arange(-8, 8, 2))
    ax.set_ylim(-8, 6)

plt.legend()

plt.savefig('librosa-890_a.png', dpi=500, bbox_inches='tight')


fig, axs = plt.subplots(1, len(biases), figsize=(16, 4))

for bias, ax in zip(biases, axs.ravel()):
    ax.plot(np.log10(powers), np.log10(out_rel_diffs), '-s', color='g', label='pow vs. expm1')
    ax.set_xlabel('log10(power)')
    ax.set_ylabel('log10(relative difference)')
    ax.set_title("Bias = {:5.2f}".format(bias))
    ax.set_xticks(np.arange(-12, 2, 2))
    ax.set_yticks(np.arange(-8, 8, 2))
    ax.set_ylim(-8, 6)

plt.legend()

plt.savefig('librosa-890_b.png', dpi=500, bbox_inches='tight')