CoCalc -- 1508.sagews

All published worksheets from http://sagenb.org

Project: sagenb.org published worksheets

Views: ¹⁶⁸⁷⁴²
Image: ubuntu2004

%cython
#clib m

cdef extern from "math.h":
    float expf(float x)

# This function computes the dissonance of two partials
# at frequencies f1 and f2 (f1 < f2) with amplitudes of
# a1 and a2, respectively.
# This is using the Plomp-Levelt dissonance curve, as described in
# http://eceserv0.ece.wisc.edu/~sethares/consemi.html
def diss(float f1, float f2, float a1, float a2):
    cdef float s = (f2 - f1)*0.24/(0.0207*f1 + 18.96)
    return 5*a1*a2*(expf(-3.51*s) - expf(-5.75*s))

# discrete dissonance matrices
from itertools import product

# Bohlen Pierce just and et, with odd harmonic partials
# the outcome of this algorithm heavily depends on the parameters, see examples below
def compute_ddm(just_scale=True, number_harmonics=7, square_wave=True, base_freq=256, tritaves=2):
    # just_scale       = True  # else equal temperament
    # number_harmonics = 10    # only odd harmonics are counted
    # square_wave      = True  # else triangle
    # base_freq        = 256   # affects dissonace (plomb/levelt curve width)

    n   = tritaves*13 + 1        # range of triads
    ddm = matrix(RR, n, n)   # dissonance values

    # fill in the frequencies of the partials here
    par = [base_freq*(2*i + 1) for i in range(number_harmonics)]

    # with the respective amplitudes
    if square_wave:
        amp  = [1.0/i for i in range(1, len(par) + 1)]
    else: # triangle
        amp  = [1.0/(i**2) for i in range(1, len(par) + 1)]

    # steps, as frequency ratios
    if just_scale:
        steps = []
        _steps = [1, 27/25, 25/21, 9/7, 7/5, 75/49, 5/3, 9/5, 49/25, 15/7, 7/3, 63/25, 25/9]
        for i in range(tritaves):
            steps += [s*3**i for s in _steps]
        steps += [3**tritaves]
    else:
        steps = [float(3**(i/13)) for i in range(n)]

    pa = zip(par, amp)

    for i, j in product(range(n), repeat=2): # all triads 0, i, j
        pa2 = (pa + [(float(p*steps[i]),a) for (p,a) in pa]
                  + [(float(p*steps[j]),a) for (p,a) in pa])
        pa2.sort()
        ddm[i,j] = 0 # consonance is the absence of dissonance
        for k in range(0, len(pa2) - 1):     #  lower partial
            for l in range(k + 1, len(pa2)): # higher partial
                if (pa2[k][0] == pa2[l][0]):
                    continue
                # add up dissonance for all pairs of partials of the three notes
                ddm[i,j] += diss(pa2[k][0], pa2[l][0], pa2[k][1], pa2[l][1])
    
    return ddm

# note that the heat map values can only be interpreted
# relatively, not absolutely.
# i.e., looking at one matrix, bright/yellow is consonant
# while dark/red is more dissonant.
# but don't compare the color between one matrix and another,
# as they may be scaled differently.
#
# one could stack the matrices and plot them together, though
#
# consonance is higher for fewer harmonics with lower amplitudes
# also, the width of the Plomp/Levelt curve doesn't scale lineary,
# so the choice of the base tone frequency is important.
#
# a triad consists of three tones and denoted (i,j,k)
# where i,j & k are integers that map the steps in chromatic BP
#
# to get 2-dimensional data, we fix i at 0, the base tone of the scale
# j and k are then allowed to move between 0 and 3*13, if you choose 3 tritaves
# the value of the matrix at (j,k) then gives the dissonance of the triad (0,j,k)
#
# note that the matrix is symmetric, because the triads (0,j,k) and (0,k,j) are identical
#
# on the diagonal and the 0 row and column, two of the three tones are identical
# these are not proper triads and thus more consonant

# compute dissonance heat map for triads
#  with options:
#  just tuning
#  9 harmonics of square wave
#  base tone of triads is middle c
#  all triads, that fit in 3 tritaves above base tone
ddm  = compute_ddm(True, 9, True, 261.6, 3)
mp = matrix_plot(matrix(ddm.rows()[::-1]), cmap='hot')
mp.save("bla.png")

# compute dissonance heat map for triads
#  with options:
#  equal temperament
#  9 harmonics of square wave
#  base tone of triads is middle c
#  all triads, that fit in 3 tritaves above base tone
ddm  = compute_ddm(False, 9, True, 261.6, 3)
mp = matrix_plot(matrix(ddm.rows()[::-1]), cmap='hot')
mp.save("bla.png")

# compute dissonance heat map for triads
#  with options:
#  equal temperament
#  9 harmonics of triangle wave
#  base tone of triads is middle c
#  all triads, that fit in 3 tritaves above base tone
ddm  = compute_ddm(True, 7, True, 261.6, 3)
mp = matrix_plot(matrix(ddm.rows()[::-1]), cmap='hot')
mp.save("bla.png")

# compute dissonance heat map for triads
#  with options:
#  equal temperament
#  9 harmonics of triangle wave
#  base tone of triads is at two octaves below middle c
#  all triads, that fit in 3 tritaves above base tone
ddm  = compute_ddm(True, 7, True, 261.6/4, 3)
mp = matrix_plot(matrix(ddm.rows()[::-1]), cmap='hot')
mp.save("bla.png")

