Sound synthesis on micropython

[Rumidom]

Hi, I'm working on a synthesizer project and wanted to write a sound synthesis library for Micropython, I went as far as writing a oscillator in pure python to produce a bytearray to be passed to a DAC over I2S. The idea being that on a dual core microcontroller (eg: pi pico) i could control the parameters of the synthesizer on one core and on a different thread I would continuously play a wavetable on the I2C DAC. I went as far as writing a small test to see if it would work. and it does! but it currently takes more time to produce the wavetable then to actually play it:

import math
import time

# Define note frequencies (middle octave - 4)
BASE_NOTES = {
    'C': 262,
    'C#': 277,
    'D': 294,
    'D#': 311,
    'E': 330,
    'F': 349,
    'F#': 370,
    'G': 392,
    'G#': 415,
    'A': 440,
    'A#': 466,
    'B': 494
}

def sine_wave(freq,x):
    return math.sin(freq * 2 * math.pi * x) 

def saw_wave(freq,x):
    value = math.atan(math.tan(freq*math.pi*x))
    return  value
    
def tri_wave(freq,x):
    return 0.6366*math.asin(math.sin(freq*2*math.pi*x))

def sqr_wave(freq,x):
    if (x * freq) % 1 < 0.5:
        return 1.0
    else:
        return -1.0
    
def generateWaveTable(sampleRate,duration,func,freq,vol):
    n_samples = sampleRate*duration
    step =  duration / (n_samples - 1)
    f_output = []
    maxint16 = 2**15-1
    for i in range(n_samples):
        t = step * i
        out = func(freq,t)*vol
        f_output.append(int(out*maxint16))
    return f_output

def generateByteWaveTable(sampleRate,n_samples,func,freq,vol):
    #n_samples = sampleRate*duration
    step =  duration / (n_samples - 1)
    output = bytearray()
    maxint16 = 2**15-1
    for i in range(n_samples):
        t = step * i
        out = func(freq,t)*vol
        output += int(out*maxint16).to_bytes(2,'little')
    return output

sampleRate = 44100
frequency = BASE_NOTES['C']
duration = 1
wavelenght = 1/frequency
n_samples = int(wavelenght*sampleRate)
waveform = sqr_wave
volume = 1

initTime = time.ticks_ms()
y_bytes = generateByteWaveTable(sampleRate,n_samples,waveform,frequency,volume)
print("Time to produce WaveTable:",time.ticks_ms()-initTime,"ms")
print("Sample Duration:",n_samples/sampleRate*1000," ms")

Output on a Pi Pico (rp2040)

Time to produce WaveTable: 27 ms
Sample Duration: 3.80952368  ms

Are there any tricks I could use to speed this up? or this is the wrong approach? is sending the I2S words directly to the DAC better? I know circuit-python has a quite extensive C library for sound synthesis called synthio but would be nice if I could do this in mycropython

[ricksorensen]

What processor are you using? Does it have a floating point processor (unit)?

I would suggest looking at I2S for the DAC process. See: https://github.com/miketeachman/micropython-i2s-examples

Also to get a more accurate time tick difference try:

time.ticks_diff(time.ticks_ms() , inittime)

this will work around overflow in the tick timer.
https://docs.micropython.org/en/latest/library/time.html#time.ticks_diff

Also https://hackaday.com/2024/08/08/tulip-is-a-micropython-synth-workstation-in-an-esp32/

Also for curiosity see:
https://hackaday.com/2024/08/08/tulip-is-a-micropython-synth-workstation-in-an-esp32/

https://www.youtube.com/watch?v=pktBirYcHD0

https://github.com/Lana-chan/Cardputer-MicroHydra/blob/main/MicroHydra/lib/beeper.py

[Rumidom]

Im using a rp4020 to test things, but I'm waiting on a rp2350 and a UDA1334A that I ordered, thanks! I'll check those links

[pi-mst]

In case you're not aware, the Tulip Creative Computer is playing in a similar space that might have a lot you can investigate. More specifically, the creator, @bwhitman, has also created the AMY music synthesizer library to power Tulip.

[bwhitman]

Thanks @pi-mst . Yes, we have spent a lot of time getting synths working in Micropython. Like Circuitpython's synthio, AMY is mostly in C with micropython wrappers. I really don't think you're going to get acceptable results in pure Micropython on current MCU hardware. If you want to write your own oscillators in Python, a shortcut could be ulab, which works great in Micropython to vectorize many operations (and hands them off to C). But AMY and Synthio are good higher-level solutions.

[Rumidom]

@bwhitman @pi-mst Thanks for the tips! I see that the ulab repository has pre-built binaries, I'll check it out. I might just end up using synthio with circuit-python if I can't get the performance on micropython

