Dissonance calculation for up to 19 harmonics with microtonic scales shown -- also see https://gist.github.com/endolith/3066664

_image.md

calculate.py

"""
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")