JS80P Waveshapers and Distortions¶

import math as m
import matplotlib.pyplot as plt


def plot_wavefolder():
    FOLD_MIN = 0.0
    FOLD_TRANSITION = 0.5
    FOLD_MAX = 5.0 + FOLD_TRANSITION

    PI = m.pi
    PI_DOUBLE = PI * 2.0
    PI_HALF = PI / 2.0
    PI_SQR = PI * PI

    TRIANGLE_SCALE = 8.0 / PI_SQR

    S0 = TRIANGLE_SCALE
    S1 = PI_HALF
    S2 = TRIANGLE_SCALE / 9.0
    S3 = PI_HALF * 3.0
    S4 = TRIANGLE_SCALE / 25.0
    S5 = PI_HALF * 5.0
    S6 = TRIANGLE_SCALE * 2.0 / PI
    S7 = TRIANGLE_SCALE / (27.0 * PI)
    S8 = TRIANGLE_SCALE / (125.0 * PI)

    TRIG_OFFSET = PI_DOUBLE * m.ceil(FOLD_MAX * S5)

    def f(x):
        return (
            S0 * m.sin(S1 * x + TRIG_OFFSET)
            - S2 * m.sin(S3 * x + TRIG_OFFSET)
            + S4 * m.sin(S5 * x + TRIG_OFFSET)
        )

    def F0(x):
        return (
            -S6 * m.cos(S1 * x + TRIG_OFFSET)
            + S7 * m.cos(S3 * x + TRIG_OFFSET)
            - S8 * m.cos(S5 * x + TRIG_OFFSET)
        )

    N = 2000
    folding = 5.0
    prev = -folding
    F0_prev = F0(prev)
    approximation = (F0(-folding + 1.0 / N) - F0(-folding)) * N

    xs = []
    fs = []
    approximations = []

    for i in range(N):
        x = (2.0 * (i / N)) - 1.0
        x = folding * x
        delta = x - prev

        xs.append(x)
        fs.append(f(x))

        if abs(delta) > 0.0:
            F0_x = F0(x)
            approximation = (F0_x - F0_prev) / delta
            approximations.append(approximation)
            F0_prev = F0_x
        else:
            approximations.append(approximation)

        prev = x

    plt.figure(figsize=(15, 6))
    plt.plot(xs, fs, label="true f(x)")
    plt.plot(xs, approximations, label="AAD approx.")
    plt.legend()
    plt.title("Wavefolder")
    plt.show()


plot_wavefolder()

import math as m
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt


INPUT_MAX = 3.0
INPUT_MIN = - INPUT_MAX

PI2 = m.pi * 2.0


def tanh_f(x, steepness):
    return m.tanh(steepness * x * 0.5)


def tanh_F0(x, steepness, c):
    return x + (2.0 / steepness) * m.log1p(m.exp(-steepness * x)) + c


def join(alpha, gamma, h, dh, H):
    def combine_branches(funcs):
        return lambda x: (
            funcs[0](x) if x < gamma
                else (funcs[1](x) if x < 1.0 else funcs[2](x))
        )

    E = (alpha - 1.0) / 4.0
    cg = -9.0 * E

    g = lambda x: E * (x - 3.0) ** 2.0 + 1.0
    dg = lambda x: 2.0 * E * x - 6.0 * E
    G = lambda x: (((E / 3.0) * x - 3.0 * E) * x + 9.0 * E + 1.0) * x + cg

    L1 = gamma
    L2 = gamma ** 2.0
    L3 = gamma ** 3.0
    L4 = gamma ** 4.0

    eq_mtx = np.array(
        [
            [   3.0 * L2,   2.0 * L1,        1.0,         0.0,   0.0,     0.0],
            [         L3,         L2,         L1,         1.0,   0.0,     0.0],
            [   3.0 * L4,   4.0 * L3,   6.0 * L2,   12.0 * L1,  12.0,   -12.0],
            [        3.0,        2.0,        1.0,         0.0,   0.0,     0.0],
            [        1.0,        1.0,        1.0,         1.0,   0.0,     0.0],
            [        3.0,        4.0,        6.0,        12.0,  12.0,     0.0],
        ]
    )
    eq_v = np.array(
        [
                  dh(gamma),
                   h(gamma),
            12.0 * H(gamma),
                1.0 - alpha,
                      alpha,
                12.0 * G(1),
        ]
    )
    A, B, C, D, cf, ch = la.inv(eq_mtx).dot(eq_v)

    H_p_ch = lambda x: H(x) + ch

    f = lambda x: ((A * x + B) * x + C) * x + D
    df = lambda x: (3.0 * A * x + 2.0 * B) * x + C
    F = lambda x: ((((A / 4.0) * x + B / 3.0) * x + C / 2.0) * x + D) * x + cf

    return (
        (combine_branches((h, f, g)), combine_branches((H_p_ch, F, G))),
        (A, B, C, D, cf, ch)
    )


