1407.sagews

%cython
#clib m

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

def diss_curve(par, amp, float ran=4.0, int res=1200):
    x = [0.0]*res
    y = [0.0]*res
    cdef float step = (ran - 1.0)/res
    cdef int n = len(par)
    cdef int xx
    cdef int i
    cdef int j
    cdef float dif
    cdef float s
    cdef float a
    cdef float pai
    cdef float paj
    cdef float xxx
    cdef float yxx
    for xx in range(res):
        xxx = 1.0 + xx*step
        yxx = 0.0
        for i in range(2*n):
            if i < n:
                pai = par[i]
            else:
                pai = xxx*par[i-n]
            for j in range(2*n):
                if j < n:
                    paj = par[j]
                else:
                    paj = xxx*par[j-n]
                dif = paj - pai
                if dif >= 0:
                    s = dif*0.24/(0.0207*pai + 18.96)
                    a = amp[i%n]*amp[j%n]
                    yxx += 5*a*(expf(-3.51*s) - expf(-5.75*s))
        x[xx] = xxx
        y[xx] = yxx
    return x, y

Auto-generated code...  

def plot_diss_curve(x, y, steps_12=False, steps_bp=False):
    DC = list_plot(zip(x,y), plotjoined=True)
    xmax = x[-1]
    ymax = max(y)
    # also plot the steps of 12-tet
    if steps_12:
        N = int(12*log(xmax)/log(2)) + 1
        for i in range(N):
            xx = 2**(i/12)
            DC += line([(xx, 0), (xx,ymax)], color='green')
    # or the steps of bohlen-pierce-et
    if steps_bp:
        N = int(13*log(xmax)/log(3)) + 1
        for i in range(N):
            xx = 3**(i/13)
            DC += line([(xx, 0), (xx,ymax)], color='red')
    return DC

def plot_partials(par, amp):
    n = int(max(par)/min(par) + 1)
    P = list_plot(zip(par, amp))
    for i in range(1, n+1):
        P += line([(i*par[0], 0), (i*par[0],1)], color='gray')
    return P

r = factorial(6)
print r

720

# 12 tet, with harmonic partials

# fill in the frequencies of the partials here
par = [440*i for i in [1, 2, 3, 4, 5, 6]]

# with the respective amplitudes
amp  = [0.8**i for i in range(len(par))]

# compute coordinates of curve
x, y = diss_curve(par, amp, ran=2.2, res=r)

# show partials as points, with harmonics as reference
plot_partials(par, amp).show()

# and plot dissonance curve with scale steps
plot_diss_curve(x, y, steps_12=True, steps_bp=False).show()

# 12-tet, with non-harmonic partials

# fill in the frequencies of the partials here
par = [440*i for i in [1, 2, 2**(19/12), 4, 2**(28/12), 2**(31/12)]]

# with the respective amplitudes
amp  = [0.8**i for i in range(len(par))]

# compute coordinates of curve
x, y = diss_curve(par, amp, ran=2.2, res=r)

# show partials as points, with harmonics as reference
plot_partials(par, amp).show()

# and plot dissonance curve with scale steps
plot_diss_curve(x, y, steps_12=True, steps_bp=False).show()

# Bohlen Pierce et, with odd harmonic partials

# fill in the frequencies of the partials here
par = [440*i for i in [1, 3, 5, 7, 9]]

# with the respective amplitudes
amp  = [0.8**i for i in range(len(par))]

# compute coordinates of curve
x, y = diss_curve(par, amp, ran=3.2, res=r)

# show partials as points, with harmonics as reference
plot_partials(par, amp).show()

# and plot dissonance curve with scale steps
plot_diss_curve(x, y, steps_12=False, steps_bp=True).show()

# Bohlen Pierce et, with non-harmonic partials

# fill in the frequencies of the partials here
par = [440*i for i in [1, 3, 3**(19/13), 3**(23/13), 9]]

# with the respective amplitudes
amp  = [0.8**i for i in range(len(par))]

# compute coordinates of curve
x, y = diss_curve(par, amp, ran=3.2, res=r)

# show partials as points, with harmonics as reference
plot_partials(par, amp).show()

# and plot dissonance curve with scale steps
plot_diss_curve(x, y, steps_12=False, steps_bp=True).show()