白雜訊White Noise跟隨機行走Random Walk [FORTH語言]－早安，苦命工程師的胡言亂語

前言

最近因為工作上會用到，所以做了一些研究。當然啦，這篇偏向統計物理的心得筆記。但是因為有用 FORTH 來模擬了一下，應該也是不錯 FORTH 語言的參考，所以把它歸類在 FORTH 語言裡面。

一般而言，這種模擬可以用任何語言， R 語言， Python，MatLab ... 都行。因為 gFORTH 有很完整跟非常標準的浮點運算，所以我都把 gFORTH 當數值分析語言來用，也蠻方便的。做個心得筆記記錄，以後要用時查閱方便。

白雜訊

摘自維基百科，

「白雜訊，是一種功率譜密度為常數的隨機訊號或隨機過程。即此訊號在各個頻段上的功率一致。由於白光是由各種頻率（顏色）的單色光混合而成，因而此訊號的平坦功率譜性質稱為「白色」，此訊號也因此得名為白雜訊。相對的，其他不具有這一性質的噪音訊號則稱為有色雜訊。」

所以白雜訊，是一種隨機訊號，特色是在不同在不同頻率下具有相同的強度。這種聲音聽起來很平穩，有規律、一致性，呈均勻狀態。

進一步數學上，白雜訊的定義如下，其必須滿足三個條件

今有一 { Xi } 的隨機數列集合，

(1) E(Xi) = 0 其期望值為零

(2) Var(Xi) = sigma^2 其變異數為定值

(3) 任意的 Xi, Xj 無相關 Cov(Xi, Xj) = 0

這樣的隨機亂數的序列則被稱為白雜訊隨機序列。

white noise pic.png

電腦語言上的白雜訊

很顯然的，一般電腦語言上的虛擬亂數，大概就符合白雜訊的定義。

真的亂數是非決定性的，但是數位電腦是決定性的東西。怎樣的輸入，經過數位運算後就是怎樣的輸出。重複一樣的輸入，運算後就是重複一樣的結果。所以電腦亂數只能透過數學演算法模擬出一系列很接近真實亂數的序列。但是只要用來產生亂數的 seed 設為相同，這樣產出來的亂數序列也會一模一樣的。所以電腦語言上的亂數產生器被稱為虛擬亂數產生器。

為了模擬白雜訊，所以從 FORTH Scientific Libaray 上找到一個浮點亂數產生器，跟大多數電腦語言一樣，它可以產生 0 - 1 之間的一個浮點亂數。讓我們用它來測試模擬一些白雜訊的特性吧。

這個亂數產生器的使用方式很簡單，

(1) 首先你需要一個雙整數的 seed，可以任意選的雙整數如 1478.789 然後用 2! 存入 seed 雙整數變數。

(2) 然後用 prng 即可產生一個 0 - 1 之間的隨機浮點亂數至浮點堆疊之中供取用。

我們就用這個亂數產生器為基礎，來探索熟悉一下白雜訊的特性吧。

prng 會產生一個 0 到 1 之間的白雜訊，假定我們要的雜訊從 - max-step 到 + max-step 之間。所以

10e0 fconstant max-step

: rwalk ( -- F: pos) \ generate a random walk
max-step 2e0 f* prng 0.5e0 f- f*
;

rwalk 每執行一次，就會得到 - max-step 到 + max-step 之間的一個隨機雜訊。

就 prng 減去 0.5 ，讓它變成 -0.5 到 0.5 之間的亂數，再乘上 2*max-step 就可以得到 - max-step 到 + max-step 之間的亂數囉。

這裡取名叫 rwalk ，就像有一個醉漢，投一個 - max-step 到 max-step 的一個骰子，來決定這次是要往前或往後走多遠。這種叫做醉漢隨機行走 Random Walk。所以 rwalk 除了是每次白雜訊的雜訊外。也可以看做是醉漢走路，每次走路的步數跟方向。

如果我們把白雜訊對時間積分起來，就是把一系列n 個隨時間推進所產生的雜訊加起來。

或者是，醉漢隨機行走，一開始假設位置為0，走了n 次之後它的最終位置我們把它叫做 n walks 。所以

: walks ( n -- F: pos) \ generate many and walks
0e0 ( start-pos)
0
do rwalk f+
loop
;

就執行一個 n 次的迴圈，透過 rwalk 產生每次向前或向後的步數 (白雜訊)，然後把它們加起來回傳結果。

白雜訊的平均值跟變異數

接下來我們來算算醉漢隨機行走 Random Walk (或是白雜訊) 的平均值跟標準差。來看看它是否符合白雜訊的定義。

: E(rwalk) ( n -- F: avg.rwalk)
dup walks s>d d>f f/
;

平均值的計算，超簡單的，就用. n walks 叫它走 n 步後除上 n 的步數就是平均值 (期望值) 囉。

: Var(rwalk) ( n -- F: var.rwalk)
0Var!
0 do rwalk Var+! loop
Var@
;

變異數的計算，先用 0Var! 起始變異數計算。