def bit_crush(type_const, k, alpha, gamma):
    PI2k = PI2 * k
    W = 1.0 / (1.4 * m.pi * k)

    h_k = lambda x: x - W * m.sin(PI2k * x)
    dh_k = lambda x: 1 - W * m.cos(PI2k * x)
    H_k = lambda x: x ** 2.0 / 2.0 + W * m.cos(PI2k * x) / PI2k

    funcs, (A, B, C, D, cf, ch) = join(alpha, gamma, h_k, dh_k, H_k)

    print(f"""\

    initialize_bit_crush_tables(
        {type_const},
        {k},
        {alpha},
        {gamma},
        {A},
        {B},
        {C},
        {D},
        {cf},
        {ch}
    );""")

    return alpha, gamma, funcs


def harmonics(type_const, alpha, gamma, w1, w3, w5, norm=1.0):
    w1 *= norm
    w3 *= norm
    w5 *= norm

    # Chebyshev polynomials and their derivatives and integrals
    t1 = lambda x: x
    t3 = lambda x: 4.0 * x ** 3.0 - 3.0 * x
    t5 = lambda x: 16.0 * x ** 5.0 - 20.0 * x ** 3.0 + 5.0 * x

    dt1 = lambda x: 1.0
    dt3 = lambda x: 12.0 * x ** 2.0 - 3.0
    dt5 = lambda x: 80.0 * x ** 4.0 - 60.0 * x ** 2.0 + 5.0

    T1 = lambda x: (1.0 / 2.0) * x ** 2.0
    T3 = lambda x: x ** 4.0 - (3.0 / 2.0) * x ** 2.0
    T5 = lambda x: (16.0 / 6.0) * x ** 6.0 - 5.0 * x ** 4.0 + (5.0 / 2.0) * x ** 2.0

    h = lambda x: w1 * t1(x) + w3 * t3(x) + w5 * t5(x)
    dh = lambda x: w1 * dt1(x) + w3 * dt3(x) + w5 * dt5(x)
    H = lambda x: w1 * T1(x) + w3 * T3(x) + w5 * T5(x)

    funcs, (A, B, C, D, cf, ch) = join(alpha, gamma, h, dh, H)

    print(f"""\

    initialize_harmonic_tables(
        {type_const},
        {w1},
        {w3},
        {w5},
        {alpha},
        {gamma},
        {A},
        {B},
        {C},
        {D},
        {cf},
        {ch}
    );""")

    return alpha, gamma, funcs


def expression(type_const, alpha, gamma, h_expr, H_expr, dh_expr):
    def cpp_to_py(expr):
        return (
            expr
                .replace("std::", "m.")
                .replace("Math::PI_HALF", "(m.pi / 2.0)")
                .replace("Math::PI", "m.pi")
        )

    h = eval("lambda x: " + cpp_to_py(h_expr))
    dh = eval("lambda x: " + cpp_to_py(dh_expr))
    H = eval("lambda x: " + cpp_to_py(H_expr))
    
    funcs, (A, B, C, D, cf, ch) = join(alpha, gamma, h, dh, H)

    print(f"""\

    initialize_spline_tables(
        {type_const},
        [&](Number const x, Number& y, Number& antiderivative) {{
            y = {h_expr};
            antiderivative = {ch} + {H_expr};
        }},
        {alpha},
        {gamma},
        {A},
        {B},
        {C},
        {D},
        {cf}
    );""")

    return alpha, gamma, funcs


