public
Last active

Sethares dissmeasure function in Python

  • Download Gist
readme.md
Markdown

Adaptation of Sethares' dissonance measurement function to Python

Example is meant to match the curve in Figure 3:

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.

sethares.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
from __future__ import division
from numpy import exp, asarray, argsort, sort, sum, minimum
 
 
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.
"""
fvec = asarray(fvec)
amp = asarray(amp)
 
# 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
 
ams = amp[argsort(fvec)]
fvec = sort(fvec)
 
D = 0
for i in range(1, len(fvec)):
Fmin = fvec[:-i]
S = Dstar / (S1 * Fmin + S2)
Fdif = fvec[i:] - fvec[:-i]
if model == 'min':
a = minimum(ams[:-i], ams[i:])
elif model == 'product':
a = ams[i:] * ams[:-i] # Older model
else:
raise ValueError('model should be "min" or "product"')
Dnew = a * (C1 * exp(A1 * S * Fdif) + C2 * exp(A2 * S * Fdif))
D += sum(Dnew)
 
return D
 
 
if __name__ == '__main__':
from numpy import array, linspace, append
import matplotlib.pyplot as plt
 
# Similar to 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])
alpharange = 2.3
method = 'min'
 
# # Davide Verotta Figure 4 example
# freq = 261.63 * array([1, 2, 3, 4, 5, 6])
# amp = 1 / array([1, 2, 3, 4, 5, 6])
# alpharange = 2.0
# method = 'product'
 
# call function dissmeasure for each interval
diss = array([0])
for alpha in linspace(1, alpharange, 1000):
f = append(freq, alpha*freq)
a = append(amp, amp)
d = dissmeasure(f, a, method)
diss = append(diss, d)
 
plt.plot(linspace(1, alpharange, len(diss)), diss)
plt.xscale('log')
plt.xlim(1, alpharange)
 
plt.xlabel('frequency ratio')
plt.ylabel('dissonance')
 
for n, d in [(2,1), (3,2), (5,3), (4,3), (5,4), (6,5)]:
plt.axvline(n/d, color='silver')
plt.annotate(str(n) + ':' + str(d),
(n/d, max(diss)*2/3),
horizontalalignment='center')
 
plt.show()

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.