Numpy wavetable "synthesizer", with a basic shapes wavetable + spectral time skew effect

numpy-wavetable-synth-spectral-time-skew.py

# Based on https://www.youtube.com/watch?v=oVi9WfoA0QI and the code https://github.com/mtytel/vital/blob/c0694a193777fc97853a598f86378bea625a6d81/src/synthesis/producers/spectral_morph.h#L132-L178

#%%
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
from scipy.fft import rfft, irfft
from IPython.display import display, Audio
import sounddevice as sd

# interpolate from sine to saw to triangle to square
def basic_shapes_wavetable(num_waveforms, waveform_length, freq=1):
    wavetable_data = np.zeros((num_waveforms, waveform_length))

    for i in range(num_waveforms):
        phase = np.arange(waveform_length) / waveform_length * freq * 2 * np.pi
        sine_wave = np.sin(phase)

        saw_wave = np.linspace(-1, 1, num=waveform_length)
        triangle_wave = np.abs(np.linspace(-1, 1, num=waveform_length * 2)[0:waveform_length]) * 2 - 1
        square_wave = np.sign(saw_wave)

        if i < num_waveforms // 3:
            wt_pos = i / (num_waveforms // 3 - 1)
            waveform = (1 - wt_pos) * sine_wave + wt_pos * saw_wave
        elif i < 2 * num_waveforms // 3:
            wt_pos = (i - num_waveforms // 3) / (num_waveforms // 3 - 1)
            waveform = (1 - wt_pos) * saw_wave + wt_pos * triangle_wave
        else:
            wt_pos = (i - 2 * num_waveforms // 3) / (num_waveforms // 3 - 1)
            waveform = (1 - wt_pos) * triangle_wave + wt_pos * square_wave

        waveform -= np.mean(waveform)  # Remove DC offset
        waveform /= np.max(np.abs(waveform))  # Normalize the waveform to the range [-1, 1]
        wavetable_data[i] = waveform

    return wavetable_data

# 'Spectral Time Skew' oscillator effect that scrolls through the wavetable a different amount for every harmonic. 
def apply_spectral_skew(wavetable, wavetable_pos, skew, freq=1):
    num_waveforms, waveform_length = wavetable.shape
    fft_size = waveform_length * 2

    wt_index = int(wavetable_pos * (num_waveforms - 1))
    max_skew_index = num_waveforms - 1 - wt_index
    skew_index = int(skew * max_skew_index)
    skewed_index = min(wt_index + skew_index, num_waveforms - 1)

    base_waveform = wavetable[wt_index]
    padding_size = waveform_length // 2
    zeros = np.zeros((wavetable.shape[0], padding_size))
    fft_wavetable = np.concatenate((zeros, wavetable, zeros), axis=-1)
    fft_wavetable = rfft(fft_wavetable, axis=-1)

    skewed_waveform = np.zeros_like(fft_wavetable[0])
    skewed_waveform[0] = base_waveform[0]

    num_harmonics = 256
    max_harmonic = np.min([fft_wavetable.shape[-1]-1, num_harmonics])

    for i in range(max_harmonic):
        if i == 0:
            continue
        skew_amount = skew_index * i / (waveform_length // 2)
        skewed_index = int(min(wt_index + skew_amount, num_waveforms - 1))
        from_amplitudes = fft_wavetable[skewed_index-1, :] if skewed_index > 0 else fft_wavetable[wt_index]
        to_amplitudes = fft_wavetable[skewed_index, :]
        t = skew_amount % 1.0

        real = np.interp(t, [0, 1], [from_amplitudes.real[i], to_amplitudes.real[i]])
        imag = np.interp(t, [0, 1], [from_amplitudes.imag[i], to_amplitudes.imag[i]])

        skewed_waveform[i] = real + 1j*imag
    
    skewed_waveform[max_harmonic:] = 0

    time_domain_waveform = irfft(skewed_waveform)[padding_size:-padding_size]
    time_domain_waveform /= np.max(np.abs(time_domain_waveform))

    # truncate
    time_domain_waveform = time_domain_waveform[:waveform_length]

    return time_domain_waveform


def plot_wavetable(wavetable, waveform_index):
    fig, ax = plt.subplots()
    wave = wavetable[waveform_index]
    ax.set_ylim(-1, 1)
    ax.plot(wave)
    plt.show()

def update(val):
    wt_pos = wt_pos_slider.val
    skew = skew_slider.val
    frequency = freq_slider.val # added line

    wave = apply_spectral_skew(wavetable_data, wt_pos, skew, frequency) # modified line
    line.set_ydata(wave)
    fig.canvas.draw_idle()

    play_audio(wave, sample_rate, frequency)

def play_waveform(waveform, frequency, duration):
    num_samples_in_waveform = len(waveform)
    increment = (num_samples_in_waveform * frequency) / sample_rate

    # Generate samples for the desired duration
    num_output_samples = int(duration * sample_rate)

    positions = (np.arange(num_output_samples) * increment) % num_samples_in_waveform
    positions_int = positions.astype(int)
    positions_frac = positions - positions_int

    sample_1 = waveform[positions_int]
    sample_2 = waveform[(positions_int + 1) % num_samples_in_waveform]

    output = (1 - positions_frac) * sample_1 + positions_frac * sample_2

    output /= np.max(np.abs(output))

    return output

def play_audio(waveform, sample_rate, frequency=110):
    duration = 1.0  # seconds
    repeated_waveform = np.concatenate([waveform] * int(sample_rate * duration / len(waveform)))
    repeated_waveform = play_waveform(waveform, frequency, duration)
    audio = Audio(repeated_waveform, rate=sample_rate)
    sd.play(repeated_waveform, sample_rate)
    display(audio)

fig, ax = plt.subplots()
ax.set_ylim(-1, 1)
plt.subplots_adjust(bottom=0.25)

num_waveforms = 128
waveform_length = 256
frequency = 55
wavetable_pos = 0.3 # wavetable progress
skew = 0.4 # how far to scan the wavetable for non-fundamental harmonics

wavetable_data = basic_shapes_wavetable(num_waveforms, waveform_length)

wave = apply_spectral_skew(wavetable_data, wavetable_pos=wavetable_pos, skew=skew)

line, = ax.plot(wave)

wt_pos_slider_ax = plt.axes([0.2, 0.1, 0.65, 0.03])
wt_pos_slider = Slider(wt_pos_slider_ax, 'Wavetable', 0.0, 1.0, valinit=0)

skew_slider_ax = plt.axes([0.2, 0.05, 0.65, 0.03])
skew_slider = Slider(skew_slider_ax, 'Skew', 0.0, 1.0, valinit=0)

freq_slider_ax = plt.axes([0.2, 0.15, 0.65, 0.03])
freq_slider = Slider(freq_slider_ax, 'Frequency', 20, 1000, valinit=frequency)

wt_pos_slider.on_changed(update)
skew_slider.on_changed(update)
freq_slider.on_changed(update)

sample_rate = 44100
play_audio(wave, sample_rate, frequency)

plt.show()