カオスでチャンスオペレーション：ダフィング方程式 編

#はじめに
カオス(力学系)からチャンスオペレーション(偶然性の音楽)ができないか？という試みをSiv3Dで行っています。カオスを用いることで、単純なランダムと違って曲を数式で一意に定義可能で、初期値鋭敏性により偶然性の音楽を生成することができ、有界性による限られた範囲での動きを音階に当てはめやすいというメリットが考えられます。

#ダフィング方程式での例
ここではダフィング方程式のサンプルをご紹介します。
https://twitter.com/hat_a2c/status/1308292986116005889

※リンク先の表題では「ダフィングアトラクタ」としておりましたが、ダフィング方程式のアトラクタは『ジャパニーズアトラクタ』が一般名称でした。

#使い方
ダフィング方程式は二重井戸ポテンシャルなので２状態の行き来を表現することができます。
そこで今回はx<0とx>0で平行調の関係になるように音階を配置しました。
dx/dtは音量としてHigh(音量大)、Low(音量小)の二値を配置しました。
方程式は⊿tを小さくして4次のルンゲクッタ法で求めています。

参考：
https://dn-voice.info/music-theory/godoken/
https://ja.wikipedia.org/wiki/Tonnetz
https://musicsounds.art/key_decision/
https://datt-music.com/ongaku-riron/key-insyo-tokutyo-matome/
http://www.u.tsukuba.ac.jp/~watanabe.norio.fw/diff2sol.htm

#ソースコード
OpenSiv3D(0.4.3)を前提としています。ここでは『ダフィングアトラクタ × Eマイナー近親調』のコードとしております。素人コードなので冗長な部分があるかと思いますが、ご活用頂ければと思います。

Main.cpp

# include <Siv3D.hpp>

void func(double n[10], double m)
{
    n[9] = n[8];
    n[8] = n[7];
    n[7] = n[6];
    n[6] = n[5];
    n[5] = n[4];
    n[4] = n[3];
    n[3] = n[2];
    n[2] = n[1];
    n[1] = n[0];
    n[0] = m;
}

LineString funcLine(double m[10], double n[10], Vec2 p, int d)
{
    Vec2 l0 = Vec2(m[9] * d + p.x, p.y - n[9] * d);
    Vec2 l1 = Vec2(m[8] * d + p.x, p.y - n[8] * d);
    Vec2 l2 = Vec2(m[7] * d + p.x, p.y - n[7] * d);
    Vec2 l3 = Vec2(m[6] * d + p.x, p.y - n[6] * d);
    Vec2 l4 = Vec2(m[5] * d + p.x, p.y - n[5] * d);
    Vec2 l5 = Vec2(m[4] * d + p.x, p.y - n[4] * d);
    Vec2 l6 = Vec2(m[3] * d + p.x, p.y - n[3] * d);
    Vec2 l7 = Vec2(m[2] * d + p.x, p.y - n[2] * d);
    Vec2 l8 = Vec2(m[1] * d + p.x, p.y - n[1] * d);
    Vec2 l9 = Vec2(m[0] * d + p.x, p.y - n[0] * d);
    return LineString
    {
        l0,l1,l2,l3,l4,l5,l6,l7,l8,l9
    };
}