SQR = (
    4.0 / m.pi, 4.0 / (3.0 * m.pi),
    4.0 / (5.0 * m.pi),
)
TRI = (
    8.0 / (1.0 * m.pi) ** 2.0,
    -8.0 / (3.0 * m.pi) ** 2.0,
    8.0 / (5.0 * m.pi) ** 2.0,
)

TANH_F0_3_C = INPUT_MAX - tanh_F0(INPUT_MAX, 3.0, 0.0)
TANH_F0_5_C = INPUT_MAX - tanh_F0(INPUT_MAX, 5.0, 0.0)
TANH_F0_10_C = INPUT_MAX - tanh_F0(INPUT_MAX, 10.0, 0.0)

DISTORTIONS = (
    ("tanh(3x/2)",  (1.0, 1.0, (lambda x: tanh_f(x,  3.0), lambda x: tanh_F0(x,  3.0,  TANH_F0_3_C)))),
    ("tanh(5x/2)",  (1.0, 1.0, (lambda x: tanh_f(x,  5.0), lambda x: tanh_F0(x,  5.0,  TANH_F0_5_C)))),
    ("tanh(10x/2)", (1.0, 1.0, (lambda x: tanh_f(x, 10.0), lambda x: tanh_F0(x, 10.0, TANH_F0_10_C)))),

    (
        "sin(x*pi/2)",
        expression(
            "TYPE_SIN",
            0.99,
            0.9,
            h_expr="std::sin(x * Math::PI_HALF)",
            H_expr="(-2.0 / Math::PI) * std::cos(x * Math::PI_HALF)",
            dh_expr="Math::PI_HALF * std::cos(x * Math::PI_HALF)",
        ),
    ),

    (
        "x^2",
        expression(
            "TYPE_SQR",
            0.95,
            0.95,
            h_expr="x * x",
            H_expr="x * x * x / 3.0",
            dh_expr="2.0 * x",
        ),
    ),

    (
        "sqrt(x)",
        expression(
            "TYPE_SQRT",
            0.97,
            0.9,
            h_expr="std::sqrt(x)",
            H_expr="(2.0 / 3.0) * std::pow(x, 3.0 / 2.0)",
            dh_expr="(1.0 / 2.0) * std::pow(x + 0.000001, -1.0 / 2.0)",
        ),
    ),

    (
        "x^3",
        expression(
            "TYPE_CUBE",
            0.95,
            0.97,
            h_expr="x * x * x",
            H_expr="x * x * x * x / 4.0",
            dh_expr="3.0 * x * x",
        ),
    ),

    (
        "cbrt(x)",
        expression(
            "TYPE_CBRT",
            0.97,
            0.9,
            h_expr="std::pow(x, 1.0 / 3.0)",
            H_expr="(3.0 / 4.0) * std::pow(x, 4.0 / 3.0)",
            dh_expr="(1.0 / 3.0) * std::pow(x + 0.000001, -2.0 / 3.0)",
        ),
    ),

    ("harmonic 1+3", harmonics("TYPE_HARMONIC_13",    0.90, 0.80, 1.0, 1.0, 0.0, norm=1.10)),
    ("harmonic 1+5", harmonics("TYPE_HARMONIC_15",    0.90, 0.80, 1.0, 0.0, 1.0, norm=0.67)),
    ("harmonic 1+3+5", harmonics("TYPE_HARMONIC_135", 0.90, 0.80, 1.0, 1.0, 1.0, norm=1.20)),
    ("harmonic SQR", harmonics("TYPE_HARMONIC_SQR",   0.90, 0.80, *SQR, norm=1.10)),
    ("harmonic TRI", harmonics("TYPE_HARMONIC_TRI",   0.90, 0.80, *TRI, norm=1.20)),

    ("bit crush 1", bit_crush("TYPE_BIT_CRUSH_1",       1, 0.980, 0.720)),
    ("bit crush 2", bit_crush("TYPE_BIT_CRUSH_2",       2, 0.980, 0.630)),
    ("bit crush 3", bit_crush("TYPE_BIT_CRUSH_3",       4, 0.950, 0.810)),
    ("bit crush 4", bit_crush("TYPE_BIT_CRUSH_4",       8, 0.970, 0.903)),
    ("bit crush 4.6", bit_crush("TYPE_BIT_CRUSH_4_6",  12, 0.975, 0.936)),
    ("bit crush 5", bit_crush("TYPE_BIT_CRUSH_5",      16, 0.980, 0.953)),
    ("bit crush 5.6", bit_crush("TYPE_BIT_CRUSH_5_6",  24, 0.986, 0.968)),
    ("bit crush 6", bit_crush("TYPE_BIT_CRUSH_6",      32, 0.990, 0.976)),
    ("bit crush 6.6", bit_crush("TYPE_BIT_CRUSH_6_6",  48, 0.990, 0.984)),
    ("bit crush 7", bit_crush("TYPE_BIT_CRUSH_7",      64, 0.993, 0.988)),
    ("bit crush 7.6", bit_crush("TYPE_BIT_CRUSH_7_6",  96, 0.997, 0.992)),
    ("bit crush 8", bit_crush("TYPE_BIT_CRUSH_8",     128, 0.999, 0.994)),
    ("bit crush 8.6", bit_crush("TYPE_BIT_CRUSH_8_6", 192, 0.999, 0.993)),
    ("bit crush 9", bit_crush("TYPE_BIT_CRUSH_9",     256, 0.999, 0.995)),
)


