Last active
September 3, 2024 15:04
-
-
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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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