1/fゆらぎの続き

前回書いた1/fゆらぎのコードを使ってハードウェア化しようと整数化をしていたところ、この間欠カオス法に面白い特徴があることが分かった。

まず、真っ当に単精度浮動小数で実装した間欠カオス法で、1024個の値を生成してパワースペクトルを調べたのがこちら。

FFTの窓関数を適用していないので、512以上の部分はイメージが出てスペクトルが強くなってたり、前回紹介した端っこの処理方法を変えているので低周波側でスペクトルが頭打ちになってたりしてるが、ざっくり1/fになっている。

これをそのまま24bit固定小数で実装しなおしたものがこちら。

floatの仮数部は24bit分なので、0.0～1.0区間では24bit固定小数と同じ精度になる。当然、ここではほとんど同じパターンになる（カオス法なので細かい違いは出る）。

次に16bit固定小数にしたものがこちら。

1/fのパワースペクトルにならず、1/fっぽいグランドノイズと特定の周波数で強いピークの2つに分かれている。

さらに8bit固定小数にする。

ピークそのもののスペクトルは1/fっぽく並んでいるものの、特定の周波数（ここでは49の倍数になっている）のみの成分だけになり、ほぼ周期性のあるパターンになってしまっている。
こういうパターンでたとえばロウソクの炎の模擬をすると、安いキャンドルICのような周期的な点滅が目立つ挙動になってしまう。

そこで、どれぐらいの精度があれば1/fゆらぎが得られるかを調べたところ、18bit固定小数でピークが目立たないパターンになった。

18bit長はハードウェア実装する場合に18bit×18bitの乗算器1個でカバーできるため都合が良い。実際これがギリギリ下限の精度だろう。MCUでソフト実装する場合は18bitに削るメリットは無いので、24bit固定小数を使う方がいいだろう。

●結論
間欠カオス法には18bit以上の演算精度が必要である。

おまけ。
MCU向けに擬似ランダムをXorshift32に置き換えた、24bit固定小数間欠カオス法のソースコード

#include <stdint.h>

uint32_t xorshift32(void)
{
  static uint32_t x = 2463534242;

  x ^= (x << 13);
  x ^= (x >> 17);
  x ^= (x << 15);
  return x;
}

uint32_t pow_u24(uint32_t a, int s)
{
  uint64_t m, ans;

  m = a & 0xffffff;
  ans = (m * m) >> s;
  return (uint32_t)ans;
}

uint32_t fluctuator(void)
{
  static uint32_t f = 0;

  if ((f & 0xf00000) == 0 || (f & 0xf00000) == 0xf00000) {
    f = (xorshift32() & 0x7fffff) + 0x400000;
  } else {
    if (f & (1UL<<23)) {
      f = (f - pow_u24((1UL<<24) - f, 23)) & 0xffffff;
    } else {
      f = (f + pow_u24(f, 23)) & 0xffffff;
    }
  }
  return f;
}