def plot_distortion(dist_type):
    if not (0 <= dist_type < len(DISTORTIONS)):
        raise ValueError(
            f"Distortion type must be between 0 and {len(DISTORTIONS) - 1}, got: {dist_type}",
        )

    title, (alpha, gamma, (f_pos, F0_pos)) = DISTORTIONS[dist_type]

    f = lambda x: f_pos(x) if x >= 0.0 else -f_pos(-x)
    F0 = lambda x: F0_pos(x) if x >= 0.0 else F0_pos(-x)

    begin = -3.5
    end = 3.5
    range_ = end - begin
    N = 2000
    prev = begin
    F0_prev = F0(prev)

    xs = []
    fs = []
    F0s = []
    guide = []
    approximations = []
    F0_x = 0.0
    f_x = 0.0
    M = 0.0
    gM = 0.0

    for i in range(N):
        x = begin + range_ * (i / N)
        delta = x - prev

        if x < INPUT_MIN:
            F0_x = -x
            f_x = -1.0
        elif x >= INPUT_MAX:
            F0_x = x
            f_x = 1.0
        else:
            F0_x = F0(x)
            f_x = f(x)

        abs_f_x = abs(f_x)

        if abs_f_x > gM:
            gM = abs_f_x

        if x < -1.0:
            guide.append(-1.0)
        elif x > 1.0:
            guide.append(1.0)
        else:
            if 0 <= x <= gamma and abs_f_x > M:
                M = abs_f_x

            guide.append(x)

        xs.append(x)
        fs.append(f_x)
        F0s.append(F0_x)

        if abs(x - INPUT_MIN) < 0.0000001:
            approximations.append(-1.0)
        elif abs(x - INPUT_MAX) < 0.0000001:
            approximations.append(1.0)
        elif abs(delta) > 0.00000001:
            approximations.append((F0_x - F0_prev) / delta)
        else:
            approximations.append(f_x)

        prev = x
        F0_prev = F0_x

    print(f"global_max={gM}")
    print(f"max={M}")
    print(f"norm={1.0 / M}")
    print(f"norm_alpha={alpha / M}")

    plt.figure(figsize=(15, 7.5))
    plt.plot(xs, guide, label="guide")
    plt.plot(xs, fs, label="true f(x)")
    plt.plot(xs, approximations, label="AAD approx.")
    plt.plot(xs, F0s, label="F0(x)")
    plt.title(f"Distortion {title}")
    plt.legend()
    plt.show()

