Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Band-limited waveform generation in Python/SciPy

Returns periodic band-limited square, triangle, or sawtooth waveform with period 2π

Meant to be called the same way as scipy.signal.waveforms.square, but without the aliasing effects produced by the existing "naïve" implementation. Main goal is numerical accuracy for doing simulations, so the wave is generated by additive synthesis.

At low frequencies, the difference from the naive version is minimal, and the computational cost of the additive synthesis is high.

At high frequencies, the difference is very obvious, and the computational cost of additive synthesis is not too bad.

It currently requires a constant sampling frequency (no chirps) and doesn't support width or duty cycle.

Usage:

t = linspace(0, 1, num = 1000, endpoint = False)
plot(bl_square(2 * pi * 3 * t))

Waveform:

Waveform

Band-limited version in blue, naïve version in red (equivalent to sampling the ideal mathematical function without putting it through an antialiasing filter first)

Spectra:

Spectra

The extra spikes are the aliased harmonics, which are even more obvious when log-scaled:

dB-scaled spectra:

dB-scaled spectra

Sound examples:

TODO:

Something like BLEPs might be better for this? (Probably not; they are not as accurate due to truncation of the sinc.)

Or the discrete summation methods would be more precise? Perhaps, but they cannot generate sawtooth/triangle/square directly, because the harmonic fall-off is at a different rate. They can generate BLITs, which can then be integrated to form the others, but this accumulates numerical error and requires leaky integration. Not sure how accurate that is.

With sweeps or chirps, there is also the problem of harmonics just disappearing suddenly when they hit fs/2 and producing sudden amplitude changes, which can be fixed with a roll-off filter, but which roll-off filter? It should probably be a parameter? But how to specify it as a parameter? I could also just cut it off suddenly, and let the user do their own filter, or oversample and decimate, if they want to. That's probably best.

No, you can't cut off the sweep frequencies so that one sample has a frequency and the next does not and then filter after the fact. The transition will cause a glitch with infinite harmonics that will alias, and post-filtering won't help. So it needs some kind of filter= argument for specifying the roll-off filter. Would it specify a filter object? and then harmonics are scaled by its frequency response? And filter=None should be possible which represents a brickwall! Then it will click for sweeps or go right up to Nyquist for non-sweeps if that's what's desired.

Also should have a phase={'zero', 'minimum', 'filter'} option.