然後用一個 n 次的迴圈，用 rwalk 產生一個隨機行走步數 (或白雜訊) 然後用 Var+! 輸入此資料，進行變異數計算。最後迴圈執行結束後，用 Var@ 取出變異數計算結果。

運算的指令

: Calc.rwalk ( #walks -- )
cr ." No. E(rwalk) Var(rwalk)" cr
10 0
do i 1+ 2 .r 3 spaces dup E(rwalk) ?space f. 2 spaces dup Var(rwalk) f. cr
loop drop
;

給定 #walks 醉漢隨機行走的次數(或白雜訊所產生的次數)

連續計算 10次，這樣狀況下每次隨機行走步數(或白雜訊)的平均值，跟變異數。測試結果如下

calc rwalk results.png

分別用 n=10, 100, 1000, 10000 來進行測試，

很明顯的，隨著 n 的增加，數據逐漸收斂

每次隨機行走或白雜訊的平均值趨近於零，而變異數趨近於 33.3 (標準差約 5.77) 是個定值。所以基本上，符合上面我們所述所謂白雜訊的定義 (1) 和 (2)。

變異數的計算

變異數數學上的定義

所以

所以對於一系列的數列資料，只要算出每個數列資料平方的平均值減去數列資料平均值的平方就是變異數囉。

variable NN \ no. of data

fvariable E(x^2)
fvariable E(x)

NN 用來記錄有多少資料已經輸入計算變異數了

E(x^2) 這個變數用來記錄目前資料平方的平均值

E(x) 這個變數用來記錄目前資料平均值

: 0Var! ( --) \ rest variance
0e0 E(x^2) f!
0e0 E(x) f!
0 NN !
;

0Var! 用來歸零這些計算變數

: Var+! ( F: data --) \ raw data input for variance
1 NN +!
NN @ fdup fdup f* E(x^2) f@ fswap avg E(x^2) f!
NN @ E(x) f@ fswap avg E(x) f!
;

Var+! 用來輸入最新資料，並更新 E(x^2) 跟 E(x)

就很簡單，將變數值取出來，給上目前的 NN 跟最新資料，利用 avg 算出新的平均值後 f! 回存結束。

: avg ( n F: avg data -- F: avg')
dup 0=
if fswap fdrop exit then

s>d d>f ( F: avg data n )
fswap fover ( F: avg n data n)
f/ ( F: avg n data/n)
frot frot ( F: data/n avg n)
fdup 1e0 f- ( F: data/n avg n n-1)
fswap f/ f*
f+
;

就給定 n 及舊 avg 值跟新 data 值，然後利用下面加權平均公式算出新 avg' 平均值

avg' = data/n + (n-1)/n * avg

最後

: Var@ ( F: -- variance) \ calculation variance
E(x^2) f@ E(x) f@ fdup f* f-
;

Var@ 利用 Var(x) = E(x^2) - E(x)^2 算出真的變異數出來

隨機行走最終位置，白雜訊的積分，陀螺儀的角度偏移

這裡，我們稍微拉回來一下白雜訊的積分，跟陀螺儀之間的關係。

來聊一下陀螺儀，

目前的陀螺儀都是微機電 MEMS 的陀螺儀，它們是一種 rate 陀螺。也就是說它們是一種角速度的感測器 (利用科氏力來偵測旋轉的角速度)，提供的是角速度。

所以如果你想要知道目前的移動角度，你必須對角速度做時間的積分才能得到角度的資訊。

偏偏，這個 MEMS 的角速度感測器會受到內部白雜訊訊號的干擾跟影響。如影隨行，各種頻率都有無法消除的白雜訊訊號，最後混進了角速度對時間的積分裡面去了，最後造成一種本質上由白雜訊所造成，無法消除的角度誤差。

這也很像醉漢隨機行走，最後我們發現醉漢最終位置一直在偏移改變，而形成一種無法預測的偏移。且這個偏移的變異程度會隨時間的增加而越來越擴大。

要玩陀螺儀的話，就要稍微了解一下這個由白雜訊訊號所造成的角度誤差的一些機制。

來定義一下，角度隨時間的誤差 S(t)，為白雜訊對時間的積分

Image 2-1-0.png

定義白雜訊的隨機變數為

所以這個積分可以變成

Image 2-1.png

不同時間的角度隨時間的誤差 S(t) 形成隨機變數序列，來求看看它的平均值跟變異數吧

平均值

Image 2-2.png

為零

變異數

Image 2-4.png

Image 2-5.png

Image 2-6.png

Bingo！我們找到陀螺儀奇特的特性了。

陀螺儀經過一段時間後，它的變異程度 (變異數) 會隨著時間的經過不斷的被放大。透過對陀螺儀角速度時間積分所獲得的角度，會因為白雜訊累積的關係，角度資訊會跳動的越來越厲害。

假如我們已經知道某一時刻下陀螺儀的標準差 sigma, 則之後所有的標準差如下之公式所表示，是時間的函數

Image 2-9.png

實際應用把第一小時的標準差稱作 ARW (Angle Random Walk) 角度隨機行走，這是陀螺儀的一個規格

透過剛剛的分析， ARW 的本質來自於白雜訊對時間的積分。其實它就是白雜訊的另外一種資訊的呈現。

當一個陀螺，它的 ARW 數據很差(數值很大)，其實就表示它的本質白雜訊很高，所以最後時間積分後角度變異很大。反之就是白雜訊低，所以時間積分後角度變異很小。

但是無論白雜訊的訊號有多低，陀螺經過長時間的使用後，它的角度變異數會隨時間的經過成正比。最後終究是會垮掉而不準確的。 😭

這就是陀螺儀目前技術的難處了，單純只靠陀螺儀本身，要長時間永久精確的知道物體角度的改變，就物理上而言，那是不可能的。一定要有第三者的另外一個感測器，提供另外一個數據例如目前的絕對角度，將已經發散的時間積分給拉回來，才能夠永久維持住陀螺儀姿態感測的準確度。

而當我們知道一個陀螺儀的 ARW 數據時，我們可以快速估計它不同的時間下，角度的標準差

例如，某牌的陀囉，它的 ARW = 0.11 deg/sqrt(h)

這代表

第 1 小時，它的角度標準差是 0.11 degree

第 2 小時，它的角度標準差擴散到 sqrt(2)*0.11 = 0.1556 degree

第 3 小時，它的角度標準差繼續擴散到 sqrt(3)*0.11 = 0.1905 degree

...

數值模擬

持續回到 FORTH 語言，簡單的測試一下啦

S(t) 就是 walks，來算一下對時間的積分

: E(walks) ( #try n -- F: avg.walks)
0e0 ( sum)
over 0
do dup walks f+
loop drop
s>d d>f f/
;

#try 就是計算的次數， n 就是 S(t) 裡面的 t

就讓它 walks 走 n 步，然後做 #try 次取平均

: Var(walks) ( #try n -- F: var.walks)
0Var!
swap 0
do dup walks Var+!
loop drop
Var@
;

#try 就是計算的次數， n 就是 S(t) 裡面的 t

就讓它 walks 走 n 步，然後做 #try 次計算變異數

計算指令

: Calc.walks ( #try #walks --)
cr ." No. E(walks) Var(walks)" cr
10 0
do i 1+ 2 .r 3 spaces 2dup E(walks) ?space f. 2 spaces 2dup Var(walks) f. cr
loop 2drop
;

最終結果

calc walks results 2.png

可以看到，隨著時間 n 的增加 10 倍，變異數也一直增加 10 倍！

(啊，這不是廢話嗎，太淺薄了啊！) 😜

FORTH 程式原始碼列表

這裡用了一個 FORTH 語言的浮點亂數函式庫，根據版權宣告，是需要列出作者的版權宣告的。這裡索性完整列出整個浮點亂數函式庫原始檔

\ Pseudo random number generator in ANS Forth
\
\ Leaves a pseudo random number in the range (0,1)
\ on fp stack.
\
\
\
\ ---------------------------------------------------
\ (c) Copyright 1998 Julian V. Noble. \
\ Permission is granted by the author to \
\ use this software for any application pro- \
\ vided this copyright notice is preserved. \
\ ---------------------------------------------------
\
\ Based on GGUBS algorithm: s’ = 16807*s mod (2^31-1)
\ P. Bratley, B.L. Fox and L.E. Schrage, A guide to simulation
\ (Springer, Berlin, 1983).
\
\ To simplify transport to 16-bit machines the 32-bit
\ modular division is performed by synthetic division:
\
\ Let b = d * m1 + m2 (= 2^31 - 1)
\
\ so that ( [n] means "largest integer <= n" )
\
\
\ s’ = (s*m1) MOD b = s*m1 - [s*m1/b]*b
\ = m1 * (s - [s/d]*d) - m2 * [s/d]
\ Environmental dependence:
\ 1. assumes at least 32-bit DOUBLEs
\ 2. needs FLOATING and DOUBLE wordsets
\ 3. assumes separate floating point stack
\

-rand
MARKER -rand

2VARIABLE seed

21474.83647 D>F FCONSTANT bigdiv \ 2^31-1
1277.73 D>F FCONSTANT divis
16807. D>F FCONSTANT m1
2836. D>F FCONSTANT m2

: (rand) ( adr --) ( f: -- seed')
DUP 2@ D>F divis FOVER FOVER ( f: s d s d)
F/ F>D 2DUP D>F ( [s/d]) ( f: s d [s/d])
F* F- ( f: s-d*[s/d])
m1 F* ( f: m1*[s mod d])
D>F m2 F* F- FDUP F>D ( adr d) ROT 2! ;

: prng ( f: -- random#)
seed (rand) bigdiv
FSWAP FDUP F0< ( -- f) ( f: -- 2**31-1 seed)
IF FOVER F+ THEN FSWAP F/ ;

: test 0.1 seed 2! 1000 0 DO prng FDROP LOOP seed 2@ D. ;
\ TEST 522329230 ok

再來才是我們的程式碼如下

\
\ Random Walk Simulation
\ Frank Lin 2022.2.19
\

--rand--
marker --rand--