[Rumidom]

Quick update, using @micropython.native decorator and ulab's numpy trig functions, helped bring the time down to 16 ms and thats the best case scenario with a square wave, unfortunately ulab's scipy.signal doesn't have functions for signal generation, I don't think its feasible doing it in micropython.

import math
import time
from ulab import numpy as np
from ulab import scipy as sci

# Define note frequencies (middle octave - 4)
BASE_NOTES = {
    'C': 262,
    'C#': 277,
    'D': 294,
    'D#': 311,
    'E': 330,
    'F': 349,
    'F#': 370,
    'G': 392,
    'G#': 415,
    'A': 440,
    'A#': 466,
    'B': 494
}

@micropython.native
def sine_wave(freq,x):
    return math.sin(freq * 2 * np.pi * x)  #  Has frequency of 440Hz

@micropython.native
def saw_wave(freq,x):
    value = np.atan(np.tan(freq*np.pi*x))
    return  value

@micropython.native
def tri_wave(freq,x):
    return 0.6366*np.asin(np.sin(freq*2*np.pi*x))

@micropython.native
def sqr_wave(freq:int,x:ptr32):
    if (x * freq) % 1 < 0.5:
        return 1.0
    else:
        return -1.0
    
@micropython.native
def generateWaveTable(sampleRate,duration,func,freq,vol):
    n_samples = sampleRate*duration
    step =  duration / (n_samples - 1)
    f_output = []
    maxint16 = 2**15-1
    for i in range(n_samples):
        t = step * i
        out = func(freq,t)*vol
        f_output.append(int(out*maxint16))
    return f_output

@micropython.native
def generateByteWaveTable(sampleRate,n_samples,func,freq,vol):
    #n_samples = sampleRate*duration
    step =  duration / (n_samples - 1)
    output = bytearray()
    maxint16 = 2**15-1
    for i in range(n_samples):
        t = step * i
        out = func(freq,t)*vol
        output += int(out*maxint16).to_bytes(2,'little')
    return output

sampleRate = 44100
frequency = BASE_NOTES['C']
duration = 1
wavelenght = 1/frequency
n_samples = int(wavelenght*sampleRate)
waveform = sqr_wave
volume = 1

inittime = time.ticks_ms()
y_bytes = generateByteWaveTable(sampleRate,n_samples,waveform,frequency,volume)
print("Time to produce WaveTable:",time.ticks_diff(time.ticks_ms() , inittime),"ms")
print("Sample Duration:",n_samples/sampleRate*1000," ms")

[Rumidom]

ulab is cool, but kinda defeats the purpose of having a a synth library on vanila micropython.

[peterhinch]