plot_distortion(0)

plot_distortion(1)

plot_distortion(2)

plot_distortion(3)

plot_distortion(4)

plot_distortion(5)

plot_distortion(6)

plot_distortion(7)

plot_distortion(8)

plot_distortion(9)

plot_distortion(10)

plot_distortion(11)

plot_distortion(12)

plot_distortion(13)

plot_distortion(14)

plot_distortion(15)

plot_distortion(16)

plot_distortion(17)

plot_distortion(18)

plot_distortion(19)

plot_distortion(20)

plot_distortion(21)

plot_distortion(22)

plot_distortion(23)

plot_distortion(24)

plot_distortion(25)

plot_distortion(26)

import math as m
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt


# The spline will be linear on the [0, gamma] interval.
gamma = 1.0 / 16.0

# Parameter: f'(gamma) := beta
beta = 5.0 / 8.0

# Parameter: f(1) = alpha
alpha = 899.0 / 1024.0


alpha_m_1_o_4 = (alpha - 1.0) / 4.0


# These will be set by solve_spline()
A = B = C = D = cf = ch = 0.0

cg = -9.0 * alpha_m_1_o_4


h = lambda x: beta * x
H = lambda x: (1.0 / 2.0) * beta * x ** 2.0 + ch


f = lambda x: ((A * x + B) * x + C) * x + D
F = lambda x: ((((A / 4.0) * x + (B / 3.0)) * x + (C / 2.0)) * x + D) * x + cf


g = lambda x: (
    (alpha_m_1_o_4 * x - 6.0 * alpha_m_1_o_4) * x + 9.0 * alpha_m_1_o_4 + 1.0
)
G = lambda x: (
    (
        ((alpha_m_1_o_4 / 3.0) * x - 3.0 * alpha_m_1_o_4) * x
        + 9.0 * alpha_m_1_o_4 + 1.0
    ) * x + cg
)


INPUT_MAX = 3.0


def solve_spline():
    global A, B, C, D, cf, ch

    A1 = 3.0 * gamma ** 4.0
    B1 = 4.0 * gamma ** 3.0
    C1 = 6.0 * gamma ** 2.0
    D1 = 12.0 * gamma

    A5 = 3.0 * gamma ** 2.0
    B5 = 2.0 * gamma

    A6 = gamma ** 3.0 - 1.0
    B6 = gamma ** 2.0 - 1.0
    C6 = gamma - 1.0

    eq_mtx = np.array(
        [
            [-12.0,  12.0,    D1,    A1,    B1,    C1],
            [  0.0,  12.0,  12.0,   3.0,   4.0,   6.0],
            [  0.0,   0.0,   1.0,   1.0,   1.0,   1.0],
            [  0.0,   0.0,   0.0,   3.0,   2.0,   1.0],
            [  0.0,   0.0,   0.0,    A5,    B5,   1.0],
            [  0.0,   0.0,   0.0,    A6,    B6,    C6],
        ]
    )
    eq_v = np.array(
        [
            6.0 * beta * gamma ** 2.0,
            20.0 - 8.0 * alpha,
            alpha,
            1.0 - alpha,
            beta,
            beta * gamma - alpha,
        ]
    )
    ch, cf, D, A, B, C = la.inv(eq_mtx).dot(eq_v)


def spline(x):
    sgn = 1.0

    if x < 0.0:
        sgn = -1.0
        x = -x

    if x >= INPUT_MAX:
        return sgn

    if x >= 1.0:
        return sgn * g(x)

    if x >= gamma:
        return sgn * f(x)

    return sgn * h(x)


def spline_int(x):
    if x < 0.0:
        x = -x

    if x >= INPUT_MAX:
        return x

    if x >= 1.0:
        return G(x)

    if x >= gamma:
        return F(x)

    return H(x)


