Showcase PolyBLEP Sawtooth Oscillator

oscillator.py

import matplotlib.pyplot as plot
import numpy as np
import scipy.signal as signal


class Oscillator:
    '''
    Complex harmonic oscillator cos(phi) + 1j * sin(phi).
    '''

    phasor = 1

    def __init__(self, frequency, samplerate):

        self.frequency = frequency
        self.samplerate = samplerate

    def __call__(self, frequency=None):

        phasor = self.phasor
        omega = 2 * np.pi * (frequency or self.frequency) / self.samplerate

        self.phasor = phasor * np.exp(1j * omega)

        return phasor


class LFO(Oscillator):
    '''
    Sinusoidal wave with amplitude range [0..1].
    '''

    def __call__(self, frequency=None):

        phasor = super().__call__(frequency)

        cos = np.real(phasor)

        return (1 - cos) / 2


class OSC(Oscillator):
    '''
    The almost perfect wave, where each overtone
    amplitude is just half of the previous one.

    Designing a Pleasing Sound Mathematically
    https://www.jstor.org/stable/2690622
    '''

    def __init__(self, frequency, samplerate, shape=0.5):

        super().__init__(frequency, samplerate)

        self.shape = shape

    def __call__(self, frequency=None, shape=None):

        phasor = super().__call__(frequency)

        sin = np.imag(phasor)
        cos = np.real(phasor)

        shape = (shape or self.shape)

        return sin / (1 + (shape*shape) - (shape+shape) * cos)


class SAW(Oscillator):
    '''
    Sawtooth wave, either naive (aliased) or smoothed (antialiased).

    https://en.wikipedia.org/wiki/Sawtooth_wave

    Antialiasing Oscillators in Subtractive Synthesis
    https://ieeexplore.ieee.org/document/4117934

    Phaseshaping Oscillator Algorithms for Musical Sound Synthesis
    https://core.ac.uk/download/pdf/297014559.pdf
    http://research.spa.aalto.fi/publications/papers/smc2010-phaseshaping

    PolyBLEP Oscillators
    https://www.kvraudio.com/forum/viewtopic.php?t=375517
    '''

    t = 0

    def __init__(self, frequency, samplerate, smooth=True):

        super().__init__(frequency, samplerate)

        self.smooth = smooth

    def __call__(self, frequency=None, smooth=None):

        t = self.t
        dt = (frequency or self.frequency) / self.samplerate

        self.t = (t + dt) % 1

        smooth = (smooth or self.smooth)

        return (2 * t - 1) - (self.polyblep(t, dt) if smooth else 0)

    @staticmethod
    def polyblep(t, dt):

        assert 0 <=  t <= 1
        assert 0 <= dt <= 1

        if t < dt:
            t = t / dt
            return t+t - t*t - 1
        elif t > 1 - dt:
            t = (t - 1) / dt
            return t+t + t*t + 1
        else:
            return 0


def play(samples, samplerate):

    import pyaudio

    samples = np.atleast_1d(samples).astype(np.float32)

    audio = pyaudio.PyAudio()

    stream = audio.open(format=pyaudio.paFloat32,
                        channels=1,
                        rate=samplerate,
                        output=True)

    stream.write(samples, samples.size)

    stream.close()
    audio.terminate()


def show(samples, samplerate, spectrogram=True):

    samples = np.atleast_1d(samples)

    if spectrogram:

        from qdft import QDFT

        qdft = QDFT(samplerate, (50, samplerate/2))

        dft = qdft.qdft(np.concatenate((
            np.zeros(samplerate//2),
            samples,
            np.zeros(samplerate))))

        with np.errstate(divide='ignore'):
            db = 20 * np.log10(np.abs(dft))
            db[np.isinf(db)] = -120

        roi = (-0.5, 1 + samples.size / samplerate, 0, 1)

        args = dict(extent=roi,
                    origin='lower',
                    aspect='auto',
                    cmap='inferno',
                    interpolation='nearest')

        plot.imshow(db.T, **args)
        plot.clim(-120, 0)

        freqs = qdft.frequencies
        ticks = np.linspace(0, 1, freqs.size, endpoint=True)
        formatter = lambda tick, _: np.interp(tick, ticks, freqs).astype(int)
        plot.gca().yaxis.set_major_formatter(formatter)

        plot.xlabel('s')
        plot.ylabel('Hz')

        cbar = plot.colorbar()
        cbar.set_label('dB')

    else:

        seconds = np.arange(samples.size) / samplerate

        plot.plot(seconds, samples)
        plot.ylim(-1.1, +1.1)

        plot.xlabel('s')

    plot.tight_layout()
    plot.show()


if __name__ == '__main__':

    q = 1           # oversampling factor
    sr = 44100 * q  # sample rate in hertz
    n = 5           # signal duration in seconds

    lfo, roi  = LFO(2 / n, sr), (100, 10000)

    osc = SAW(440, sr, smooth=True)   # smooth saw
    # osc = SAW(440, sr, smooth=False)  # naive saw
    # osc = OSC(440, sr, shape=0)       # perfect sine
    # osc = OSC(440, sr, shape=0.5)     # saw like
    # osc = OSC(440, sr, shape=0.9)     # distorted sine

    y = np.array([lfo() * np.ptp(roi) + np.amin(roi) for _ in range(n * sr)])
    x = np.array([osc(_) for _ in y])

    sr = sr // q
    x = signal.decimate(x, q) if q > 1 else x
    x = 0.5 * x / np.abs(x).max()

    play(x, sr)
    show(x, sr, spectrogram=True)