# compute dissonance heat map for triads
#  with options:
#  equal temperament
#  9 harmonics of triangle wave
#  base tone of triads is at two octaves above middle c
#  all triads, that fit in 3 tritaves above base tone
ddm  = compute_ddm(True, 7, True, 261.6*4, 3)
mp = matrix_plot(matrix(ddm.rows()[::-1]), cmap='hot')
mp.save("bla.png")

# which steps are used in the most consonant way?
#  the step ranked last here could be omitted in a 
#  12-key layout.

n = 2*13 + 1
rank = [0.0]*n

for i in range(0, n):
    for j in range(0,n):
        rank[i] += ddm[i,j]
        
ranking = []
for i in range(0, 13):
    ranking.append((rank[i], i))
ranking.sort()
for r, i in ranking:
    j = (i + 13)
    print str(i).rjust(2), str(r)[:5], "\t", 
    print str(j).rjust(2), str(rank[j])[:5], "\t",
    print str(r + rank[j])[:6]

4.580 	13 4.524 	9.1047
5.512 	23 4.144 	9.6573
5.780 	17 5.823 	11.603
5.816 	21 4.654 	10.470
5.823 	19 4.312 	10.135
6.278 	16 5.898 	12.177
6.339 	22 5.300 	11.640
6.625 	24 5.911 	12.537
6.713 	20 6.703 	13.416
7.696 	18 7.231 	14.927
8.875 	15 6.594 	15.469
10.31 	25 5.282 	15.601
24.92 	14 9.769 	34.693

# what proper triads are most/least consonant?

ddl = [(ddm[i,j],i,j) for (i,j) in product(range(n),repeat=2) if (0 < i) and (i < j)]
ddl.sort()
print "All", len(ddl), "proper triads in 9:1 region,"
print "sorted by ascending dissonance."
print

for (v,i,j) in ddl:
    print 0, str(i).rjust(2), str(j).rjust(2), ":", str(v)[:6]

WARNING: Output truncated!

full_output.txt

All 325 proper triads in 9:1 region,
sorted by ascending dissonance.

13 26 : 0.0020
13 23 : 0.0234
10 23 : 0.0262
23 26 : 0.0308
13 19 : 0.0317
19 23 : 0.0326
13 21 : 0.0335
6 19 : 0.0362
4 23 : 0.0374
8 21 : 0.0382
3 26 : 0.0476
19 21 : 0.0594
10 13 : 0.0597
19 26 : 0.0611
21 23 : 0.0655
10 26 : 0.0679
16 26 : 0.0680
3 13 : 0.0723
7 26 : 0.0744
19 25 : 0.0758
22 26 : 0.0767
19 22 : 0.0771
4 26 : 0.0774
13 16 : 0.0788
6 25 : 0.0812
21 26 : 0.0833
3 19 : 0.0843
4 25 : 0.0858
17 23 : 0.0910
10 19 : 0.0929
6 10 : 0.0930
3 22 : 0.0934
3 16 : 0.0936
4 10 : 0.0938
6 23 : 0.0942
13 17 : 0.0947
21 24 : 0.0951
9 19 : 0.0957
10 21 : 0.0967
11 19 : 0.0983
17 21 : 0.0986
8 19 : 0.0987
4 17 : 0.0989
5 26 : 0.0994
6 13 : 0.1015
4 21 : 0.1018
6 21 : 0.1022
4  8 : 0.1025
8 16 : 0.1035
11 21 : 0.1047
7 13 : 0.1052
23 25 : 0.1052
9 13 : 0.1053
15 23 : 0.1061
4 13 : 0.1062
8 23 : 0.1072

...

9 12 : 0.3434
2 16 : 0.3526
12 24 : 0.3532
8 14 : 0.3589
20 21 : 0.3590
12 22 : 0.3612
2 12 : 0.3634
25 26 : 0.3652
7 12 : 0.3684
18 19 : 0.3796
14 18 : 0.3843
23 24 : 0.3902
2 21 : 0.4033
5 12 : 0.4045
16 17 : 0.4083
17 18 : 0.4133
5 24 : 0.4212
2 14 : 0.4413
12 18 : 0.4424
19 20 : 0.4483
12 14 : 0.4648
8  9 : 0.4786
14 15 : 0.4855
9 10 : 0.4948
15 16 : 0.5269
11 12 : 0.5357
7  8 : 0.5443
10 11 : 0.5811
4  5 : 0.5959
5  6 : 0.6131
13 14 : 0.6167
3  4 : 0.6220
12 13 : 0.6490
6  7 : 0.6711
1 26 : 0.7826
1 22 : 0.7869
1 23 : 0.7958
1 24 : 0.7967
1 21 : 0.8066
1 25 : 0.8126
1 20 : 0.8251
1  9 : 0.8328
1 17 : 0.8360
1 11 : 0.8390
1 10 : 0.8422
1  7 : 0.8440
1 19 : 0.8462
1  4 : 0.8496
2  3 : 0.8541
1 16 : 0.8611
1  6 : 0.8706
1  8 : 0.8777
1 18 : 0.9031
1  3 : 0.9085
1 15 : 0.9154
1 14 : 0.9239
1  5 : 0.9355
1 13 : 0.9356
1 12 : 1.0073
1  2 : 1.3379

# the circle of fifths/fourths is the unique generator
# for the chromatic 12 tone system, (besides the half-tone)

# but 13 is prime, so any stepsize is suitable, which one to choose
# for "related" tones

# tune violine strings in 7:5, i.e. 4 steps
step = 4# et half-steps

gd = {}
for i in range(13):
    #print i, (i*step)%13
    gd[(i*step)%13] = [((i+1)*step)%13]
    
circ = Graph(gd)
circ.plot().show()
circ.plot(layout="circular").show()