ノイズ生成

1. 色付きノイズの分類

ノイズはパワースペクトル密度 $S(f)$ の周波数依存性によって分類される。 一般に $S(f) \propto 1/f^\alpha$ と表され、指数 $\alpha$ によって「色」が決まる。

2. ホワイトノイズ

ホワイトノイズは、パワースペクトル密度が全周波数にわたって一定である理想的な雑音である。 離散時間では各サンプル $x[n]$ が独立同分布 (i.i.d.) に従う系列として定義される。

統計的性質

• 自己相関関数: $R_{xx}[m] = \sigma^2 \delta[m]$
• パワースペクトル密度: $S_{xx}(e^{j\omega}) = \sigma^2$（一定）
• 確率分布: ガウス分布が最も一般的だが、一様分布なども使用される

生成方法

• 一様分布ホワイトノイズ — 疑似乱数生成器（メルセンヌ・ツイスタ、xoshiro など）を使用
• ガウスホワイトノイズ — Box-Muller 法や Ziggurat 法で一様乱数を正規分布に変換

Box-Muller 法

独立な一様乱数 $U_1, U_2 \sim \text{Uniform}(0,1)$ から、独立な標準正規乱数の組 $(Z_1, Z_2)$ を生成する：

3. ピンクノイズ（$1/f$ ノイズ）

ピンクノイズはパワースペクトル密度が周波数に反比例する雑音であり、 オクターブあたりのエネルギーが等しいという特徴を持つ。 音響測定、音楽制作、自然現象のモデリングに広く使われる。 このようなスペクトルを持つノイズは「$1/f$ 揺らぎ」とも呼ばれ、自然界に広く見られる。

オクターブ等パワー性の証明

ピンクノイズのパワースペクトル密度を $E[|X(f)|^2] = f^{-1}$ とする。 任意の中心周波数 $f_0$ から 1 オクターブ上の $2f_0$ までのパワーは：

結果は $f_0$ に依存せず一定である。 すなわち、どのオクターブ帯域でもパワーは $\log 2$ で等しい。 これがピンクノイズが音響測定に用いられる理由であり、 オクターブバンド分析器で測定すると全帯域のレベルが揃う。

主要な生成手法

3.1 理想的なピンクノイズフィルタ

連続時間においてホワイトノイズをピンクノイズに変換するフィルタの伝達関数は：

これは分数階微積分（fractional calculus）の演算子であり、 有理関数では正確に表現できない。 したがって、有限次の IIR フィルタで近似する必要がある。

近似の基本方針は、$1/\sqrt{s}$ を1次項の部分分数に展開することである：

極 $p_k$ を対数等間隔に配置し、留数 $\beta_k$ を調整すると $-3$ dB/oct（$= -10$ dB/dec）のスロープを広帯域にわたって近似できる。 Oustaloup の再帰近似法はこの方針を体系化したものであり、 極と零点を厳密に対数等間隔に配置する。

3.2 Paul Kellet の IIR フィルタ

Paul Kellet が musicdsp.org で公開したフィルタは、6段の1次 IIR フィルタを並列に合成して $-3$ dB/oct 特性を近似する。 係数は理論的な最適化ではなく、Visual Basic のプログラムで インパルス応答をプロットしながら手動で調整されたものである。

3段階の精度のバリエーションがあり、ここでは最高精度の Refined（Instrumentation grade）版を分析する。

差分方程式

ホワイトノイズ $x[n]$ を入力とし、各段のフィルタ状態 $b_k[n]$ を以下で更新する：

加えて1サンプル遅延の項 $b_6$ があり、出力は：

ここで $b_6[n] = 0.115926 \cdot x[n-1]$（$b_6$ は出力に加算された後に更新される）。

フィルタ係数

極周波数の概算値は $f_p \approx f_s(1 - |\alpha_k|) / (2\pi)$ で得られる。 隣接する極の周波数比は $\times 5.9$, $\times 4.6$, $\times 4.3$, $\times 3.4$ であり、 概ね対数等間隔だが一定ではない。これは手動チューニングの痕跡である。

C++ 実装例