def plot_reduce_dist():
    N = 2000
    width = 3.5

    solve_spline()

    funcs = (
        ("0", lambda x: 0, (([], 1),)),
        ("x", lambda x: x, (([], 1),)),
        ("int S(x)", lambda x: spline_int(x), (([], 1),)),
        (
            "S(x) +{} iter(s)",
            lambda x: spline(x),
            (([], 1), ([], 5), ([], 35), ([], 90), ([], 10), ([], 5), ([], 3)),
        ),
    )

    xs = []

    min_d = 999.0

    for i in range(N):
        x = (2.0 * (i / N)) - 1.0
        x = width * x

        if x > 0.0:
            s = spline(x)
            d = x - s

            if d < min_d:
                min_d = d

        xs.append(x)

        for label_tpl, f, iters in funcs:
            fx = x

            for values, n in iters:
                for i in range(n):
                    fx = f(fx)

                values.append(fx)

    values = [
        (label_tpl.format(n), v)
        for label_tpl, f, iters in funcs
        for v, n in iters
    ]

    print(f"""
min_d = {min_d}

alpha = {alpha}
beta = {beta}
gamma = {gamma}

cg = {cg}

constexpr Number A = {A};
constexpr Number B = {B};
constexpr Number C = {C};
constexpr Number D = {D};
constexpr Number cf = {cf};
constexpr Number ch = {ch};

""")

    plt.figure(figsize=(15, 15))

    for label, v in values:
        plt.plot(xs, v, label=label)

    plt.title("Reduce")
    plt.legend()
    plt.show()


plot_reduce_dist()

import math as m
import matplotlib.pyplot as plt


def smooth_smooth_steeper(x):
    a = -252.0
    b = 1386.0
    c = 3080.0
    d = 3465.0
    e = 1980.0
    f = 462.0

    return ((((((((((a * x + b) * x - c) * x + d) * x - e) * x + f) * x) * x) * x) * x) * x) * x


def sharp_sharp_steeper(x):
    a = 1024.0
    b = 5632.0
    c = 14080.0
    d = 21120.0
    e = 21120.0
    f = 14784.0
    g = 7392.0
    h = 2640.0
    i = 660.0
    j = 110.0
    k = 11.0

    return ((((((((((a * x - b) * x + c) * x - d) * x + e) * x - f) * x + g) * x - h) * x + i) * x - j) * x + k) * x


def plot_curves():
    curves = (
        ("smooth-smooth",         lambda x: ((-2.0 * x + 3.0) * x) * x),
        ("smooth-smooth steep",   lambda x: ((((6.0 * x - 15.0) * x + 10.0) * x) * x) * x),
        ("smooth-smooth steeper", smooth_smooth_steeper),
        ("sharp-sharp",           lambda x: ((4.0 * x - 6.0) * x + 3.0) * x),
        ("sharp-sharp steep",     lambda x: ((((16.0 * x - 40.0) * x + 40.0) * x - 20.0) * x + 5.0) * x),
        ("sharp-sharp steeper",   sharp_sharp_steeper),
        ("smooth-sharp",          lambda x: x ** 2.0),
        ("smooth-sharp steep",    lambda x: x ** 3.0),
        ("smooth-sharp steeper",  lambda x: x ** 5.0),
        ("sharp-smooth",          lambda x: x * (1.0 - m.log(x + 0.001)) / (1.0 - m.log(1.001))),
        ("sharp-smooth steep",    lambda x: (x * (1.0 - m.log(x + 0.001)) / (1.0 - m.log(1.001))) ** (2.0 / 3.0)),
        ("sharp-smooth steeper",  lambda x: (x * (1.0 - m.log(x + 0.001)) / (1.0 - m.log(1.001))) ** (1.0 / 3.0)),
    )

    N = 2000
    xs = [i / N for i in range(N + 1)]

    plt.figure(figsize=(15, 15))
    
    for name, f in curves:
        fxs = [f(x) for x in xs]
        plt.plot(xs, fxs, label=name)
    
    plt.title("Curves")
    plt.legend()
    plt.show()


plot_curves()