「真の」と、まで言って良いものかどうか……。コンピュータでよく使われる、擬似乱数列ではなくて、自然界のノイズ(雑音)などの揺らぎを取り込んで、乱数を作成しようとする試みにトライしてみます。本ページでは、擬似乱数列に対して、ノイズを利用した乱数のことを真乱数と称することにします。
擬似乱数列と、その種については、「#6-7 Arduinoはサイコロを振る」をご参照ください。
今回の真乱数は、用途によっては、乱数を使用したい場面で毎回真乱数を発生させて、そのまま使用することもできるかもしれません。ただ、真乱数を発生させるのに一定の物理的時間が必要なケースが想定される(生成が遅い、時間がかかる)ので、基本的には、擬似乱数列の種を生成することを目標とします。
また、耐タンパ性については不問としますので、暗号などの用途には不向きであることも、書き添えておきます。
耐タンパ性について検索する(別ウインドウが開きます)
乱数とは、「規則性が無く、予測不能な数値」のことを指すそうです。今回の目的は、自然現象の中から、「規則性が無く、予測不能な数値」を取り出して、真乱数を生成することです。
では、Arduinoの周辺で、「規則性が無く、予測不能な数値」を持った自然現象を探してみましょう。
何か、違いますね？確かにどれも、測定値の不確実性はあります。周囲の状況(明るさや温度)によって、測定される値が変わるのですから。
では、短い時間間隔で、複数回、同じ事象を観測した場合は、どうでしょうか？周囲の状況が変化する早さに比べて、ずっと早い周期で同じ事象を観測すると、毎回観測される値は一定で、同じ値となります。これでは、次の値を予測できてしまいます。
今回必要なのは、観測する度に状態が変化して、違う値が読める、ノイズ(雑音)の様な自然現象なのです。
ノイズと言えば、「#6-4 オペアンプを使って増幅回路を作る」でサウンドセンサ(コンデンサマイク)のアナログ出力を反転増幅回路で増幅して、イヤホンで聞いた時は、すごいノイズでした。増幅率をもっと上げて、このノイズ成分に着目すると、乱数の原料に使えるでしょうか？
何だか良さそうです。次の問題は、真乱数を発生させるのに、わざわざ外付けにマイクと反転増幅回路(オペアンプ)を繋ぎますか？っていう、使い勝手の問題。こんな面倒な方法だったら、誰も使いたくないですものね？
もっと簡単で手軽な方法で真乱数を生成したい。
よく言われるのは、「フローティングになっているアナログ端子をanalogRead()する」方法です。これなら、未使用のアナログ端子を1本、フローティングにしておけば良いだけなので、手軽に使用できます。
では、フローティングノードのanalogRead()で、どの程度、値がばらつくのか、簡単なプログラムを組んで、測定してみましょう。
void setup (void)
{
// serial
Serial.begin (9600);
while (!Serial)
;
delay (10);
// measure
int idx;
for (idx=0; idx<1000; idx++)
Serial.println (analogRead (0));
// stop
Serial.println ("OK: finished and halted");
while (true)
delay (100);
}
以下図が測定結果をシリアルモニタからコピーして、グラフ化したもの。
確かに、analogRead()する度に、異なる値が読めているようですが、マクロに見てみると、370辺りの数値から始まって、徐々に数値が低下、280辺りで安定しているように見えます。安定しているとは言え、ミクロにみると、±4程度の差を持って、analogRead()する度に値がばらついているようです。
マクロには規則的で、ミクロにはノイズが観測できる、といった具合でしょうか？
この状態でリセットして、もう一度測定してみたところ、観測結果は似たようなものでした。今度は390あたりから始まって、280前後に収束します。大まかな挙動の再現性はある模様。
※フローティングノードのAD変換値と言う、不確定な事象を扱っているので、個体差があるかもしれません。皆様の機体では、違った挙動をするかもしれないので、ご注意ください。
今回真乱数の原料として使いたいのは、ミクロに見た時のノイズ成分です。±4なら、約2ビットのランダムな数値を、1回のanalogRead()で取得できるということなので、これを何回か、所望のビット数になるまで繰り返せば、必要なビット数の真乱数を生成することができるのではないでしょうか？
前節の測定を、対象のアナログ端子をGNDに落として実行すると、見事にall 0が観測できます。アナログ端子をVccに吊った時は、1,000個すべてが1,023と読み出されました。
「耐タンパ性については不問」なので、故意に対象のアナログ端子をGNDに落としたり、Vccに吊ったりするケースは無視します。しかし、誤ってそのような状態になってしまった時に、真乱数として、常に0x00000000や0xffffffffが生成されてしまう様では、困ってしまいます。
ノイズ源が期待通りに機能しているかどうか？を評価する手段が必要です。
もし、対象アナログ端子が、何らかの事故でGNDに落ちてしまった場合、K[i]＝0となるので、D[i]＝0であり、従って、B[i]＝0なので、ΣB[i]＝0＜Mとなり、エラーを検出できます。Vccに吊られたり、その他の電圧に固定された場合も同様です。
この評価を、真乱数の生成時にArduinoのプログラムでも実行して、真乱数の有効性(ΣB[i]≧M)を判定するものとします。
※そろそろ、くどくなってきますが、「耐タンパ性については不問」です。上の評価方法によって、確かに、固定値に対してはエラーを検出できますが、外部から都合の良い電圧波形を印加されると、生成する真乱数を制御することも、理論上は可能なことに注意してください。
ここまでに、analogRead()によって、1回当たり数ビットのノイズ成分(ランダム値)が観測可能なことがわかり、それを所望のビット数分求まるまで繰り返せば良さそうなことがわかりました。
測定したノイズ成分の「積み重ね方」が、本節のテーマです。
一番単純な方法は、以下図の様に、測定したノイズ成分1ビットを次々とシフトしながら積み重ねる方法です。
まあ、これでも良いでしょう……。ただ、今回の方法では、もう一工夫必要です。
前節の評価方法を思い出して、現在の測定値K[i]が、前回測定値K[i−1]と同じだった場合を考えましょう。この場合、D[i]＝0、従ってB[i]＝0となるので、最後のΣB[i]評価に対して、何ら、影響を与えません。上の例では、M＝4ビットであり、シフト可能な回数は4回に限られてしまうのですが、これでは、N＝4回繰り返した結果、ΣB[i]＜M＝4となるエラーケースが多発してしまいます。
そう、今回の手法では、求めるビット数Mの何倍かの回数Nを繰り返す必要があり、単なるシフトレジスタでは、この繰り返しに耐えられません。
そこで、以下図に示すような、LFSR(Linear Feedback Shift Register)というものを使います。
○に＋マークの素子は、EORゲートです。この構成なら、繰り返し回数がM(この例では＝4)回を上回ったとしても、最初に取り込んだノイズ成分も失うことなく、積み重ねることが可能です。
更に、ノイズ成分が多ビットだった場合を考慮して、以下図のように拡張します。
上図は、簡略化した、M＝4の場合の模式図ですが、これをM＝32に拡張した物をソフトウェアで実現して、使用します。
ここで、LFSRについて、簡単に特徴を説明しておきます。前節のLFSRの図にx4＋x3＋1と言う数式がありましたが、ガロア体(GF(2))の数式なので、「xの4乗と、xの3乗と、1を加算する」と、考えても、ダメです。筆者も、あんまし詳しくは無いので、以下より調べてみてください（１＾＾；）。
ガロア体について調べる(別ウインドウが開きます)
※スミマセン。エンジニアって、理屈を知らなくても、応用の仕方だけ知っていればOK、みたいなところがあるので……ｍ（１＿＿）ｍ。広く浅い知識でやらせていただいてます。
「x4＋x3＋1」のLFSRが、前節で示したような回路で実現できる、と言う所からスタートします。図には、4つのフリップフロップ(FF、1ビット記憶素子)がありますので、LFSRが保持するのは4ビットのデータです。「ノイズ」が常に0だったケースで、クロックを供給すると、LFSRの値がどのように変化するか、回路動作を追いかけてみましょう。
| クロック | LFSRの値 | 備考 |
| 1 | 0000 |
0000から始めると、何回クロックを入力しても、 0000のまま、変わらず。 |
| 2 | 0000 | |
| ： | ： | |
| 1 | 0001 | 0001からスタートすると、 |
| 2 | 1100 | 4ビットで表現可能な数値、 |
| 3 | 0110 | 1〜15までが満遍なく登場しながら |
| 4 | 0011 | 値が変化して行って、 |
| 5 | 1101 | ↓ |
| 6 | 1010 | ↓ |
| 7 | 0101 | ↓ |
| 8 | 1110 | ↓ |
| 9 | 0111 | ↓ |
| 10 | 1111 | ↓ |
| 11 | 1011 | ↓ |
| 12 | 1001 | ↓ |
| 13 | 1000 | ↓ |
| 14 | 0100 | ↓ |
| 15 | 0010 | 漸く最後に、 |
| 16 | 0001 | 元(0001)に戻った |
| 17 | ： | クロック2からの繰り返し |
LFSRと言うのは、上表のような動作をします。0からスタートすると、ずっと0のまま、と言う、マイナループを1つ持っていますが、それ以外の値からスタートすると、満遍なく、全ての値(1〜15)を出力した後、最初のスタート値に戻り、その動作を繰り返します。もうちょっと正確に言うと、そういう動作になるような数式「x4＋x3＋1」、すなわち回路、を選んで使います。
この4ビットLFSRの周期は、24−1＝15です(あっ、今度のは、普通に算数のべき乗)。8ビットのLFSRなら、例えば「x8＋x6＋x5＋x4＋1」が、周期28−1＝255となるLFSRになります。この数式、経験上、周期が最長となる式(これを原始多項式と呼んでいます)が知られているのですが、規則的に決まっているものではない様なので、未知のビット数だと、数式を見つけるのが一苦労です。
幸い、今回使いたいと思っている32ビットの数式は、いくつか知られているので、その中から以下の式を使います。
x32＋x22＋x2＋x1＋1
周期は232−1です。
※もし、上の多項式が原始多項式であることを確かめるには、0以外の値Xから始めて、次にXが登場するまでのクロック数が232−1クロックであること、を確かめるしかありません。ハードウェアならまだしも、ソフトウェア(論理シミュレータ含む)で確認するのは一苦労ですよ。
先ほどのシミュレーションは、「ノイズ」が0だった場合でした。次は、多ビット入力回路(下図、再掲)において、ノイズが1→2→3の順で3回出現した場合と、3→2→1の順で出現した場合でシミュレーションしてみます。
| クロック | ノイズ | LFSR値 | ノイズ | LFSR値 | 備考 | |||
| 1 | 0000 | 0001 | 0000 | 0001 |
ノイズが0000なら、 双方とも、同じ値に変化する。 |
|||
| 2 | 0000 | 1100 | 0000 | 1100 | ||||
| 3 | 0000 | 0110 | 0000 | 0110 | ||||
| 4 | 0001 | 1111 | 0011 | 1110 |
ノイズが発生 1→2→3の順(左)と、 3→2→1の順(右) |
|||
| 5 | 0010 | 1010 | 0010 | 0110 | ||||
| 6 | 0011 | 1000 | 0001 | 1110 | ||||
| 7 | 0000 | 0100 | 0000 | 1011 | ||||
| 8 | 0000 | 0010 | 0000 | 1001 | 最終結果は異なる。 |
ここで注目すべきなのは、最終結果が異なるということです。例えば、計測したノイズを、単純に加算した場合を考えましょう。計測した順が1→2→3でも、3→2→1でも、加算結果はどちらも6と、同じ結果しか残りません。しかし、LFSRを使うと、計測したデータそのもの、だけでなく、その発生順番も、状態に影響を与えるため、最終結果が異なるのです。
今回使う範囲だと、この位の情報があれば十分でしょう。
※実は、ノイズ成分を足し込んだ(EOR)事によって、演算結果がたまたま0000になってしまうと、マイナループに落ちてしまうので、その場合の処置をしておくべきなのかもしれませんね……。M＝32ビットにもなると、超々レアケースなので、無視してしまいました。
真乱数の生成コードを示します。
int idx;
int anlg;
int anlg_prev;
int dif;
uint32 trnd;
float dbit;
float sbit;
const int n_bits=8*sizeof (trnd);
const int n_iter=2*n_bits;
// generate
trnd=0;
(*nbit)=0;
anlg_prev=0;
sbit=0.;
for (idx=0; idx<=n_iter && sbit<n_bits; idx++)
{
// read noise
anlg=analogRead (ana_TRAND);
if (idx==0)
{
anlg_prev=anlg;
continue;
}
// lfsr: x^32+x^22+x^2+x^1+1
trnd ^= anlg;
trnd=((trnd>>1)&0x7fffffffUL) ^ ((trnd&0x01)? 0x80200003UL: 0);
// noise bits
dif=abs (anlg - anlg_prev);
dbit=(dif==0)? 0: log (dif)/log (2.);
if (dbit>=2.)
sbit += 2.;
else
sbit += dbit;
}
前節までに解説した通りですので、特に説明すべきことはありません。本文中で、MやNと言う記号を使っていたところが、n_bits、n_iterなどに変わっているので、その辺りを解読しながら読んでいただけば、理解可能だと思われます。
データとして取り込んだ、ノイズ源の評価もコード中で実施しているので、ループを抜けた後で、ΣB[i]＜Mとなったエラーケースを検知して、ロバスト処理を挟んであります。
// noise evaluation
(*nbit)=int (floor (sbit));
if (sbit<n_bits)
{
// ERROR: robustness
trnd ^= (uint32 (random (UINT16_MAX + 1UL))<<16) +
random (UINT16_MAX + 1UL);
}
期待したビット数Mに満たなかった場合、擬似乱数列random()の出力を足し込むことによって、見かけのランダム性を担保しています。
あくまで、ロバスト対応の処置です。真の乱数ではなくなっているので、本来は、エラーではじくべきです。
完成版のプログラムをリンクしておきます。
完成版プログラム：true_random.ino
コンパイルしてArduinoに転送すると、シリアルモニタから以下のコマンドを受け付けています。
実験サンプル：【 】内はシリアルモニタからの送信データ
【trnd,20】
OK+ true_random(), 1061433184, 58, 32.75, 11.296
OK+ true_random(), 1922511873, 25, 33.34, 6.212
OK+ true_random(), 3453499108, 25, 33.75, 6.864
ERR+ true_random(), 4265084378, 65, 14.34, 11.332
OK+ true_random(), 3751756738, 23, 33.58, 6.056
OK+ true_random(), 3399059579, 45, 32.75, 10.032
OK+ true_random(), 1319806150, 57, 32.43, 13.732
OK+ true_random(), 1161211969, 27, 33.75, 7.012
ERR+ true_random(), 3709696576, 65, 5.00, 11.244
OK+ true_random(), 147991446, 23, 32.17, 6.132
ERR+ true_random(), 2536686247, 65, 19.17, 11.472
OK+ true_random(), 92374271, 23, 33.75, 6.392
OK+ true_random(), 2604479380, 58, 33.75, 11.572
OK: trnd (err=3/10)
アナログ端子が壊れていたり、VccやGNDなど、固定電圧が印加されていたりすると、永遠にエラーが出続けて、制御が戻ってこないので、その場合は、あきらめて電源を切ってください。
そのような異常が無ければ、指定した個数の真乱数を生成して、表示します。表示データは、真乱数(uint32)、繰り返し回数、ΣB[i]、生成に要した時間[ms]、をカンマ(，)区切りで出力しています。
「ΣB[i]÷<繰り返し回数>」の計算により、analogRead() 1回分の平均ノイズビット数が求められます。
表示が"ERR+"で始まるデータは、繰り返し回数Nを超えてもΣB[i]≧Mとならなかった、エラーが発生したことを示しています。
生成したデータ(1,000個)をExcelに読み込んで、解析してみました。まずは、真乱数の発生頻度(ヒストグラム)です。
ヒストグラムと、右肩上がりの直線は累積%です。32ビットのデータが、一様に発生していることがわかります。気になる偏りは見られません。
次に、analogRead() 1回当たりの平均ノイズビット数ave(B[i])について、同様にヒストグラムを取ってみました。
変な形……。山が2か所ありますね。何で山が2つできるのか？は、わかりませんでした。しかし、1つ目の山が気になります。ave(B[i])≒0.45辺りにピークがあります。
この時、analogRead()の繰り返し回数N＝2×Mとしていました。真乱数のビット数M＝32ですので、N＝64です。ave(B[i])＝0.45だとすると、N回繰り返しても、ave(B[i])×N≒28.8にしかならず、M＝32ビットに届かないため、エラーと判定されてしまいます。
逆算すると、N＝64で、ΣB[i]≧Mとなるためのave(B[i])の条件は、ave(B[i])≧M/N＝0.5なので、上図のave(B[i])＜0.5となる発生頻度(累積%)を求めると、約22%でした。つまり、5回に1回以上の割合で、エラーが発生してしまいます。ちょっと多いですね。
min(ave(B[i]))≒0.36だったので、N＝M/min(ave(B[i]))≒2.8×Mより、N＝3×Mに設定することにしました。この結果、1,000回生成程度では、エラーの発生は認められなくなりました。ただし、確率の話になるので、何万回か繰り返すと、それなりの確率でエラーが発生します。生成時のエラーチェックとエラーへの対応(リトライなど)は、必ず必要となります。
※ここ、考え所で、繰り返し回数Nを大きくとると、回路に何らかの異常が生じた場合の、異常検知能力が低下してしまいます。本当は、N＝2×Mくらいに設定しておいて、適切にエラー処理(リトライなど)をした方が良いのかもしれません。
真乱数の生成にかかる時間ですが、analogRead() 1回には、約100μsの時間がかかります(「#4-4 デジタルオシロスコープ」参照)。それを最大N回繰り返すので、真乱数の生成には、そこそこ時間がかかります。1,000個の真乱数生成で要した1個当たり生成時間の平均値は約9.5ms、最小約5.4ms、最大約17.4msでした。
ただでさえ、フローティングノードのAD変換と言う、不確定な事象を扱うので、機体によってばらつきが見えてしまう恐れがあります。前節までの実験は、Arduino UNOで実施したのですが、手持ちの他の機体でも一通り真乱数を生成してみて、個体差が出ないか？検証してみました。
|
Arduino |
真乱数 一様性 |
繰り返し回数N | ave(B[i]) | 生成時間[ms] | |||||||
| MAX | MIN | AVE | MAX | MIN | AVE | MAX | MIN | AVE | |||
| UNO | ATmega328P-PU | ○ | 91 | 20 | 43 | 1.68 | 0.36 | 0.98 | 17.4 | 5.4 | 9.5 |
| Pro Mini | ATmega328P-AU | ○ | 59 | 20 | 33 | 1.61 | 0.57 | 1.09 | 14.3 | 5.7 | 8.3 |
| Micro | ATmega32U4-MU | ○ | 37 | 18 | 25 | 1.83 | 0.90 | 1.40 | 8.8 | 5.2 | 6.6 |
| STM32 | STM32F103C08 | ○ | 20 | 17 | 17 | 1.88 | 1.60 | 1.86 | 1.5 | 1.1 | 1.3 |
なるほど、そうですよね。正確には、「そりゃ、個体差は有るでしょ？これだけ不確定な事象を扱っているのだから」でした。個体差は有っても、結果として生成した真乱数が、それなりに使えそうなものかどうか？が大事ですね。
その点、肝心な真乱数の評価は、ヒストグラムによる一様性の確認だけなので、ちょっと不足かもしれません……。気になる方は、「擬似乱数列の評価方法」などで検索して、検証してみてください。
ここまでやっておいて言うのも何ですが……、AD変換器のノイズって、そもそも、MCUが動作したことに由来する電気的ノイズが主要因だったりします。ということは、MCUに全く同じ動作をさせると、毎回同じ様なノイズパタンが発生し、結果として同じ様な真乱数を生成してしまう可能性も、0ではないです。一応、そのようなケースに起因する「真乱数の偏りがないか？」をヒストグラムで確かめたつもりなので、心配いらないかもしれませんが、念のため、書き添えておきます。