void Main()
{
    //ウィンドウサイズ
    Window::Resize(450, 500);

    //発音データ（Eマイナー）
    const Audio pC(GMInstrument::Piano1, PianoKey::C3, 1.0s);
    const Audio pD(GMInstrument::Piano1, PianoKey::D3, 1.0s);
    const Audio pE(GMInstrument::Piano1, PianoKey::E3, 1.0s);
    const Audio pFs(GMInstrument::Piano1, PianoKey::FS3, 1.0s);
    const Audio pG(GMInstrument::Piano1, PianoKey::G3, 1.0s);
    const Audio pA(GMInstrument::Piano1, PianoKey::A3, 1.0s);
    const Audio pB(GMInstrument::Piano1, PianoKey::B3, 1.0s);

    //ダフィング方程式パラメータ
    double delta = 0.32;
    double gamma = 0.32;
    double omega = 0.60;

    //初期値
    double xo = 0.0;
    double yo = 0.0;
    double zo = 0.0;
    double x = xo;
    double y = yo;
    double z = zo;

    //閾値
    double left1 = -1.375;
    double left2 = -1.125;
    double left3 = -0.875;
    double left4 = -0.625;
    double left5 = -0.375;
    double left6 = -0.125;
    double right1 = 0.125;
    double right2 = 0.375;
    double right3 = 0.625;
    double right4 = 0.875;
    double right5 = 1.125;
    double right6 = 1.375;
    double upper = 0.25;
    double lower = -0.25;

    //表示用
    const Font font(20);
    const Font font2(15);
    double d = 100; //拡大倍率
    Vec2 px = Vec2(225, 260); //表示ポジション
    double lx[10] = { xo,xo,xo,xo,xo,xo,xo,xo,xo,xo };
    double ly[10] = { yo,yo,yo,yo,yo,yo,yo,yo,yo,yo };
    double lz[10] = { zo,zo,zo,zo,zo,zo,zo,zo,zo,zo };

    //処理用
    bool play = false; //再生用
    uint64 time = 0; //更新時間管理用
    uint64 timer = 0; //発音時間管理用
    int count = 0; //発音回数

    double speed = 10.5; //更新速度[倍速]
    int span = 150; //発音間隔[msec]
    
    while (System::Update())
    {
        if (KeySpace.down()) //スペースキーで起動
        {
            if (!play) //軌道
            {
                play = true;
                time = Time::GetMicrosecSinceEpoch();
                timer = Time::GetMicrosecSinceEpoch();
                count = 1;
            }
            else //リセット
            {
                xo = 0.0;
                yo = 0.0;
                zo = 0.0;
                x = xo;
                y = yo;
                z = zo;
                for (int i = 0; i < 10; i++)
                {
                    lx[i] = x;
                    ly[i] = y;
                    lz[i] = z;
                }
                play = false;
            }
        }

        if (play) //再生中処理
        {
            uint64 value = (Time::GetMicrosecSinceEpoch() - time); //更新時間間隔[usec]
            time = Time::GetMicrosecSinceEpoch();
            double h = (double)value * 0.001 * 0.001 * (speed / 30.0); //更新時間間隔[sec] ※30分割して繰返計算

            for (int i = 0; i < 30; i++) //ルンゲクッタ繰返
            {
                double k1x = h * y;
                double k1y = h * (- delta * y + x - Pow(x,3) + gamma * cos(omega * z));
                double sx = x + k1x / 2;
                double sy = y + k1y / 2;
                double sz = z + h / 2;
                double k2x = h * sy;
                double k2y = h * (- delta * sy + sx - Pow(sx,3) + gamma * cos(omega * sz));
                sx = x + k2x / 2;
                sy = y + k2y / 2;
                sz = z + h / 2;
                double k3x = h * sy;
                double k3y = h * (- delta * sy + sx - Pow(sx,3) + gamma * cos(omega * sz));
                sx = x + k3x;
                sy = y + k3y;
                sz = z + h;
                double k4x = h * sy;
                double k4y = h * (- delta * sy + sx - Pow(sx,3) + gamma * cos(omega * sz));
                x = x + (k1x + 2 * k2x + 2 * k3x * k4x) / 6;
                y = y + (k1y + 2 * k2y + 2 * k3y * k4y) / 6;
                z = z + h;
            }
            func(lx, x);
            func(ly, y);
            func(lz, z);

            value = (Time::GetMicrosecSinceEpoch() - timer); //再生時間[usec]
            uint64 next = span * 1000 * count; //発音時間間隔判定[usec]
            if (value >= next)
            {
                double volume = 0.5;
                if (y > upper || y < lower)
                {
                    volume = 0.5;
                }
                else
                {
                   volume = 0.15;
                }
                if (x < left1)
                {
                    pA.playOneShot(volume);
                }
                else if (x < left2)
                {
                    pC.playOneShot(volume);
                }
                else if (x < left3)
                {
                    pE.playOneShot(volume);
                }
                else if (x < left4)
                {
                    pG.playOneShot(volume);
                }
                else if (x < left5)
                {
                    pB.playOneShot(volume);
                }
                else if (x < left6)
                {
                    pD.playOneShot(volume);
                }
                else if (x < right1)
                {
                    pFs.playOneShot(volume);
                }
                else if (x < right2)
                {
                    pA.playOneShot(volume);
                }
                else if (x < right3)
                {
                    pC.playOneShot(volume);
                }
                else if (x < right4)
                {
                    pE.playOneShot(volume);
                }
                else if (x < right5)
                {
                    pG.playOneShot(volume);
                }
                else if (x < right6)
                {
                    pB.playOneShot(volume);
                }
                else
                {
                    pD.playOneShot(volume);
                }
                count += 1;
            }
        }

        if (1) //以下、表示用
        {
            LineString lineX = funcLine(lx, ly, px, d); //軌道表示
            lineX.draw(1, Palette::Skyblue);
            font2(U"A　C　E　G　B　D   F#  A　C　E　G　B　D").draw(px.x - 152, px.y + 105, Palette::White);
            Line(px.x - 160, px.y, px.x + 160, px.y).draw(1, Palette::White);;
            Line(px.x - 160, px.y - upper * d, px.x + 160, px.y - upper * d).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x - 160, px.y - lower * d, px.x + 160, px.y - lower * d).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x, px.y + 100, px.x, px.y - 100).draw(1, Palette::White);;
            Line(px.x + left1 * d, px.y + 100, px.x + left1 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + left2 * d, px.y + 100, px.x + left2 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + left3 * d, px.y + 100, px.x + left3 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + left4 * d, px.y + 100, px.x + left4 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + left5 * d, px.y + 100, px.x + left5 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + left6 * d, px.y + 100, px.x + left6 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + right1 * d, px.y + 100, px.x + right1 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + right2 * d, px.y + 100, px.x + right2 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + right3 * d, px.y + 100, px.x + right3 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + right4 * d, px.y + 100, px.x + right4 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + right5 * d, px.y + 100, px.x + right5 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            Line(px.x + right6 * d, px.y + 100, px.x + right6 * d, px.y - 100).draw(LineStyle::SquareDot, 1, Palette::White);
            font(U"x\"= -δx'+x-x +γcos(ωt)").draw(30, 10);
            font2(U"3").draw(150, 10);
            font(U"δ=",delta, U" , γ=",gamma, U", ω=",omega).draw(30, 35);
            font2(U"dx/dt").draw(px.x - 20, px.y - 125);
            font2(U"x").draw(px.x + 165, px.y - 12);
            font2(U"High").draw(px.x - 200, px.y - upper * d - 20 - 20);
            font2(U"Low").draw(px.x - 200, px.y - (upper + lower) / 2.0 * d - 15);
            font2(U"High").draw(px.x - 200, px.y - lower * d + 20);
        }
    }
}