1
\$\begingroup\$

I've written a working and complete proof of concept that shows a spectrograph in Matplotlib. I want to nail down this proof of concept before I continue with development, and I'm not thrilled with the way I've had to hack around using numpy.

The graph shows the spectrum as divided into five harmonics. Each harmonic section is shown centered around a harmonic multiple of the fundamental, with deviation from that center measured in cents. The crux of my problem is that these harmonic sections are of variable length. All sections are pulled as slices from a single, flat FFT output array. This output array is currently represented as a global (though that will change when I continue development), and its reference must remain intact and it must not be reallocated: this is necessary due to the way that the FFTW library works. The planner runs on a specific section of memory that must be mutated as needed.

My first approach was to get as close to a jagged array representation as numpy allows, which is masked arrays. I use a masked array for cents, the horizontal axes, so that the linspace, *= and log statements can be fully vectorized over all harmonics. It works, but doesn't seem to benefit me much since the masked array still needs to be sliced in animate before being passed to plot.set_data(). animate is more performance-critical than set_note as I plan for the former to be called for each animation frame, and the latter to only be called on key presses.

from math import sqrt, log
import matplotlib.pyplot as plt
import numpy as np

n_harmonics = 5    # Number of harmonics to display
f_upper = 24_000   # Nyquist frequency
n_fft_out = 32769  # FFT output size for an input of 64k
h_indices = np.arange(1, n_harmonics + 1)       # Harmonic indices
coefficients = np.empty(n_harmonics + 1)        # Frequency coefficients to find harmonic section bounds
coefficients[0] = n_fft_out/f_upper / sqrt(2)
coefficients[1:] = n_fft_out/f_upper * np.sqrt(h_indices*(h_indices + 1))

# A spectrum that is linear from its lower to upper frequency
# This will eventually receive an actual spectrum
fft_out = np.linspace(0, f_upper, n_fft_out, dtype=np.complex64)


def f_to_fft(f: float) -> int:
    return round(f / f_upper * n_fft_out)


def fft_to_f(i: int) -> float:
    return i / n_fft_out * f_upper


def set_note(f_tune_exact: float):
    """
    will eventually be invoked via keyboard on an irregular basis (once every few seconds)
    """
    bounds_flat = np.empty(n_harmonics + 1, dtype=np.uint32)
    np.rint(f_tune_exact * coefficients, casting='unsafe', out=bounds_flat)
    bounds = np.vstack((bounds_flat[:-1], bounds_flat[1:])).T

    sizes = (bounds[:, 1] - bounds[:, 0])[..., np.newaxis]
    longest = np.max(sizes)
    cents = np.ma.empty((n_harmonics, longest))
    cents[:, :] = np.linspace(bounds[:, 0], bounds[:, 0] + longest - 1, longest).T
    past_end = np.arange(longest)[np.newaxis, ...] >= sizes
    cents[past_end] = np.ma.masked

    cents *= (f_upper / f_tune_exact / n_fft_out / h_indices)[..., np.newaxis]
    cents = 1_200/log(2) * np.log(cents)
    return cents, bounds


def plot():
    fig, ax = plt.subplots()
    plots = [
        ax.plot([], [], label=str(h))[0]
        for h in range(1, n_harmonics + 1)
    ]

    ax.set_title('Harmonic spectrogram')
    ax.legend(title='Harmonic')
    ax.grid()
    ax.set_xlim(-600, 600)
    ax.set_ylim(0, 3_000)
    ax.set_xlabel('Deviation, cents')
    ax.set_ylabel('Spectral power')

    return plots


def animate():
    """
    Will eventually be invoked often, at the framerate of the matplotlib animation.
    fft_out needs to remain outside of this scope because FFT routines need to
    refer to the same segment of memory without it being reallocated.
    """
    for plot, x, (left, right) in zip(plots, cents, bounds):
        x = x[:right-left]
        y = np.abs(fft_out[left: right])
        plot.set_data(x, y)


if __name__ == '__main__':
    cents, bounds = set_note(440)
    plots = plot()
    animate()
    plt.show()

This shows the right thing:

spectrum

My second attempt is a trade-off: it has less vectorisation, but also a simpler animate that does not need to slice anything, since harmonics is completely stored by-reference:

from math import sqrt, log
import matplotlib.pyplot as plt
import numpy as np

n_harmonics = 5    # Number of harmonics to display
f_upper = 24_000   # Nyquist frequency
n_fft_out = 32769  # FFT output size for an input of 64k
h_indices = np.arange(1, n_harmonics + 1)       # Harmonic indices
coefficients = np.empty(n_harmonics + 1)        # Frequency coefficients to find harmonic section bounds
coefficients[0] = n_fft_out/f_upper / sqrt(2)
coefficients[1:] = n_fft_out/f_upper * np.sqrt(h_indices*(h_indices + 1))

# A spectrum that is linear from its lower to upper frequency
# This will eventually receive an actual spectrum
fft_out = np.linspace(0, f_upper, n_fft_out, dtype=np.complex64)


def f_to_fft(f: float) -> int:
    return round(f / f_upper * n_fft_out)


def fft_to_f(i: int) -> float:
    return i / n_fft_out * f_upper


def set_note(f_tune_exact: float):
    """
    will eventually be invoked via keyboard on an irregular basis (once every few seconds)
    """
    bounds_flat = np.empty(n_harmonics + 1, dtype=np.uint32)
    np.rint(f_tune_exact * coefficients, casting='unsafe', out=bounds_flat)
    bounds = np.vstack((bounds_flat[:-1], bounds_flat[1:])).T

    sizes = (bounds[:, 1] - bounds[:, 0])[..., np.newaxis]
    longest = np.max(sizes)
    cents = np.linspace(bounds[:, 0], bounds[:, 0] + longest - 1, longest).T
    cents *= (f_upper / f_tune_exact / n_fft_out / h_indices)[..., np.newaxis]
    cents = 1_200/log(2) * np.log(cents)
    cents = [
        cent[:size[0]]
        for cent, size in zip(cents, sizes)
    ]

    harmonics = [
        fft_out[left: right]
        for left, right in bounds
    ]

    return cents, harmonics


def plot():
    fig, ax = plt.subplots()
    plots = [
        ax.plot([], [], label=str(h))[0]
        for h in range(1, n_harmonics + 1)
    ]

    ax.set_title('Harmonic spectrogram')
    ax.legend(title='Harmonic')
    ax.grid()
    ax.set_xlim(-600, 600)
    ax.set_ylim(0, 3_000)
    ax.set_xlabel('Deviation, cents')
    ax.set_ylabel('Spectral power')

    return plots


def animate():
    """
    Will eventually be invoked often, at the framerate of the matplotlib animation.
    fft_out needs to remain outside of this scope because FFT routines need to
    refer to the same segment of memory without it being reallocated.
    """
    for plot, x, y in zip(plots, cents, harmonics):
        plot.set_data(x, np.abs(y))


if __name__ == '__main__':
    cents, harmonics = set_note(440)
    plots = plot()
    animate()
    plt.show()

My opinion is that the second, mask-less, by-reference, non-vectorised approach to forming harmonics and cents is more reasonable. I'm open to all feedback, but specifically what I'm really looking for is feedback on the benefits and drawbacks to these two approaches; if there is a way to use numpy/matplotlib better; and especially if there's a way to have my cake and eat it too - to have a truly jagged numpy array, which I have not been able to find.

\$\endgroup\$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.