// Paul Kellet pink noise filter (Refined / Instrumentation grade)
// 入力: white — ホワイトノイズサンプル ([-1, +1])
// 状態: b0〜b6 は呼び出し間で保持する (初期値 0)
double white = rng.nextDouble() * 2.0 - 1.0;  // 一様乱数
b0 = 0.99886 * b0 + white * 0.0555179;
b1 = 0.99332 * b1 + white * 0.0750759;
b2 = 0.96900 * b2 + white * 0.1538520;
b3 = 0.86650 * b3 + white * 0.3104856;
b4 = 0.55000 * b4 + white * 0.5329522;
b5 = -0.7616 * b5 - white * 0.0168980;
double pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362;
b6 = white * 0.115926;
// pink を出力として使用（必要に応じて正規化）

3.3 伝達関数 $H(z)$ の導出

各段の伝達関数

差分方程式 $b_k[n] = \alpha_k \cdot b_k[n-1] + \beta_k \cdot x[n]$ の Z 変換を取ると：

したがって各段の伝達関数は：

$b_6$ の伝達関数

$b_6$ は出力の加算に使用された後に $x[n]$ で更新される。 すなわち、加算時点での $b_6$ の値は前サンプルの $x[n-1]$ に基づく：

全体の伝達関数

出力 $y[n] = \sum_{k=0}^{5} b_k[n] + b_6[n] + 0.5362 \cdot x[n]$ より：

各項を展開すると：

3.4 周波数特性 $H(\omega)$ の導出

$z = e^{j\omega}$ を代入する。各 IIR 段の分母は：

各段の実部・虚部

$H_k(e^{j\omega}) = \frac{\beta_k}{1 - \alpha_k e^{-j\omega}}$ の実部と虚部を分離すると：

$b_6$ + 直達項

全体の応答

振幅特性・位相特性

ここで $\omega = 2\pi f / f_s$（$f_s$: サンプリング周波数）。

各段の振幅応答

個々の IIR 段の振幅応答は以下の閉じた形で表される：

ただし全体の振幅応答は各段の振幅の単純な和ではない。 位相の異なる複素数の和であるため、実部と虚部を別々に合算してから 絶対値を取る必要がある。

周波数特性のプロット

$f_s = 44100\;\text{Hz}$ における Kellet フィルタの振幅特性と、理想的な $-3$ dB/oct からの偏差を以下に示す。

• 10 Hz – 4 kHz: 偏差 $\pm 0.1$ dB 以内 — ほぼ完璧な近似
• 4 kHz – 10 kHz: 最大 $-1$ dB 程度 — わずかに減衰が大きい
• 10 kHz 以上: Nyquist に向かって急落（IIR の限界）

6つの1次フィルタのみで可聴帯域のほぼ全域をカバーしている、非常に優れた近似である。

3.5 サンプリング周波数依存性

Kellet フィルタの係数は固定であるため、 サンプリング周波数 $f_s$ が変わると極周波数がすべてシフトし、 近似の精度が変化する。 極周波数は $f_p \approx f_s(1 - |\alpha_k|) / (2\pi)$ に比例するため：

Lubomir によるサンプルレート別係数（2009年）

musicdsp.org の同ページに Lubomir が投稿した拡張版では、 各段の目標極周波数 $f_k$ をサンプルレート帯域ごとに定義し、 $\alpha_k = e^{-2\pi f_k / f_s}$ で極係数を算出する方式が示されている。 ただしゲイン係数 $\beta_k$ は提供されておらず、 出力は各段の単純な加算 b0+b1+b2+b3+b4+b5+white-b6 となっている。 Kellet のような精密なゲイン調整が行われていないため、 スペクトルの平坦性は Kellet オリジナルより劣る可能性がある。

高サンプルレートほど極周波数が密集しており、Lubomir 自身も 「高サンプルレートではフィルタの品質が劣化する」と指摘している。 これは 6 段の IIR フィルタという構造的な限界に起因するものであり、 極の数を増やさない限り本質的な改善は難しい。

対策の選択肢

1. 係数テーブル切替: 上記の Lubomir 方式で $f_s$ 帯域ごとに係数を切り替える
2. 動的係数計算: 目標極周波数を固定し、$\alpha_k = e^{-2\pi f_k / f_s}$ で任意の $f_s$ に対応する
3. IFFT 法: 周波数領域で直接スペクトルを構成し、$f_s$ 依存性を完全に回避する（次節で解説）

3.6 IFFT 法 — サンプリング周波数非依存の生成

IIR フィルタのサンプリング周波数依存性を完全に回避する方法として、 周波数領域でスペクトルを直接構成し、逆 FFT する手法がある。

原理

