はじめに
炎のゆらぎを始めとして，自然界に多く存在するといわれる1/fゆらぎをデジタル上で再現します．今回は実際にLEDを光らせて，1/fゆらぎを表現します．
1/fゆらぎとは
1/fゆらぎの正体とは，-10dB/decつまり周波数がパワースペクトルを持つノイズのことです．ピンクノイズとも呼ばれたりします．
一方ホワイトノイズ(一様乱数)は，周波数に依らず一様のパワーを持っています．
1/fゆらぎの生成方法
ホワイトノイズを生成するライブラリは多く存在するかと思いますが，ピンクノイズを生成してくれるライブラリは多くありません．
たとえば組み込み用途で考えると，ホワイトノイズの生成に，Cの場合はstdlib.hのrand()やC++の場合はstd::randomなどの標準ライブラリが利用できます．標準ライブラリが使用できない場合でも，ホワイトノイズをせ生成できることは多いかと思います．
ですが，ピンクノイズの生成したいとなると，自前で用意する必要が出てきます．
ホワイトノイズからフィルタによる生成
今回は，ホワイトノイズを基にピンクノイズを生成します．
やることとしては簡単で，ホワイトノイズを-10dB/decの特性を持つフィルタに通してあげることでピンクノイズを生成します．
-10dB/decの特性を持つフィルタ
位相遅れフィルタ
$$
P(s)=\dfrac{T_1s+1}{T_2s+1}
$$
を直列に繋ぐごとによって，
$$
P(s)=\prod_{k}^n\dfrac{4^ks+1}{2\cdot4^ks+1}
$$
-10dB/decの特性を持つフィルタを作成します．
このとき，$T1$と$T2$の間は1oct，各フィルタ間は2octごとに並べます．
$-4\leq k \leq 1$とした際の特性が以下のようになります．
$k$の値は，特性を得たい周波数帯域でフィルタが有効になるように設定します．
この際，無駄に帯域を増やすと，フィルタ次数とそれに伴い計算量も増えることに注意が必要です．今回は，目で見て気になりそうだなと思われる周波数帯域0.1~100Hzあたりを対象としています，
このフィルタにホワイトノイズを通した出力が以下のようになります．
当然ですが，出力の周波数成分に着目するとフィルタ同様の-10dB/decのスペクトル特性を持っています．
以下は，上記のフィルタ作成に使用したPythonスクリプトになります．
import matplotlib.pyplot as plt
import numpy as np
from control import tf,c2d
from control.matlab import lsim, bode, tfdata
# サンプリング時間
Ts = 0.005
# フィルタ作成
P = [tf([1],[1])]
for k in range(-4, 2):
tmp = tf([4**k, 1], [2*4**k, 1])
P.append(tmp)
P[0] = P[0]*tmp
gain, phase, w = bode(P, np.logspace(-3, 3)*(2*np.pi), Hz=True)
plt.grid()
plt.savefig("Filter_Bode.png")
# 入力：ホワイトノイズ(一様乱数)
t = np.arange(1, 20, Ts)
size = len(t)
u = np.random.rand(len(t))
u = u - np.ones(size)*0.5
plt.figure()
plt.plot(t, u)
plt.xlabel("Time[s]")
plt.ylabel("Input[-]")
plt.grid()
plt.savefig("Input.png")
freq = np.fft.fftfreq(size, Ts)
sp = np.fft.fft(u)/size
freq = freq[1:size//2]
sp = sp[1:size//2]
pw = 20*np.log10(np.abs(sp))
plt.figure()
plt.plot(freq, pw)
plt.xscale('log')
plt.xlabel("Freqency[Hz]")
plt.ylabel("Magnitude[dB]")
plt.grid()
plt.savefig("Input_FFT.png")
# 出力：フィルタに対する応答
y, t_o, x_0 = lsim(P[0], u, t)
plt.figure()
plt.plot(t, y)
plt.xlabel("Time[s]")
plt.ylabel("Output[-]")
plt.grid()
plt.savefig("Output.png")
freq = np.fft.fftfreq(size, Ts)
sp = np.fft.fft(y)/size
freq = freq[1:size//2]
sp = sp[1:size//2]
pw = 20*np.log10(np.abs(sp))
plt.figure()
plt.plot(freq, pw)
plt.xscale('log')
plt.xlabel("Freqency[Hz]")
plt.ylabel("Magnitude[dB]")
plt.grid()
plt.savefig("Output_FFT.png")
LEDを1/fで揺らがせる
電気回路
マイコンはESP32を用います．ESP32のPWM生成器を26pinと繋ぎ，LEDの光量をPWMによって制御します．また，GPIOの出力では心持たないこともあり，トランジスタを介して，LEDをON/OFFします．
マイコン実装
マイコンに実装するため，Python上で上記のフィルタを離散化し，離散化フィルタ(IIRフィルタ)のパラメータを計算します，
P_d = c2d(P[0] , Ts, method='tustin')
num, den = tfdata(P_d)
num = num[0][0]
den = den[0][0]
print(P_d)
print(num)
print(den)
0.0214 z^6 - 0.1036 z^5 + 0.2039 z^4 - 0.2073 z^3 + 0.1132 z^2 - 0.03073 z + 0.003058
-------------------------------------------------------------------------------------
z^6 - 5.315 z^5 + 11.68 z^4 - 13.56 z^3 + 8.768 z^2 - 2.983 z + 0.4161
dt = 0.005
[ 0.03340178 -0.11369083 0.13390645 -0.05034986 -0.01533875 0.0141659 -0.00209462]
[ 1. -4.19253354 6.83645287 -5.29849454 1.80667153 -0.10096111 -0.05113514]
ESP32では，rand()で一様乱数を生成し，先ほど求めた離散化フィルタに通すことで，ピンクノイズを生成します．それに比例してduty比を変化させることで，LEDの光量を1/fで揺らがせます．
#include <Arduino.h>
#include "Time.h"
#define NUM 7
// -----------------------------------------------------------
// Paramter
// - PWM
const uint8_t bit_pwm = 8;
const double duty_max = pow(2, bit_pwm);
const double freq_pwm = 1200;
const uint8_t pwm_ch = 0;
// - LED Pin
const uint8_t led_pin = A19;
// - Time
const unsigned long time_int = 5000;
unsigned long time_prv = 0;
// - Filter
const double a[] = { 0.03340178, -0.11369083, 0.13390645, -0.05034986, -0.01533875, 0.0141659 , -0.00209462};
const double b[] = { 1.0 , -4.19253354, 6.83645287, -5.29849454, 1.80667153, -0.10096111, -0.05113514};
double u[NUM];
double y[NUM];
// - Random
const long rand_max = 100000;
// -----------------------------------------------------------
// Setup
void setup() {
// PWMの初期設定
ledcSetup(0, freq_pwm, bit_pwm);
ledcAttachPin(led_pin, pwm_ch);
// フィルタの初期化
for (size_t i = 0; i < NUM; i++)
{
u[i] = 0.0;
y[i] = 0.0;
}
}
// -----------------------------------------------------------
// Loop
void loop() {
// 時間周期の調整
while((micros() - time_prv) < time_int){
}
time_prv = micros();
// フィルタの更新
// - 内部変数配列を右シフト
for(size_t i = (NUM - 1); i > 0; i--){
u[i] = u[i - 1];
y[i] = y[i - 1];
}
// - 入力uに一様乱数をセットする
u[0] = (double)(random(rand_max) - rand_max/2)/(double)rand_max;
// - 出力yの計算
y[0] = 0;
for(size_t i = 0; i < NUM; i++){
y[0] += b[i]*u[i] - a[i]*y[i];
}
y[0] = y[0]/a[0];
// 生成されたピンクノイズに合わせてduty比を変化させる
double gain = 10.0;
double duty = (y[0]*gain + 0.5)*duty_max;
ledcWrite(pwm_ch, duty);
}