There is a lot that can be done to improve the code in the original post. Use pre-allocated integer arrays rather than variable size lists. Use integer processing with the Viper code emitter. On RP2 consider the use of DMA to transfer the integer array contents to the I2S DAC (I'm not sure if this is possible).

Find faster algorithms!

def saw_wave(freq,x):
    value = math.atan(math.tan(freq*math.pi*x))
    return  value

Cunning, but slow. It's easy to generate saw, triangle and square waves with integer processing. It may be faster to generate an integer sin with an integer lookup table and interpolation.

However the CircuitPython library looks highly developed and is written in C. It depends whether you want an out-of-the-box solution or a project!

[GitHubsSilverBullet]

Interesting I'm just using 44100 samples/sec because I read online that this is CD quality, but makes sense that for simple tones you dont need high resolution, I will try to modify the oscillator output to with a LFO, evelope and effects tho, so playing from flash doesn't work

No, no. You'll use the full 44.1ksps. It's just that you don't need to fill that huge one second buffer.
Sometimes a few full waveforms will suffice. And then you repeat them over and over again without CPU interaction.
For instance, for ›A'‹ you'll need just four full sine waveforms adding up to 401 samples. The frequency will be 439.900 Hz and that's only 0.393 ct below target.
Even those audiophile ›golden ears‹ can't hear that.

[Rumidom]

oh I see, i'm already doing that, I wrote 'Sample Duration' but I meant 'Wave-table Duration' (I don't know much about music, I'm learning quite a lot actually), but I'm using just one full wave, the calculation you mentioned is to match the beginning and the ending of the wave-table?

[GitHubsSilverBullet]

the calculation you mentioned is to match the beginning and the ending of the wave-table?

It's an iterative process.

• Start at tone ›c0‹
• Determine the optimal length of one full waveform (most likely not an integer)
• Round to the nearest integer.
• Compute resulting frequency.
• Compute deviation in ct.
• If above 1ct, recompute with increasing number of waveforms until deviation is less than 1ct.
• Proceed with next tone until ›h8‹ has been reached

[Rumidom]

Now that you mentioned there is also a mistake on the code I posted, indeed duration is set to 1s as a global when it should be just the waveform duration, sorry for the confusion

[GitHubsSilverBullet]

Just be careful. If you stick to just one full waveform, that only works with very low frequencies.
The higher the frequencies, the more significant the deviation will be .
Up to the point were even someone with a bad hearing aid will notice the odd tones.

#waves note #oct   frequency -> frequency  #samples deviation
   1     e    8  5274.041 Hz  5512.500 Hz         8  +76.558 ct
   1     f    8  5587.652 Hz  5512.500 Hz         8  -23.442 ct
   1     f#   8  5919.911 Hz  6300.000 Hz         7 +107.732 ct
   1     g    8  6271.927 Hz  6300.000 Hz         7   +7.732 ct
   1     g#   8  6644.875 Hz  6300.000 Hz         7  -92.268 ct
   1     a    8  7040.000 Hz  7350.000 Hz         6  +74.603 ct
   1     a#   8  7458.620 Hz  7350.000 Hz         6  -25.397 ct
   1     h    8  7902.133 Hz  7350.000 Hz         6 -125.397 ct

Whereas a variable number of waveforms will produce magnitudes better notes.

#waves note #oct   frequency -> frequency  #samples deviation
  11     e    8  5274.041 Hz  5272.826 Hz        92   -0.399 ct
   9     f    8  5587.652 Hz  5590.141 Hz        71   +0.771 ct
  20     f#   8  5919.911 Hz  5919.463 Hz       149   -0.131 ct
  29     g    8  6271.927 Hz  6269.118 Hz       204   -0.776 ct
  11     g#   8  6644.875 Hz  6645.205 Hz        73   +0.086 ct
  15     a    8  7040.000 Hz  7037.234 Hz        94   -0.680 ct
  23     a#   8  7458.620 Hz  7458.088 Hz       136   -0.123 ct
  12     h    8  7902.133 Hz  7898.507 Hz        67   -0.794 ct

[Rumidom]

The code I posted previously has a mistake in it, the step size is calculated for 1 second where it should calculated for 3.8 ms, but for the proposes of comparing performance i think its ok. I did a bare minimal example with viper optimization and got down to 7 ms

import math
import time

# Define note frequencies (middle octave - 4) Hz
BASE_NOTES_HERTZ = {
    'C': 262,
    'C#': 277,
    'D': 294,
    'D#': 311,
    'E': 330,
    'F': 349,
    'F#': 370,
    'G': 392,
    'G#': 415,
    'A': 440,
    'A#': 466,
    'B': 494
}

@micropython.viper
def sqr_wave(half_w_length:int,x:int) -> int:
    if x < half_w_length:
        return int(32767)
    else:
        return int(-32767)

@micropython.native
def generateByteWaveTable(sampleRate,func,freq,vol):
    w_len = 1/frequency
    w_len_us = w_len*1000000
    half_w_len = int(round(w_len_us/2))
    n_samples = int(w_len*sampleRate)
    step =  w_len / (n_samples - 1)
    step_us = round(step*1000000)
    output = bytearray(n_samples*2)
    #output = bytearray()
    for i in range(n_samples):
        t = int(step_us * i)
        out = func(half_w_len,t)*vol
        out_bytes = out.to_bytes(2,'little')
        output[i*2] = out_bytes[0]
        output[(i*2)+1] = out_bytes[1]
    return output

sampleRate = 44100
frequency = BASE_NOTES_HERTZ['C']
wavelenght = 1/frequency
n_samples = int(wavelenght*sampleRate)
waveform = sqr_wave
volume = 1
inittime = time.ticks_us()
WaveTable_bytes = generateByteWaveTable(sampleRate,waveform,frequency,volume)
print("Time to produce WaveTable:",time.ticks_diff(time.ticks_us(), inittime)/1000,"ms")
print("Wavetable Duration:",(n_samples/sampleRate)*1000," ms")
print(WaveTable_bytes)

Time to produce WaveTable: 6.885 ms
Wavetable Duration: 3.809524  ms

I plotted the WaveTable on Jupyter notebook, and seems to be correct:

As @GitHubsSilverBullet mentioned 1/(n_samples/sampleRate) is not quite the right frequency and needs to be sampled a little longer to match, but I think this serves as a good indication. Its still not fast enough to produce the waveform faster than playing it, Its almost there and can probably be optimized further. However doesn't leave much room for adding effects on top of these simple tones