パワースペクトル密度 $S(f) \propto 1/f$ を持つ信号の振幅スペクトルは $|X(f)| \propto 1/\sqrt{f}$ である。 ここに一様ランダムな位相 $\theta_k \in [0, 2\pi)$ を与え、逆 FFT すると 理想的なピンクノイズが得られる：

周期性と必要な長さ

IFFT で生成された信号は周期 $N$ の周期信号であるため、 ループ再生しても接続点で不連続が生じない。 ただし人間の聴覚は雑音の繰り返しに対して意外と敏感であり、 周期が短いと周期性を知覚する可能性がある。

$N = 2^{20}$（約24秒）で 4 MB — 現代のシステムでは問題にならないサイズである。

IIR フィルタ法との比較

C++ 実装例（IFFT 法）

// IFFT 法によるピンクノイズバッファ生成
const int N = 1 << 20;  // 2^20 ≈ 24秒 @44.1kHz
std::vector<std::complex<double>> spectrum(N);
std::mt19937 rng(42);  // シード固定で再現性確保
std::uniform_real_distribution<double> dist(0.0, 2.0 * M_PI);

// DC = 0
spectrum[0] = 0.0;

// 正の周波数: 振幅 1/√k, ランダム位相
for (int k = 1; k < N / 2; ++k) {
    double amp = 1.0 / std::sqrt((double)k);
    double phase = dist(rng);
    spectrum[k] = std::polar(amp, phase);
}

// Nyquist (実数)
spectrum[N / 2] = 1.0 / std::sqrt((double)(N / 2));

// 共役対称 → 実数出力を保証
for (int k = 1; k < N / 2; ++k)
    spectrum[N - k] = std::conj(spectrum[k]);

// 逆 FFT → 時間領域信号
std::vector<double> buffer(N);
ifft(spectrum.data(), buffer.data(), N);

// ピーク正規化
double peak = *std::max_element(buffer.begin(), buffer.end(),
    [](double a, double b){ return std::abs(a) < std::abs(b); });
double norm = 1.0 / std::abs(peak);
for (auto& s : buffer) s *= norm;

// 使用時: buffer[n % N] でループ再生

任意の色ノイズへの拡張

IFFT 法の大きな利点は、振幅スペクトルを $|X[k]| \propto 1/k^{\alpha/2}$ とするだけで 任意の $\alpha$ を持つ色付きノイズを同一の手法で生成できることである。 ピンク ($\alpha = 1$)、ブラウン ($\alpha = 2$)、ブルー ($\alpha = -1$)、 バイオレット ($\alpha = -2$) のいずれも統一的に扱える。

3.7 シェルビング型フィルタの縦続接続法

Kellet フィルタが並列合成であるのに対し、 シェルビング型フィルタを縦続接続して $-3$ dB/oct 特性を得る方法がある。 初期値をオクターブ間隔のシェルビング特性から構成し、 非線形最適化によって高精度な近似を実現する。

基本的な考え方

高域で $-3$ dB 減衰するシェルビング型フィルタを 11 段用意し、 各段のカットオフ周波数を Nyquist 周波数の $1/2, 1/4, 1/8, \ldots$ とオクターブごとに配置する。 これらを縦続接続すると、オクターブ上がるごとに $-3$ dB のシェルビング特性が 1 段分ずつ追加されるため、全体として概ね $-3$ dB/oct のスロープが得られる。

ただし、シェルビングフィルタの遷移帯域幅が有限であるため、 このままでは高域の近似誤差が大きい。 そこで初期値から出発して非線形最適化を行い、 20–20,000 Hz の範囲で $f^{-1/2}$ 特性に対する近似精度を高める。

最適化手法

非線形最適化には Hooke-Jeeves の直接探索法と Nelder-Mead のシンプレックス法を交互に適用する。 どちらか一方ではローカル・ミニマムに陥りやすいため、 二つの手法を交互に用いることで脱出しやすくしている。

フィルタ構成

図3 の 1 次 IIR フィルタを 11 段縦続接続して構成する。 各段の伝達関数は：

全体の伝達関数は各段の積となる：

フィルタ係数

サンプリング周波数 44.1 kHz 時の最適化された係数を以下に示す。

20 Hz でゲインが 1 になるよう、1 段目の $a_0, a_1$ で調整したため、 この段のみ係数の絶対値が 1 を超えている。 固定小数点で実現する場合は、ゲインを全段に分散すべきである。

