ホワイトノイズを生成する

最近まで、ホワイトノイズを作るには、一様乱数を使えばよいと思っていた。
たしかに一様乱数でも、ちゃんと周波数特性が均一なホワイトノイズになる。シンセやAudacity等でホワイトノイズを生成すると、たいてい一様乱数が出てくる。

しかし、信号処理の世界では正規乱数（値が正規分布的に分布する乱数）を使うほうが一般的らしい。
正規乱数によるホワイトノイズのことを「ホワイトガウシアンノイズ」と呼ぶこともある。

なぜ一般的な信号処理で「正規乱数によるノイズ信号」（ガウシアンノイズ）を使いがちなのかというと、正規乱数は「最も平均的なノイズ信号」といえるからだ。
あらゆる雑音信号をたくさん集めて平均していくと、集めれば集めるほど正規乱数に近づいていく（これは中心極限定理から言える）。
また、ガウシアンノイズは数学的に扱いやすい性質がいろいろとあり（ここでは詳細は省く）、数式上の利便性もある。

調べたところ、正規乱数のメジャーな生成方法には、おおよそ3種類があるらしい。


1. 一様乱数を12回足して6を引く（近似）

中心極限定理をそのまま利用した方法。
一様乱数をたくさん足し合わせて平均すれば、正規乱数に近いものが得られる。

近似なので精度は高くないが、ちょっとした用途ならこれで十分かもしれない。
また、乱数を1つ作るのにループを要するため、大量に生成する場合には遅いという話も聞く。

亜種として「16回足して8を引く」実装もある（例: Web上で音声合成ができる PinkTrombone）。
要するに $2N$ 回足して $N$ を引けばよい。


2. ボックス・ミュラー法

2つの一様乱数 $ 0 < a, b < 1 $ をもとに、2つの正規乱数 $ x, y $ を生成する方法。
$ x = \sqrt{-2 \log a} \cos 2 \pi b $
$ y = \sqrt{-2 \log a} \sin 2 \pi b $

おそらく最も教科書的でよく知られている方法で、実装も簡単。
三角関数・対数・平方根といった数学関数を使うことになるので、オーディオ信号生成のように大量の乱数が必要な用途ではもっと高速化したくなるケースがあるかもしれない。


3. ジッグラト・アルゴリズム

超ざっくりいうと、ルックアップテーブルを使う方法。
この名前は、古代メソポタミアの「ジッグラト」(Ziggurat) のように、確率密度関数を階段状に区切ったテーブルを使うことに由来するらしい。

ボックス・ミュラー法と比べて、高速に動作するのがウリらしい。
Julia の randn 関数は、このアルゴリズムで実装されている（ソースコード）。

上述の2手法と比べると少し複雑なので、説明がてら試しに実装してみる。

基礎: 棄却サンプリング

ジッグラド・アルゴリズムの前に、まずはその基礎をなす「棄却サンプリング (Rejection sampling)」という手法について説明する。

これは、重み付き乱数（確率的な重みがついて偏った値が出てくる乱数）を生成する方法のひとつで、「棄却」に「サンプリング」などという仰々しい名前をしているが、考え方と実装方法はとてもシンプルだ。

重み付き乱数 $x$ の分布を $p(x)$ として、この曲線をグラフ上で考えてみる。
この関数 $p(x)$ は、出てきやすい $x$ のところで値が大きくなるような関数だとする。
たとえば、$p(x)$ が下のような曲線を描くとすれば、乱数 $x$ は 0 付近の値が出がちな乱数になる。
正規乱数の場合は、関数 $p(x)$ をガウス関数（たとえば$\exp(-x^2)$など）にすればよい。

さて、このグラフ上にランダムな点 $(x, y)$ を打ってみる。
このとき、点の位置と $p(x)$ の曲線を比べて「点が下にある場合」と「点が上にある場合」の2パターンが出てくる。

打った点が $p(x)$ より「下」にある場合は、その点の $x$ 座標を重み付き乱数として、アルゴリズムを終える（受理）。
逆に、点が $p(x)$ より「上」にあった場合、その点はなかったことにして（棄却）、もういちどランダムな点を打ちなおす。

このような処理をおこなえば、うまいこと確率的な重みがついた乱数 $x$ が生成できる。

ジッグラト・アルゴリズム

棄却サンプリングをもっと効率よくおこなうのが、ジッグラト・アルゴリズムである。

実は棄却サンプリングには無駄がある。
ランダムに打った点が「棄却」と判定された場合、点を何度も打ち直すはめになる。条件が悪ければ、何度でもループする可能性がある。そのような場合、関数 $p(x)$ の値を何度も計算しなおす必要があり、あまり効率が良くない。

（すみません！ 書きかけです！）

おまけ: 一様乱数と正規乱数の性質

値の範囲

言うまでもなく、一様乱数は一定の範囲内（たとえば -1 ～ +1）に収まる。

正規乱数は、理論上は無限小～無限大の範囲で散らばる。
とはいえ、標準正規分布（平均 $0$、分散 $1$ の正規分布）の場合、-3 ～ +3 の範囲内に約 99.7% が収まる計算になる。
つまり標準正規乱数を 3 ～ 4 程度の値で割れば（およそ -9.5 〜 -12 dB 程度）、ほとんどの信号値が -1 〜 +1 の範囲内に収まってくれる。

RMS

オーディオ用途では、信号の平均的な音量を知りたいことがある。

平均的な音量の指標として使われるのが Root Mean Square (RMS) である。
分布の平均が $0$ の場合、RMSと標準偏差の計算式が一致する。この性質を使って、ノイズ信号のRMSの理想値が計算できる。

