Instantly share code, notes, and snippets.

Last active May 7, 2024 08:38
Show Gist options
• Save endolith/3066664 to your computer and use it in GitHub Desktop.
Sethares dissmeasure function in Python

Adaptation of Sethares' dissonance measurement function to Python

Example is meant to match the curve in Figure 3:

Original model used products of the two amplitudes a1⋅a2, but this was changed to minimum of the two amplitudes min(a1, a2), as explained in G: Analysis of the Time Domain Model appendix of Tuning, Timbre, Spectrum, Scale.

This weighting is incorporated into the dissonance model (E.2) by assuming that the roughness is proportional to the loudness of the beating. ... Thus, the amplitude of the beating is given by the minimum of the two amplitudes.

With the first 6 harmonics at amplitudes 1/n starting at 261.63 Hz, using the product model, it also perfectly matches Figure 4 of Davide Verotta - Dissonance & Composition, so it should be trustworthy.

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 """ import numpy as np 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 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 = 500 * array([1, 2, 3, 4, 5, 6]) amp = 0.88**array([0, 1, 2, 3, 4, 5]) r_low = 1 alpharange = 2.3 method = 'product' # # Davide Verotta Figure 4 example # freq = 261.63 * array([1, 2, 3, 4, 5, 6]) # amp = 1 / array([1, 2, 3, 4, 5, 6]) # r_low = 1 # alpharange = 2.0 # method = 'product' n = 3000 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 plt.figure(figsize=(7, 3)) plt.plot(linspace(r_low, alpharange, len(diss)), diss) plt.xscale('log') plt.xlim(r_low, alpharange) plt.xlabel('frequency ratio') plt.ylabel('sensory dissonance') intervals = [(1, 1), (6, 5), (5, 4), (4, 3), (3, 2), (5, 3), (2, 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.tight_layout() plt.show()

### mauro-belgiovine commented Sep 25, 2023 • edited

Just modified the code to calculate up to 19th harmonic and it looks great for microtonal usages. Thanks!

@lynzrand Sorry to resurrect this thread, but this is not available anymore. I'd be curious to check your implementation, thanks!

@endolith Thanks for this translated code.

### lynzrand commented Sep 26, 2023 • edited

@lynzrand Sorry to resurrect this thread, but this is not available anymore. I'd be curious to check your implementation, thanks!

Hi! The old code repository was accidentally deleted at some time by me, probably because I was cleaning up my replit account. The implementation is just a simple edit from the code in the original code.

Anyway, the code should be available again in the same repository (https://replit.com/@01010101lzy/Dissonance). I also created a gist for it: https://gist.github.com/lynzrand/e65c777d501289ae3876b59b27bd0f62

Update: I have updated my code to also contain (microtonal) scales as a reference, and also different frequency functions!

### mauro-belgiovine commented Sep 26, 2023

@lynzrand Sorry to resurrect this thread, but this is not available anymore. I'd be curious to check your implementation, thanks!

Hi! The old code repository was accidentally deleted at some time by me, probably because I was cleaning up my replit account. The implementation is just a simple edit from the code in the original code.

Anyway, the code should be available again in the same repository (https://replit.com/@01010101lzy/Dissonance). I also created a gist for it: https://gist.github.com/lynzrand/e65c777d501289ae3876b59b27bd0f62

Update: I have updated my code to also contain (microtonal) scales as a reference, and also different frequency functions!

@lynzrand thank you so much for this!!