# -*- coding: utf-8 -*-
"""
Bandlimited versions of scipy.signal.waveforms.
Intent is mathematical perfection over performance;
these use additive synthesis, so they are slow, but exact.
Less ideal methods using BLIT:
Sawtooth can be made by integrating BLIT minus a DC value to prevent integrator wandering off
Square can be made by integrating bipolar BLIT
Triangle can be made by integrating square
But I'm trying to avoid leaky integration, etc.
"""
from __future__ import division
from numpy import asarray, zeros, pi, sin, cos, amax, diff, arange, outer
# TODO: De-duplicate all this code into one generator function and then make these into wrapper functions
def bl_sawtooth(t): # , width=1
"""
Return a periodic band-limited sawtooth wave with
period 2*pi which is falling from 0 to 2*pi and rising at
2*pi (opposite phase relative to a sin)
Produces the same phase and amplitude as scipy.signal.sawtooth.
Examples
--------
>>> t = linspace(0, 1, num = 1000, endpoint = False)
>>> f = 5 # Hz
>>> plot(bl_sawtooth(2 * pi * f * t))
"""
t = asarray(t)
if abs((t[-1]-t[-2]) - (t[1]-t[0])) > .0000001:
raise ValueError("Sampling frequency must be constant")
if t.dtype.char in ['fFdD']:
ytype = t.dtype.char
else:
ytype = 'd'
y = zeros(t.shape, ytype)
# Get sampling frequency from timebase
fs = 1 / (t[1] - t[0])
# fs = 1 / amax(diff(t))
# Sum all multiple sine waves up to the Nyquist frequency
# TODO: Maybe choose between these based on number of harmonics?
# Slower, uses less memory
for h in range(1, int(fs*pi)+1):
y += 2 / pi * -sin(h * t) / h
# Faster, but runs out of memory and dies
# h = arange(1, int(fs * pi) + 1)
# phase = outer(t, h)
# y = 2 / pi * -sin(phase) / h
# y = sum(y, axis=1)
return y
def bl_triangle(t):
"""
Return a periodic band-limited triangle wave with
period 2*pi which is falling from 0 to pi and rising from
pi to 2*pi (same phase as a cos)
Produces the same phase and amplitude as scipy.signal.sawtooth(width=0.5).
Examples
--------
>>> t = linspace(0, 1, num = 1000, endpoint = False)
>>> f = 5 # Hz
>>> plot(bl_triangle(2 * pi * f * t))
"""
t = asarray(t)
if abs((t[-1]-t[-2]) - (t[1]-t[0])) > .0000001:
raise ValueError("Sampling frequency must be constant")
if t.dtype.char in ['fFdD']:
ytype = t.dtype.char
else:
ytype = 'd'
y = zeros(t.shape, ytype)
# Get sampling frequency from timebase
fs = 1 / (t[1] - t[0])
# Sum all odd multiple sine waves up to the Nyquist frequency
# Slower, uses less memory
for h in range(1, int(fs * pi) + 1, 2):
y += 8 / pi**2 * -cos(h * t) / h**2
# Faster, but runs out of memory and dies
# h = arange(1, int(fs * pi) + 1, 2)
# phase = outer(t, h)
# y = 8 / pi**2 * -cos(phase) / h**2
# y = sum(y, axis=1)
return y
def bl_square(t, duty=0.5):
"""
Return a periodic band-limited square wave with
period 2*pi which is +1 from 0 to pi and -1 from
pi to 2*pi (same phase as a sin)
Produces the same phase and amplitude as scipy.signal.square.
Similarly, duty cycle can be set, or varied over time.
Examples
--------
>>> t = linspace(0, 1, num = 10000, endpoint = False)
>>> f = 5 # Hz
>>> plot(bl_square(2 * pi * f * t))
>>> sig = np.sin(2 * np.pi * t)
>>> pwm = bl_square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
>>> plt.subplot(2, 1, 1)
>>> plt.plot(t, sig)
>>> plt.subplot(2, 1, 2)
>>> plt.plot(t, pwm)
>>> plt.ylim(-1.5, 1.5)
"""
return bl_sawtooth(t - duty*2*pi) - bl_sawtooth(t) + 2*duty-1
def blit(t):
"""
Return a periodic band-limited impulse train (Dirac comb) with
period 2*pi (same phase as a cos)
Examples
--------
>>> t = linspace(0, 1, num = 1000, endpoint = False)
>>> f = 5.4321 # Hz
>>> plot(blit(2 * pi * f * t))
References
----------
http://www.music.mcgill.ca/~gary/307/week5/bandlimited.html
"""
t = asarray(t)
if abs((t[-1]-t[-2]) - (t[1]-t[0])) > .0000001:
raise ValueError("Sampling frequency must be constant")
if t.dtype.char in ['fFdD']:
ytype = t.dtype.char
else:
ytype = 'd'
y = zeros(t.shape, ytype)
# Get sampling frequency from timebase
fs = 1 / (t[1] - t[0])
# Sum all multiple sine waves up to the Nyquist frequency
N = int(fs * pi) + 1
for h in range(1, N):
y += cos(h * t)
y /= N
# h = arange(1, int(fs * pi) + 1)
# phase = outer(t, h)
# y = 2 / pi * cos(phase)
# y = sum(y, axis=1)
return y
if __name__ == "__main__":
from numpy import linspace
import matplotlib.pyplot as plt
from scipy.signal import square, sawtooth
fs = 500
t = linspace(0, 1, num = fs, endpoint = False)
f = 5.432 # Hz
# Test that waveforms match SciPy versions
plt.figure()
plt.subplot(3, 1, 1)
plt.plot(t, square(2 * pi * f * t), color='gray')
plt.plot(t, bl_square(2 * pi * f * t))
plt.title('Square')
plt.subplot(3, 1, 2)
plt.plot(t, sawtooth(2 * pi * f * t, 0.5), color='gray')
plt.plot(t, bl_triangle(2 * pi * f * t))
plt.title('Triangle')
plt.subplot(3, 1, 3)
plt.plot(t, sawtooth(2 * pi * f * t), color='gray')
plt.plot(t, bl_sawtooth(2 * pi * f * t))
plt.title('Sawtooth')
# Square wave duty cycle test
plt.figure()
plt.subplot(3, 1, 1)
width = 0.5
plt.plot(t, square(2*pi*f*t, width), color='gray')
plt.plot(t, bl_square(2*pi*f*t, width))
plt.margins(0, 0.1)
plt.subplot(3, 1, 2)
width = 0.01
plt.plot(t, square(2*pi*f*t, width), color='gray')
plt.plot(t, bl_square(2*pi*f*t, width))
plt.margins(0, 0.1)
plt.subplot(3, 1, 3)
width = 2/3
plt.plot(t, square(2*pi*f*t, width), color='gray')
plt.plot(t, bl_square(2*pi*f*t, width))
plt.margins(0, 0.1)
@wdputro

This comment has been minimized.

Copy link

commented Mar 28, 2019

Hi,
Thanks for sharing this Gist. I'm pretty new to the topic, but do you have any advice on how to create the additive version of the Spectra image you upload with Additive Synthesis?

Thanks,

@endolith

This comment has been minimized.

Copy link
Owner Author

commented Apr 10, 2019

@wdputro I don't understand the question. This is additive synthesis.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.