カオスでチャンスオペレーション：ローレンツ方程式 編

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

#ローレンツ方程式での例
ここではローレンツ方程式のサンプルをご紹介します。
https://twitter.com/hat_a2c/status/1308263997246898176

#使い方
ローレンツ方程式は3変数ですので、今回はx軸に音階、z軸に音量を割り当てています。
音階はコード理論を参考に並べるとイイ感じになります。
音量はHigh(音量大)、Low(音量小)の二値としています。
方程式は⊿tを小さくして4次のルンゲクッタ法で求めています。

参考：
https://dn-voice.info/music-theory/godoken/
https://ja.wikipedia.org/wiki/Tonnetz
https://musicsounds.art/key_decision/

#ソースコード
OpenSiv3D(0.4.3)を前提としています。素人コードなので冗長な部分があるかと思いますが、ご活用頂ければと思います。

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);
    
    //発音データ
    constexpr GMInstrument instrument = GMInstrument::Piano1; //楽器
    const Audio pF(instrument, PianoKey::F3, 1.2s);
    const Audio pA(instrument, PianoKey::A3, 1.2s);
    const Audio pC(instrument, PianoKey::C3, 1.2s);
    const Audio pE(instrument, PianoKey::E3, 1.2s);
    const Audio pG(instrument, PianoKey::G3, 1.2s);
    const Audio pB(instrument, PianoKey::B3, 1.2s);
    const Audio pD(instrument, PianoKey::D3, 1.2s);

    //ローレンツ方程式パラメータ
    double sigma = 10.0;
    double b = 12.0 / 5.0;
    double r = 25.0;

    //初期値
    double xo = Random() * 20 - 10;
    double yo = Random() * 20 - 10;
    double zo = Random() * 40;
    double x = xo;
    double y = yo;
    double z = zo;

    //閾値
    int left1 = -10;
    int left2 = -6;
    int left3 = -2;
    int right1 = 2;
    int right2 = 6;
    int right3 = 10;
    int upper = 25;
    int lower = 15;

    //表示用
    const Font font(20);
    const Font font2(15);
    int d = 5; //拡大倍率
    Vec2 px = Vec2(225, 360); //表示ポジション
    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 on = 0; //発音判定用
    int count = 0; //発音回数
    double speed = 3; //更新速度[倍速]
    int span = 180; //発音間隔[msec]
    double change = 64.0; //平行調切替用

    while (System::Update())
    {
        if (KeySpace.down()) //スペースキーで起動
        {
            if (!play) //起動
            {
                play = true;
                time = Time::GetMicrosecSinceEpoch();
                timer = Time::GetMicrosecSinceEpoch();
                count = 1;
                on = 0;
            }
            else //リセット
            {
                x = Random() * 20 - 10;
                y = Random() * 20 - 10;
                z = Random() * 40;
                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 * 0.1; //更新時間間隔[sec] ※10分割して繰返計算

            for (int i = 0; i < 10 * speed; i++) //ルンゲクッタ繰返
            {
                double k1x = h * (sigma * (y - x));
                double k1y = h * (r * x - y - x * z);
                double k1z = h * (x * y - b * z);
                double sx = x + k1x / 2;
                double sy = y + k1y / 2;
                double sz = z + k1z / 2;
                double k2x = h * (sigma * (sy - sx));
                double k2y = h * (r * sx - sy - sx * sz);
                double k2z = h * (sx * sy - b * sz);
                sx = x + k2x / 2;
                sy = y + k2y / 2;
                sz = z + k2z / 2;
                double k3x = h * (sigma * (sy - sx));
                double k3y = h * (r * sx - sy - sx * sz);
                double k3z = h * (sx * sy - b * sz);
                sx = x + k3x;
                sy = y + k3y;
                sz = z + k3z;
                double k4x = h * (sigma * (sy - sx));
                double k4y = h * (r * sx - sy - sx * sz);
                double k4z = h * (sx * sy - b * sz);
                x = x + (k1x + 2 * k2x + 2 * k3x * k4x) / 6;
                y = y + (k1y + 2 * k2y + 2 * k3y * k4y) / 6;
                z = z + (k1z + 2 * k2z + 2 * k3z * k4z) / 6;
            }
            func(lx, x); 
            func(ly, y);
            func(lz, z);

            value = (Time::GetMicrosecSinceEpoch() - timer); //再生時間値[usec]
            uint64 next = span * 1000 * count; //発音時間間隔判定[usec]
            if (on == 0)
            {
                if ((int)((double)count / change) % 2 == 0) //平行調変換判定
                {
                    on = 1; //メジャーコード
                }
                else
                {
                    on = 2; //マイナーコード
                }
            }
            else if (on == 1 && value >= next) //メジャーコード
            {
                double volume = 0.5;
                if (z > upper || z < lower)
                {
                    volume = 0.5;
                }
                else
                {
                   volume = 0.15;
                }
                if (x < left1)
                {
                    pD.playOneShot(volume);
                }
                else if (x < left2)
                {
                    pB.playOneShot(volume);
                }
                else if (x < left3)
                {
                    pG.playOneShot(volume);
                }
                else if (x < right1)
                {
                    pE.playOneShot(volume);
                }
                else if (x < right2)
                {
                    pC.playOneShot(volume);
                }
                else if (x < right3)
                {
                    pA.playOneShot(volume);
                }
                else
                {
                    pF.playOneShot(volume);
                }
                count += 1;
                on = 0; //発音OFF
            }
            else if (on == 2 && value >= next) //マイナーコード
            {
                double volume = 0.5;
                if (z > upper || z < lower)
                {
                    volume = 0.5;
                }
                else
                {
                    volume = 0.15;
                }
                if (x < left1)
                {
                    pB.playOneShot(volume);
                }
                else if (x < left2)
                {
                    pG.playOneShot(volume);
                }
                else if (x < left3)
                {
                    pE.playOneShot(volume);
                }
                else if (x < right1)
                {
                    pC.playOneShot(volume);
                }
                else if (x < right2)
                {
                    pA.playOneShot(volume);
                }
                else if (x < right3)
                {
                    pF.playOneShot(volume);
                }
                else
                {
                    pD.playOneShot(volume);
                }
                count += 1;
                on = 0; //発音OFF
            }
        }


        LineString lineHist = funcLine(lx, lz, px, d); //軌道表示
        if ((int)((double)count / change) % 2 == 0) //平行調変換判定
        {
            lineHist.draw(1, Palette::Pink);
            font2(U"D   B   G   E   C   A   F").draw(px.x - 70, px.y + 5, Palette::Pink);            
        }
        else
        {
            lineHist.draw(1, Palette::Skyblue);
            font2(U"B   G   E   C   A   F   D").draw(px.x - 70, px.y + 5, Palette::Skyblue);
        }

        Line(px.x - 100, px.y, px.x + 100, px.y).draw(1, Palette::White);;
        Line(px.x, px.y, px.x, px.y - 200).draw(1, Palette::White);;
        Line(px.x + left1 * d, px.y, px.x + left1 * d, px.y - 200).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x + left2 * d, px.y, px.x + left2 * d, px.y - 200).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x + left3 * d, px.y, px.x + left3 * d, px.y - 200).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x + right1 * d, px.y, px.x + right1 * d, px.y - 200).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x + right2 * d, px.y, px.x + right2 * d, px.y - 200).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x + right3 * d, px.y, px.x + right3 * d, px.y - 200).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x - 100, px.y - upper * d, px.x + 100, px.y - upper * d).draw(LineStyle::SquareDot, 1, Palette::White);
        Line(px.x - 100, px.y - lower * d, px.x + 100, px.y - lower * d).draw(LineStyle::SquareDot, 1, Palette::White);

        font(U"x'= σ(y-x)").draw(30, 10);
        font(U"y'= rx - y - xz").draw(30, 35);
        font(U"z'= xy - bz").draw(30, 60);
        font(U"σ=", sigma, U", r=", r, U", b=", b, U", ").draw(30,85);
        font(U"(xo,yo,zo)=(", xo, U",", yo, U",", zo, U")").draw(30, 110);

        font2(U"z").draw(px.x - 4, px.y - 225);
        font2(U"x").draw(px.x + 104, px.y-12);
        font2(U"High").draw(px.x - 140, px.y - upper * d - 20 - 20);
        font2(U"Low").draw(px.x - 140, px.y - (upper + lower) / 2.0 * d - 15);
        font2(U"High").draw(px.x - 140, px.y - lower * d + 20);

    }
}