一様乱数（-1 ～ +1）のRMSは $ 1 / \sqrt{3} \approx 0.57735 $（参考）、すなわち約 -4.7712 dB と計算できる。

標準正規乱数の標準偏差は $ 1 $ と定義されているので、RMSは $ 1 $ 、すなわち 0 dB である。


おまけ: ベルベットノイズ

人間にとって知覚的に「なめらか」に聞こえるノイズとして、ベルベットノイズ (velvet noise) という特殊な作り方のホワイトノイズがあるらしい。

もともとはリバーブ用に開発された信号で、元論文の著者のサイト (Velvet Noise Reverbration Demos) でガウシアンノイズとベルベットノイズを聴き比べられる。
たしかに「Gaussian white noise」と「Velvet noise (略) 7000」を比べると、後者のほうがなめらかに聴こえる。


コメント

コメントを投稿

このブログの人気の投稿

基本波形の生成

イメージ
オーディオやマルチメディア関連のプログラムを書いていると、手早くノコギリ波や矩形波を生成したいことがある。 そこで、簡単にいろいろな波形を生成する方法をまとめる。 1. ノコギリ波 1サンプルごとに少しずつ値を足していくと、直線的に増加していく波形ができる。 値が 1.0 を超えた瞬間に 0.0 へ戻すことで、ノコギリ状の周期波形が生成できる。 これは、値の小数点以下のみを取り出すことで実装できる。 実装してみる。 以下のコードでは、オーディオ信号として扱いやすいよう出力値の範囲を 0〜1 から -1〜1 へ変換している。 double sawtooth(double 周波数, double サンプルレート) { static double phase = 0.0; // 初期値 phase += 周波数 / サンプルレート; phase -= floor(phase); // 整数部分を引き算 return 2.0 * (phase - 0.5); // -1.0～1.0 の範囲に広げる } 「小数点以下を取り出す」を「1.0で割ったときの剰余」として実装する方法もある。 double another_sawtooth(double 周波数, double サンプルレート) { static double phase = 0.0; // 初期値 phase += 周波数 / サンプルレート; phase = fmod(phase, 1.0); // 小数部分を取り出す return 2.0 * (phase - 0.5); // -1.0～1.0 の範囲に広げる } プログラミング言語によっては（C# など）、浮動小数点数（float や double）に対して剰余演算子 % が使えるので、この実装方法のほうがシンプルに書ける。 2. 矩形波 ノコギリ波が 0.5 を超えたら 1.0、超えなければ -1.0 を出力する。 条件を「0.75 を超えたら 1.0 を出力」に変えれば、デューティ比25%の矩形波になる。 double pulse(double 周波数, double サンプルレート) { static double phase =
続きを読む

(1)C言語で音声合成もどき ～WAVファイルを生成する～

この記事は素人が音声合成で遊んでいるだけの記事です。完全に行き当たりばったりなので、紹介している内容の保証はできません。 また、記事の内容を予告なく変更することがあります。 「C言語で音声合成もどき」ということで、このシリーズではフォルマント合成と呼ばれる手法を使った音声合成プログラムをC言語で作っていきます。 本題の音声合成に入る前に、まずは肩慣らしとして 音声ファイルを作るプログラム をつくります。 今回は  8bit, 44100Hz, 1ch  の WAV形式 のファイルを作成するプログラムを書いてみます。 ( この記事では、比較的実装が簡単なWAV形式を扱います。 ひょっとするとMP3形式等のほうが馴染み深いかもしれませんが、MP3ではデータ圧縮にやや複雑なアルゴリズムを使うことになります。) なお、WAVファイルの構造については WAVファイルフォーマット WAVE file format を参考にしました。
続きを読む

(2)C言語で音声合成もどき ～母音の生成～

イメージ
この記事は素人が音声合成で遊んでいるだけの記事です。完全に行き当たりばったりなので、紹介している内容の保証はできません。また、記事の内容を予告なく変更することがあります。 生成した「お」の波形 今回は本題の音声合成に移っていきます。 前回 は波形データをWAVファイルに書き出す部分を作ったので、ここからはさっそく母音の合成をしていきます。 ざっくりいうと人間の声は、 声帯 で基となる音がつくられ、 声道 で加工されたものです。 今回のプログラムは、この「声帯→声道」という構成を参考にして音声を合成していきます。 まずは、声帯にあたる部分から実装することにします。 声帯が振動すると、周期的な波が発生します。 つまり、なんらかの周期的な波を使えば模倣できそうです。 一般的なシンセサイザーのようにノコギリ波や三角波を使うこともできますが、ここでは声帯振動の波を模倣した Rosenberg波 とよばれる波形を使います。 以下がRosenberg波の式です。 出典： https://dspace.jaist.ac.jp/.../896paper.pdf のP14 (2.1)、一部式の変形・表記変更あり τ1は 声門開大期 *1 の長さ 、τ2は 声門閉小期 *2 の長さ を表します。 *1 声門開大期……声門が閉じた状態から開ききるまでの時間 *2 声門閉小期…… 声門が開いた状態から閉じきるまでの時間 これらの値は、τ1 + τ2 が 0〜1 の範囲におさまるように決めます。 今回は  τ1 = 0.90, τ2 = 0.05 としました。 Cで書くとこんな感じです。 double GenRosenberg(double freq) { /* Rosenberg波を生成 */ static double t = 0; double tau = 0.90; /* 声門開大期 */ double tau2 = 0.05; /* 声門閉小期 */ double sample = 0.0; t += freq / (double)サンプリング周波数; t -= floor(t); if (t <= tau) {
続きを読む