周波数特性

最適化後の振幅特性と理想 $-3$ dB/oct からの偏差を図4 に示す。

20–20,000 Hz の範囲で $\pm 0.02$ dB 以内の近似精度が得られている。 Kellet フィルタ（6 段並列、$\pm 0.1$ dB）と比べて 5 倍以上高精度である。

Kellet フィルタとの比較

3.8 2 次 IIR 縦続接続法

3.7 節のシェルビング型フィルタ（1 次 IIR）を 2 段ずつペアリングして 2 次 IIR（biquad）セクションにまとめ、 $L_p$ ノルム・ホモトピー法と交互最適化を 組み合わせることで、ほぼ等リップルな高精度近似を実現した。

設計手法

1. 初期値の構成: 3.7 節で得た 11 段の 1 次 IIR を 2 段ずつ畳み込んで 5 個の 2 次セクションにまとめ、残り 1 段を退化した 2 次セクション （$b_2 = a_2 = 0$）として扱う。計 6 セクション、フィルタ次数 12。
2. $L_p$ ノルム・ホモトピー法: 評価関数に $L_p$ ノルムを用い、$p$ を段階的に $2 \to 4 \to 8 \to \cdots \to 1024$ と増加させる。 $$\min_c \left(\frac{1}{N}\sum_{i=1}^{N} |e_i(c)|^p\right)^{1/p}$$
   • $p = 2$: 二乗誤差の最小化に相当し、目的関数が滑らかなため大域的な最適解の谷を見つけやすい
   • $p \to \infty$: 最大誤差の最小化（minimax）に近づき、ピーク誤差の等リップル化が進む
   • 前の $p$ での解を次の $p$ の初期値にすることで、滑らかに minimax 最適解へ遷移する
   • $p = 2 \to 1024$ のサイクルを複数回繰り返すと、 2 周目以降の $p = 2$ で再び滑らかな目的関数から探索し直すため、 さらに良い解の谷を発見できる
3. 交互最適化（各 $p$ の段階で実施）:
   • 分子係数 $(b_0, b_1, b_2)$ を固定して分母係数 $(a_1, a_2)$ を最適化（安定性制約付き）
   • 分母係数を固定して分子係数を最適化
   • これを繰り返す（Nelder-Mead シンプレックス法）
   全 30 変数を同時に最適化するとローカル・ミニマムに陥りやすいが、 分子（18 変数）と分母（12 変数）に分割することで収束が安定する。
4. 等リップル仕上げ: $L_p$ ホモトピーで minimax の谷に入った後、ペナルティ付きコスト関数 $$J(c) = \|e\|_\infty + \lambda\bigl(\max_k |e_k| - \min_k |e_k|\bigr)$$ で $\lambda$ を $20 \to 0.5$ に減衰させながら最適化する。 第 2 項がリップルの不均一度を罰するため、 minimax 誤差を維持しつつリップル振幅を等しくする効果がある。

フィルタ構成

各セクションの伝達関数は：

全体の伝達関数は各セクションの積：

フィルタ係数

サンプリング周波数 44.1 kHz 時の最適化された係数を以下に示す。

すべての極は実数軸上にあり、$|z| < 1$ を満たすため安定である。

周波数特性

$L_p$ ホモトピー法による最適化後の振幅特性と、 交互最適化のみの結果との比較を図5 に示す。

$L_p$ ホモトピー法で minimax の谷に入った後、 等リップルペナルティ付きコスト関数 $J(c) = \|e\|_\infty + \lambda(\max |e_k| - \min |e_k|)$ で仕上げることで、最大誤差 $\pm 0.0075$ dB と等リップル比 0.80 を両立した。 交互最適化のみ（$\pm 0.011$ dB）に対して 32% の誤差削減である。 誤差曲線のリップルはほぼ等振幅であり、 チェビシェフ近似の等リップル最適解にほぼ到達していることを示唆している。

最適化結果

$L_p$ ホモトピー法を 5 サイクル（$p = 2 \to 1024$ を 5 周）実行し、 等リップルペナルティ（$\lambda = 20 \to 0.5$）で仕上げた。 最適化の収束過程と最終結果を以下に示す。

