Generators for sawtooth waveforms using various methods

gen.py

print('************************************************')
print('pip must first be used to install the following:')
print('pip3 install numpy simpleaudio')
print('************************************************')
print()

import numpy as np
import simpleaudio as sa
import wave, struct, math, random

# Geneartor params
toneFreqHz = 1760
fs = 44100
seconds = 2

###################################
modes = ['genTime', 'genSines', 'genPolyBlep']
genMode = modes[2] # Edit this index to change the generation method
alsoPlayAudio = True # Or just write to WAV file on disk
###################################

MAX_SHORT = 2**15 - 1;
NYQUIST = fs / 2

# Write to wav file
filename = f"python_{toneFreqHz}Hz_{genMode}_tone.wav"

wav = wave.open(filename,'w')
wav.setnchannels(1) # mono
wav.setsampwidth(2)
wav.setframerate(fs)

####################################################################################

out = []
phase = 0
sampleTime = 1.0 / fs
minVal = -0.25
maxVal = minVal * -1

if (genMode == modes[0]):
	# Generate a saw wave directly in the time domain, no band limiting
	# Essentially this is a clone of SawOsc in audiograph.js:
	# https://github.com/maximecb/noisecraft/blob/main/public/audiograph.js#L342
	for n in range(seconds * fs):
		phase += sampleTime * toneFreqHz
		cyclePos = phase % 1 # Values are in range 0..1
		valFloat = minVal + cyclePos * (maxVal - minVal)
		out.append(valFloat)

		# Write to wav file on disk, but in short int format
		# Packed data range is then [-32767, 32767] because of minVal/maxVal above
		data = struct.pack('<h', (int)(valFloat * MAX_SHORT))
		wav.writeframesraw(data)

if (genMode == modes[1]):
	# Generate a saw wave by adding together sin waves of various frequencies
	# Formula per Wikipedia, Xsawtooth(t):
	# https://en.wikipedia.org/wiki/Sawtooth_wave
	n_tones = 20 # Trade-off between saw wave completeness vs. compute time
	for n in range(seconds * fs):
		phase += sampleTime * toneFreqHz
		sample = 0;
		for k in range(1, n_tones + 1):
			freq = k * toneFreqHz
			if (freq > NYQUIST): continue; # Skip frequencies above nyquist to avoid aliasing
			valFloat = math.sin(k * phase * 2 * math.pi) / k; # valFloat is in range -1..1
			valFloat *= ((-1) ** (k + 1))
			sample += valFloat
		# TODO: fix normalization to properly scale from min to max values
		sample = minVal + ((sample + 1) / 2) * (maxVal - minVal)
		out.append(sample)

		data = struct.pack('<h', (int)(sample * MAX_SHORT))
		wav.writeframesraw(data)

def poly_blep(phase, phaseInc):
	if (phase < phaseInc):
		phase /= phaseInc
		return phase+phase - (phase*phase) - 1
	elif (phase > (1 - phaseInc)):
		phase = (phase - 1) / phaseInc;
		return (phase*phase) + phase+phase + 1
	else:
		return 0

if (genMode == modes[2]):
	# Generate a saw wave by using PolyBLEP method, described here:
	# https://www.metafunction.co.uk/post/all-about-digital-oscillators-part-2-blits-bleps
	# https://www.kvraudio.com/forum/viewtopic.php?t=375517
	phaseInc = sampleTime * toneFreqHz # Equivalent to dt, or freq / sample_rate
	for n in range(seconds * fs):
		adjustedPhase = phase + 0.5
		if (adjustedPhase >= 1): adjustedPhase = adjustedPhase % 1

		output = 2 * adjustedPhase - 1 # Raw saw wave, range from -1..1
		offset = poly_blep(adjustedPhase, phaseInc)
		output = minVal + (output - offset + 1) / 2 * (maxVal - minVal)
		out.append(output)

		# Incredibly, we can see that only the two samples around the discontinuity
		# are modified at all. Small change with a big effect.
		if (offset != 0): print(f"output: {output}, offset: {offset}")
		else: print(f"output: {output}")

		data = struct.pack('<h', (int)(output * MAX_SHORT))
		wav.writeframesraw(data)

		phase += phaseInc

####################################################################################

wav.close()

if alsoPlayAudio:
	# Ensure that highest value is in 16-bit range
	audio = np.array(out) * (MAX_SHORT)

	# Convert to 16-bit data
	audio = audio.astype(np.int16)

	# Start playback
	play_wav = sa.play_buffer(audio, 1, 2, fs)

	# Wait for playback to finish before exiting
	play_wav.wait_done()