Walsh変換、フーリエ変換を用いた正規分布擬似乱数の生成手法
Walsh transform (WT) and Fourier transform (FT) generating normal/Gaussian
distribution sequences
2025.4.13 関連文献として下記の資料を見つけました． 直交変換にこのような応用があるとは知りませんでした．
直交変換による悪条件行列の改善法 小林 富士男, 福山大学工学部紀要 第10号, 89-95, 1988年3月
DFTによる悪条件行列の改善 小林 富士男, 尾関 孝史, 筒本 和広, 福山大学工学部紀要 第30巻, 191-199, 2006年12月
2025.4.11 関連文献として下記の資料を見つけました．
逆フーリエ変換を利用した任意スペクトルをもつ正規信号の発生 青島 伸治, 宮川 康彦, 計測自動制御学会論文集 15 (3), 389-394, 1979 PDFダウンロード
著者名を見て、もしやと思って調べたら TSP(Time Streched Pulse)/ATSP(Aoshima's Time Stretched Pulse) の青島氏でした．
2025.4.8 先行文献として下記の資料を見つけました． このような文献が今まで『埋もれていた』とは驚きです．
アダマール変換を利用する多数の無相関な正規分布の不規則信号発生法 泉 照之, 計測自動制御学会論文集 27 (3), 365-367, 1991 PDFダウンロード
上記文献に『一様乱数をアダマール変換して正規乱数を発生させる方法がある5)．』とありますが、参考文献の「5) 吉澤 正：Yatesの算法を応用した擬似正規乱数について、日本オペレーションズ・リサーチ学会研究会アブストラクト集、83/84
(1969)」に関する詳細情報はネット上には無いようです． 有償の科学技術文献複写サービスを利用するか、大学図書館・国会図書館等で調べるしかありません．（2025.4.9追記）
参考文献 5) の著者名で調べると下記の情報がありました． 後者のコラムの末尾に『統計にはYatesの算法といわれるものが知られているが，詩人Yeatsと統計家Frank Yates（1902-1994）ではつづりが違うことがおもしろい』とあります．（2025.4.13追記）
(株)日科技研：吉澤正先生がご逝去されました
第41話 アイルランドの山の姿から —台地分布を考える—（六一学者の千字一話）
Yatesの算法とは Yates analysis のことでしょうか？ 他にYatesの業績とHadamard変換の関連についた述べた英文の情報がいくつかありますが、統計学の方面はまったく専門外なので良くわかりません．（2025.4.16追記）
上記文献の中に「選点直交多項式」という用語が出てきますが懐かしい限りです． 東大の技官になった同級生のH君が無線関連の卒研で使っていた手法で、指導教官の先生は「選点直交関数を見つけてきて、その応用を考えたのはH君の功績だ」と話していたものです．
計測制御分野の古い関連文献を立て続けに２つ見つけましたが、時代背景として 計測用途の連続時間の試験信号をマイコンとD/Aコンバータを用いて生成する～という試みが活発におこなわれていたことが考えられます． 探せばまだまだ有用な文献をサルベージ出来るかもしれません．
- Walsh変換、フーリエ変換を用いて正規分布擬似乱数を生成出来ることが分かったので、その手法をご紹介します．（偶然の発見です） 一様分布乱数を正規分布乱数へ簡単に変換出来ます．
- 厳密な解析・特性評価等はまだおこなっていません．
- 特許出願等はおこなっていませんので、どなたでもご自由にご利用いただけます． 学会発表の予定もありません． 学生さんの卒論・修論のネタなどにもご利用ください．
| Walsh変換を用いた正規分布擬似乱数の生成（WT法） |
- Nヶの平均 0 の一様分布乱数列をWalsh変換すると、Nヶの正規分布擬似乱数列が得られます． 高速Walsh変換を用いるためにデータ長Nは2の冪（N=2M）とします．
| フーリエ変換を用いた正規分布擬似乱数の生成（FT法） |
- Nヶの平均 0 の一様分布乱数列２セットから複素数列を生成してフーリエ変換します． 変換結果の実部と虚部をそれぞれ取り出すと、Nヶの正規分布擬似乱数列が２セット得られます． 高速フーリエ変換を用いるためにデータ長Nは2の冪（N=2M）とします．
- WT法・FT法で一様分布乱数列を正規分布に変換する基本原理は中心極限定理（CLT, Central Limit Theorem）を用いた手法と同様です． 違いは単純な加算ではなく、重み付け係数付き加算をしていることです．
- Nヶの一様分布乱数の線形結合によりNヶの正規分布乱数 g[n] を作り出しているので、出力 g[n] 間には相関が生じているのではないか？（乱数としての特性が劣化しているのではないか？）という疑義が生じますが、詳細は未確認です． データ長Nが十分に大きければ実用上の問題が生じるほどではない～と期待したいところですが.
. . . .
- 擬似乱数の検定方法として離散フーリエ変換検定というのがあるそうですが、それとWT法・FT法の技術的関連性はよく分かりません．（2025/3/28追記）
- 一様分布乱数の代わりに、-1/1の２値の擬似乱数またはオフセットを除去した２値の擬似乱数列（平均0）を用いても正規分布乱数（シミュレーション結果）を生成できます．（2025/3/29追記）
- 見方を変えると直交変換法（WT法、FT法）は以下のような問題であると考えられます．（2025.3.31追記、2025.4.1修正、2025.4.3修正）
問題の記述は非常に単純ですが、厳密な数学的証明も簡単に出来るでしょうか？
- NヶのN次元直交ベクトル cn があるとき
cn={ cn0, cn1, cn2, ... , cn(N-1) } （n=0, 1, 2, ... , N-1）
cn*ck=0 （n≠k）
- 平均0の乱数unを要素に持つN次元ベクトル u と cn との内積を考える
u={ u0, u1, u2, ... , uN-1 }
u0+u1+u2+ ... +uN-1=0
- Nケの内積 gn はどのような特性を持つだろうか？（特に相関特性は？）
gn=u *cn （n=0, 1, 2, ... , N-1）
- 自明な例として直交ベクトル cn が
cnm=1 （n=m）
cnm=0 （n≠m）
であったならば gn=un となる
- では cn がWalsh変換やフーリエ変換などの直交変換の係数であった時の gn はどうなるか？
シミュレーションから gn は正規分布乱数となることが確認されていて、これは中心極限定理からも納得出来る結果ではあるが、数学的に厳密な証明や直感的に分かりやすい証明は出来るだろうか？
- 生成した正規分布乱数 gn に対してBox-Muller法と逆の変換をかけて、乱数検定をしたらどうなるだろうか？（相関特性の評価も出来るだろうか？）
- 処理を単純に周波数分析（Walsh解析・フーリエ解析）と解釈すると、入力信号がランダム・ノイズなら各周波数成分の振幅・位相もランダムで分析出力サンプル（変換後の正規分布擬似乱数サンプル）間の相関は無い～ということになりますが、擬似乱数の性質の議論としてはこれで良いのでしょうか？（2025.4.14追記）
| C言語のサンプル・プログラムとシミュレーション結果 |
| Scilabのサンプル・プログラムとシミュレーション結果 |
- 動作確認は Scilab Ver.6 でおこなっています．
- グラフはいずれも一様分布乱数のseedを変えて４回実行したシミュレーション結果を重ね書きしたものです．
- データ長はすべてN=220=1048576です．
上のグラフの形がかなり歪（左右非対称）に見えるので、一様分布乱数のseedを変えてもう一度プロットしてみました（2025.3.28追加）
（SEED='f0a50fa5', 'a5a53c3c', '0000ffff', 'ffff0000'）
一様乱数の代わりに-1/1の２値の擬似乱数を用いても正規分布乱数を生成できます．（2025/3/29追記）
左右のグラフは使用した２値の擬似乱数のseedが異なります．（Scilabプログラム pn_walsh_transform_method.sce）
- FT法のminimumなScilabシミュレーション・プログラム minimum_ft.sce （2025.4.13追加）
正規分布擬似乱数生成処理をしているのは実質的にたったの３行です． FFTが使えるなら、どのようなプログラミング言語にも移植できるでしょう．
/* Gaussian random number sequence generator by FFT */
/* (c) 2025 cepstrum.co.jp */
clear; // all clear
DATLEN=2^20; // length=2^20=1048576
rand('uniform'); // set uniform random number mode
rand('seed', hex2dec('a5a50f0f')); // set random number seed
//==========================================================================================
x=(rand(1:DATLEN)-0.5)+%i*(rand(1:DATLEN)-0.5); // generate complex uniform random number
y=fft(x, -1); // FFT
gauss_rand=[real(y), imag(y)]; // complex ---> real conversion
//==========================================================================================
scf(0); clf;
plot2d(gauss_rand(1:3000), axesflag=1);
scf(1); clf;
histplot(-1500:50:1500, gauss_rand, axesflag=1);
savewave('out.wav', 0.9*gauss_rand/max(abs(gauss_rand)));
2025.4.7追加
- とりあえず生成した正規分布擬似乱数の自己相関をプロットしてみました．
- 生成方式により非常に微妙な違いがあるものの、WT法・FT法（直交変換法）に異常な特性が見られるということは無いようです．
WT法
FT法
Box-Muller法
中心極限定理に基づいた手法
- 今まで誰も気が付かなかったのが不思議なほど単純なWT法・FT法・直交変換法を発見するにいたった経緯は次のとおりです．
- 離散時間信号処理（デジタル信号処理）の分野のある特別な試験信号を生成する手法として、周波数領域で信号の振幅周波数特性・位相特性を定義して、その逆フーリエ変換により時間信号を得る～というものがあります．
- ならば同様に周波数領域でランダムな位相特性を定義して、それを逆フーリエ変換すれば特性の良いランダム・ノイズ／擬似乱数が得られるのではないか？と考えて試してみました． 試行錯誤を繰り返しているうちに、適当な擬似乱数をフーリエ変換すれば正規分布特性が得られることに気づき、さらにフーリエ変換の代わりにWalsh変換を使っても良いではないかということになり、最終的に中心極限定理に基づく正規分布乱数生成手法と同じことをやっているのだと理解出来ました．
- このような形で一様分布乱数、正規分布とWalsh変換・フーリエ変換が一つに結びつくとは興味深いことです．（ガウス関数のフーリエ変換がガウス関数になるという既知の性質はあるそうですが）
FFTの使い方（電子科学シリーズ91）、安居院猛、中嶋正之、産報出版、1981（ISBNコードが日本で使われるようになる以前の本です）
FORTRANの算術型GOTO文／計算型GOTO文を目にしたのは何十年ぶりでしょうか. . . .
- 高速フーリエ変換はFFTWを使うのが定番です． ScilabもFFTWを利用しています．