12 ピークすべてで符号が $+, -, +, -, \ldots$ と交代しており、 チェビシェフの交替定理が要求する条件を満たしている。 ただし低域側（20–500 Hz）の 6 ピークが $\pm 0.006$ dB 前後であるのに対し、 高域側（1–20 kHz）では $\pm 0.0075$ dB とやや大きく、 等リップル比は 0.80 にとどまる。 7 段に増やすことでリップル振幅の均一性がさらに向上する（3.9 節参照）。

各手法の比較

3.9 2 次 IIR 縦続接続法（7 段、14 次）

3.8 節の 6 段 biquad（12 次）では符号交代は達成されたが、 等リップル比が 0.80 にとどまり、リップル振幅の均一性が不十分であった。 7 段目（14 次）を追加して自由度を増やすことで、 等リップル比 0.96 のより均一な等リップル近似を達成した。

設計

3.8 節の最適化済み 6 段係数を初期値とし、7 段目を通過特性（$H_7(z) = 1$）で初期化して 同じ $L_p$ ホモトピー＋等リップル仕上げを適用した。 最適化の結果、7 段目は Nyquist 付近の応答を補正するセクションとして機能している。

フィルタ係数

周波数特性

12 個のリップルすべてで符号が $+, -, +, -, \ldots$ と厳密に交代しており、 チェビシェフの交替定理が要求する条件を満たしている。 6 段でも符号交代は達成されているが、等リップル比は 0.80 にとどまる。 7 段では等リップル比が 0.96 に向上し、リップル振幅がより均一になっている。

6 段との比較

段数による精度の変化

段数を変えて同じ $L_p$ ホモトピー＋等リップル最適化を適用した結果を以下にまとめる。 5 段（10 次）では完全等リップル（ratio = 1.000）を達成したが、精度は $\pm 0.141$ dB にとどまる。 段数を増やすごとに精度が向上し、7 段（14 次）では $\pm 0.006$ dB と Kellet フィルタの約 17 倍の精度を実現する。

5 段で完全等リップル（ratio = 1.000）が得られた理由は、10 次の自由度に対して リップル数が 6（= $\lfloor N/2 \rfloor + 1$）と十分少なく、 すべてのリップル振幅と符号を制御するのに十分な極・零点があるためである。 6 段・7 段でも符号交代は達成されているが、リップル数が 12 に増加するため 振幅の均一性は低下する（等リップル比: 6 段 0.80、7 段 0.96）。 実用上は 7 段の $\pm 0.006$ dB が最も高精度な選択肢である。

4. ブラウンノイズ（$1/f^2$ ノイズ）

ブラウンノイズ（ブラウニアンノイズ、レッドノイズとも呼ばれる）は、 ホワイトノイズの累積和（ランダムウォーク）として生成される。 ロバート・ブラウンのブラウン運動に由来する。

• 生成: $x[n] = x[n-1] + w[n]$（$w[n]$: ホワイトノイズ）
• パワースペクトル密度: $S(f) \propto 1/f^2$（$-6$ dB/oct）
• 用途: 地震波モデリング、リラクゼーション音源

IFFT 法で生成する場合は、振幅スペクトルを $|X[k]| \propto 1/k$ に設定すればよい （パワースペクトル密度が $1/k^2$ に対応する）。

5. 疑似乱数生成器

ディジタルノイズの品質は疑似乱数生成器（PRNG）の性能に依存する。 高品質な PRNG はノイズ生成の基盤となる。

6. 応用

• 音響測定 — ピンクノイズによるスピーカー・ルームの周波数応答測定
• ディザリング — 量子化歪みの分散にブルーノイズやホワイトノイズを使用
• シミュレーション — モンテカルロ法における確率的モデリング
• 通信 — チャネルノイズのシミュレーション（AWGN チャネル）
• 音楽・サウンドデザイン — シンセサイザーの音源としてのノイズ波形
• インパルス応答測定 — TSP（時間引き伸ばしパルス）と組み合わせた室内音響計測

7. 参考資料

• Paul Kellet, "Pink noise filter" — musicdsp.org（Kellet フィルタの原典）
• Phil Burke, "DSP Generation of Pink Noise" — firstpr.com.au（各手法の包括的解説）
• "1st and 2nd order pink noise filters" — musicdsp.org
• Colors of noise — Wikipedia
• Pink noise — Wikipedia
• White noise — Wikipedia
• Brownian noise — Wikipedia
• Box–Muller transform — Wikipedia

関連項目

• 信号生成
• フィルタ設計
• ディジタルフィルタ設計
• フーリエ解析