ノイズ生成

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$ ノイズ）

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

主要な生成手法

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$) のいずれも統一的に扱える。

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

関連項目

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