Skip to content

Instantly share code, notes, and snippets.

@lynzrand
Last active September 3, 2024 15:04
Show Gist options
  • Save lynzrand/e65c777d501289ae3876b59b27bd0f62 to your computer and use it in GitHub Desktop.
Save lynzrand/e65c777d501289ae3876b59b27bd0f62 to your computer and use it in GitHub Desktop.
Dissonance calculation for up to 19 harmonics with microtonic scales shown -- also see https://gist.github.com/endolith/3066664
"""
Python translation of http://sethares.engr.wisc.edu/comprog.html
"""
from typing import Callable, List, Tuple
import numpy as np
import math
MAX_HARMONICS = 19
MAX_HARMONICS_SHOWN = 11
EDS_SHOWN = [(12, 2), (17, 2), (19, 2), (13, 3)] # 12edo, 17edo, 19edo, 13ed3
SAMPLES = 700
def dissmeasure(fvec, amp, model='min'):
"""
Given a list of partials in fvec, with amplitudes in amp, this routine
calculates the dissonance by summing the roughness of every sine pair
based on a model of Plomp-Levelt's roughness curve.
The older model (model='product') was based on the product of the two
amplitudes, but the newer model (model='min') is based on the minimum
of the two amplitudes, since this matches the beat frequency amplitude.
"""
# Sort by frequency
sort_idx = np.argsort(fvec)
am_sorted = np.asarray(amp)[sort_idx]
fr_sorted = np.asarray(fvec)[sort_idx]
# Used to stretch dissonance curve for different freqs:
Dstar = 0.24 # Point of maximum dissonance
S1 = 0.0207
S2 = 18.96
C1 = 5
C2 = -5
# Plomp-Levelt roughness curve:
A1 = -3.51
A2 = -5.75
# Generate all combinations of frequency components
idx = np.transpose(np.triu_indices(len(fr_sorted), 1))
fr_pairs = fr_sorted[idx]
am_pairs = am_sorted[idx]
Fmin = fr_pairs[:, 0]
S = Dstar / (S1 * Fmin + S2)
Fdif = fr_pairs[:, 1] - fr_pairs[:, 0]
if model == 'min':
a = np.amin(am_pairs, axis=1)
elif model == 'product':
a = np.prod(am_pairs, axis=1) # Older model
else:
raise ValueError('model should be "min" or "product"')
SFdif = S * Fdif
D = np.sum(a * (C1 * np.exp(A1 * SFdif) + C2 * np.exp(A2 * SFdif)))
return D
def create_amp_square(n_harmonics) -> Tuple[np.ndarray, np.ndarray]:
"""
Return the (frequency, amplitude) pairs for a square wave
"""
freq = SAMPLES * np.array(range(1, n_harmonics + 1, 2))
amp = 1 / np.array(range(1, n_harmonics + 1, 2))
return freq, amp
def create_amp_sawtooth(n_harmonics) -> Tuple[np.ndarray, np.ndarray]:
"""
Return the (frequency, amplitude) pairs for a sawtooth wave
"""
freq = SAMPLES * np.array(range(1, n_harmonics + 1))
amp = 1 / np.array(range(1, n_harmonics + 1))
return freq, amp
def create_amp_exponential(n_harmonics) -> Tuple[np.ndarray, np.ndarray]:
"""
Return the (frequency, amplitude) pairs for a waveform that has exponentially
decreasing amplitudes for harmonics
"""
freq = SAMPLES * array(range(1, n_harmonics + 2))
amp = 0.88**array(range(0, n_harmonics + 1))
return freq, amp
def calculate_dissonance(freq_amp_functions: List[Callable], r_low: float,
alpharange: float, method: str, n: int) -> np.ndarray:
diss_all = []
for func in freq_amp_functions:
freq, amp = func(MAX_HARMONICS)
diss = empty(n)
a = concatenate((amp, amp))
for i, alpha in enumerate(linspace(r_low, alpharange, n)):
f = concatenate((freq, alpha * freq))
d = dissmeasure(f, a, method)
diss[i] = d
diss_all.append(diss)
diss_all = np.array(diss_all)
diss_all = diss_all / np.max(diss_all, axis=1)[:, None]
return diss_all
if __name__ == '__main__':
from numpy import array, linspace, empty, concatenate
import matplotlib.pyplot as plt
"""
Reproduce Sethares Figure 3
http://sethares.engr.wisc.edu/consemi.html#anchor15619672
"""
freq_amp_functions = [
create_amp_square, create_amp_sawtooth, create_amp_exponential
]
r_low = 1
alpharange = 3.1
method = 'min'
n = 3000
diss_all = calculate_dissonance(freq_amp_functions, r_low, alpharange,
method, n)
fig, ax = plt.subplots(figsize=(35, 5))
for i, diss in enumerate(diss_all):
p1 = ax.plot(linspace(r_low, alpharange, len(diss)),
diss,
label=freq_amp_functions[i].__name__)
plt.xscale('log')
plt.legend()
plt.xlim(r_low, alpharange)
plt.xlabel('frequency ratio')
plt.ylabel('sensory dissonance')
intervals = [(x + y, y) for y in range(1, MAX_HARMONICS_SHOWN + 1)
for x in range(0,
int((alpharange - r_low) * y) + 1)
if math.gcd(x, y) == 1]
for n, d in intervals:
plt.axvline(n / d, color='silver')
plt.yticks([])
plt.minorticks_off()
plt.xticks([n / d for n, d in intervals],
['{}/{}'.format(n, d) for n, d in intervals])
plt.tick_params('x', labelrotation=45)
# Add a secondary linear axis at the bottom, showing edo intervals
ax2_pos = 30
for divisions, size in EDS_SHOWN:
ax2 = ax.twiny()
ticks = []
t0 = 0
while size**(t0 / divisions) < alpharange:
ticks.append(t0)
t0 += 1
ax2.set_xscale('log')
ax2.set_xticks([size**((n) / divisions) for n in ticks])
if size == 2:
ax2.set_xticklabels([f"{n}\\{divisions}" for n in ticks])
else:
ax2.set_xticklabels([f"{n}\\{divisions}ed{size}" for n in ticks])
ax2.set_xticks([n / d for n, d in intervals], minor=True)
ax2.tick_params(axis='x',
which='minor',
labelsize=0,
labelcolor='none',
color='gray',
tickdir='in')
ax2.xaxis.set_ticks_position('bottom')
ax2.set_xlim(left=r_low, right=alpharange)
# ax2.xaxis.set_label_position('left')
# ax2.set_xlabel(f"{divisions}ed{size}")
ax2_pos += 25
ax2.spines['bottom'].set_position(('outward', ax2_pos))
plt.tight_layout()
plt.savefig("out.png")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment