デジタル信号にノイズをのせよう

## はじめに

Verilogのシミュレータで表示されるデジタル信号って、きれいな形してますよね。現実世界では、シリコンもPCB上も職場も家も、ノイズばかりなのに。アナログ回路におぜん立てされた環境で、ぬくぬくとしているデジタル信号にも、私たちと同じ苦しみを味わってもらいましょう。

1. [デジタル信号をなまらせよう](https://qiita.com/ds54e/private/f821bb6010355d6cb020)
2. [デジタル信号をもっとなまらせよう](https://qiita.com/ds54e/private/a56d97306a7a08cbec70)
3. [デジタル信号をもっともっとなまらせよう](https://qiita.com/ds54e/private/7a13da3d7586f1901a7c)
4. [デジタル信号にノイズをのせよう](https://qiita.com/ds54e/private/f2b61a7ce50397087137) ← 今ここ

## ノイズを入れつつ、波形もなまらせる

先ほどは、ついかっとなって、デジタル信号にノイズを食らわせてやるなんて、ひどいことを言ってしまいました。しかし、私は知っています、憎しみは何も生みません。悲しみの連鎖はここで終わらせるべきなのです。

……ということで、ノイズをぶち込んで気持ちがスッキリしたら、自分がやったことを反省し、デジタル信号にのせてしまったノイズを1次LPFでやさしく取り除いてあげましょう。「ごめん、あの時はどうかしていた」「もう絶対しないから許してほしい」と、相手に気持ちを伝えながら、ね[^dv]。

![add_noise.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/6400a758-97f2-c03c-e956-9991c96230d6.png)


[^dv]: [デートDV110番 - デートDVについて](https://ddv110.org/aboutddv)

## 差分方程式を立てる

正規分布 $N(0, \sigma^2)$ に従うホワイトノイズ $\xi(t)$ を加えた場合、微分方程式は次のように書けます。

```math
\begin{align}
&V_\mathrm{out}(t) + RC\frac{dV_\mathrm{out}(t)}{dt}
= V_\mathrm{in}(t) + \xi(t)\\
\therefore&
dV_\mathrm{out}
= \left(-\frac{1}{RC}V_\mathrm{out} + \frac{1}{RC}V_\mathrm{in}\right)dt
+ \frac{1}{RC} dW
\end{align}
```

つづいて、微小時間 $dt$ を有限の短さの時間ステップ $Δt$ で近似します。すると、次のような、時間ステップごとの出力電圧 $V_\mathrm{out}$ の差分を表す式が得られます。

```math
\begin{align}
\Delta V_\mathrm{out}
&= \left(-\frac{1}{RC}V_\mathrm{out} + \frac{1}{RC}V_\mathrm{in}\right) \Delta t
+ \frac{1}{RC} \Delta W
\end{align}
```

上の式を使うと、出力電圧は最終的に次の差分方程式で表せます。

```math
\begin{align}
V_\mathrm{out}[k+1]
&= V_\mathrm{out}[k]
+ \left(-\frac{1}{RC}V_\mathrm{out}[k]
+ \frac{1}{RC}V_\mathrm{in}[k]\right) \Delta t
+ \frac{1}{RC} \Delta W[k]\\
+ &= V_\mathrm{out}[k]
+ \left(-\frac{1}{RC}V_\mathrm{out}[k]
+ \frac{1}{RC}V_\mathrm{in}[k]\right) \Delta t
+ q[k]
\end{align}
```

ここで、$q[k]$ は正規分布 $N \left(0, \sigma^2\frac{\Delta t}{RC}\right)$ に従う乱数です。これは確率微分方程式[^sde]の数値解法の中でも最も素朴な方法であり、オイラー丸山法[^em]と呼ばれています。

[^sde]: 「確率微分方程式」という単語に聞き覚えがない場合は、[AIcia Solid Project - 確率微分方程式](https://www.youtube.com/playlist?list=PLhDAH9aTfnxIhf-iRKYTVOSXPqDGgfRFP) がおすすめ。厳密な話をおいといて、ノリと感覚で説明してくれるので、何となく分かった気になれます。たぶん。

[^em]: [mutsumunemitsu - 確率微分方程式を数値計算してみよう](https://qiita.com/mutsumunemitsu/items/4d71d8a44d79255bd9ba)
 [Mustafa Bayram - Numerical methods for simulation of stochastic differential equations](https://advancesindifferenceequations.springeropen.com/articles/10.1186/s13662-018-1466-5)
 [Numerical SDE Simulation - Euler vs Milstein Methods](https://hautahi.com/sde_simulation)
 [Civil Engineering - Random Vibration](http://www.iitg.ac.in/budhaditya.hazra/CE%20607/Into_SDEQ.pdf)

## ボックスミューラー法を使う

ホワイトノイズを模擬するために、正規分布に従う乱数が必要です。しかし、SystemVerilogには、正規分布に従う `real` 型の乱数を生成する関数が用意されていません[^dist_normal]。

ここでは、次式で表される[ボックスミューラー法](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform)を使って、正規分布に従う `real` 型の乱数を生成します。

```math
\begin{align}
&Z_0 = \sqrt{-2 \ln U_1} \cos(2\pi U_2)\\
&Z_1 = \sqrt{-2 \ln U_1} \sin(2\pi U_2)
\end{align}
```

[^dist_normal]: SystemVerilogには、正規分布に従う乱数を生成する関数 `$dist_normal()` があります……が、実数ではなく整数を返すので、今回の用途には使えませんでした。

上式をSystemVerilogで実装すると、次のようになります[^dpic]。

[^dpic]: 今回は、ボックスミューラー法をSystemVerilogで直接実装しましたが、C言語で書いてDPI-Cで呼び出すというやり方もあるようです。[HDL Advent Calendar 2015](https://qiita.com/advent-calendar/2015/hdl)で紹介されています（今から5年前！）
[mizutomo - SystemVerilogって正規分布関数ないんだってよ](https://qiita.com/mizutomo/items/b4d368e3544cc8a39b64)

```verilog:rdist_normal.sv
module test;
  timeunit 1s;
  timeprecision 1ps;
  
  function real rdist_normal (
    input real mu,
    input real sigma
  );
    localparam real pi = 4*$atan(1);
    localparam real uint_max = $itor(2)**32 - 1;
    real u1, u2;
    real z0;
    u1 = $urandom()/uint_max;
    u2 = $urandom()/uint_max;
    z0 = mu + sigma * $sqrt(-2*$ln(u1)) * $cos(2*pi*u2);
    return z0;
  endfunction
  
  parameter int N = 200000;
  int fd;
  
  initial begin 
    fd = $fopen("rdist_normal.csv", "w");
    for (int i=0; i<N; i++) begin
      $fdisplay(fd, rdist_normal(0, 1));
    end
    $fclose(fd);
    $finish();
  end

endmodule
```

上のテストベンチで得られた乱数の確率分布を、正規分布と比較します。

```py:rdist_normal.py
import os
import pathlib
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

os.chdir(pathlib.Path("__file__").parent)
z = np.loadtxt("rdist_normal.csv")
plt.hist(
  z,
  bins=50,
  density=True,
  ec='black',
  lw=0.5,
  alpha=0.8,
  zorder=2,
  label="rdist_normal",
)

x = np.linspace(-3, 3, 301)
y = norm.pdf(x, 0, 1)
plt.plot(
  x, y,
  lw=3,
  zorder=3,
  label="N(0,1)",
)

plt.grid()
plt.legend()
plt.xlim([-3, 3])
plt.savefig("rdist_normal.png")

```
最終的に下図のような結果が得られます。

![rdist_normal.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/d046b6f9-5621-3814-6f45-12fa11d534a2.png)

# SystemVerilogで記述する

先ほど得られた差分方程式と乱数を使い、SystemVerilogで動作モデルを記述します。

```verilog:LPF1N.sv
module LPF1N #(
  parameter realtime timestep = 100ps,
  parameter real initial_vc = 0,
  parameter real R = 1e3,
  parameter real C = 500e-15,
  parameter real sigma = 0.1
)(
  input real VIN,
  input real VSS,
  output real VOUT
);
  timeunit 1s;
  timeprecision 1ps;

  real _vin;
  assign _vin = VIN - VSS;

  event initial_step;
  initial begin
    #0;
    -> initial_step;
  end

  event start_timer;
  event stop_timer;
  event next_step;

  always @(start_timer, stop_timer) begin
    if (start_timer.triggered()) begin
      fork begin
        #(timestep);
        -> next_step;
      end join_none
    end else begin // stop_timer
      disable fork;
    end
  end

  function real rdist_normal (
    input real mu,
    input real sigma
  );
    localparam real pi = 4*$atan(1);
    localparam real uint_max = $itor(2)**32 - 1;
    real u1, u2;
    real z0;
    u1 = $urandom()/uint_max;
    u2 = $urandom()/uint_max;
    z0 = mu + sigma * $sqrt(-2*$ln(u1)) * $cos(2*pi*u2);
    return z0;
  endfunction

  localparam realtime tau = R*C;
  realtime t;
  realtime last_t;
  realtime delta_t;
  real last_vin;
  real vc;
  real vn;
  real delta_vc;
  
  always @(initial_step, next_step) begin
    if (initial_step.triggered()) begin
      vc = initial_vc;
      last_vin = _vin;
      last_t = $realtime();
      -> start_timer;
    end else begin
      -> stop_timer;
      t = $realtime();
      delta_t = t - last_t;
      delta_vc = (-vc + last_vin)*delta_t/tau;
      vc += delta_vc;
      vn = rdist_normal(0, sigma*$sqrt(delta_t/tau));
      last_vin = _vin;
      last_t = t;
      -> start_timer;
    end
  end

  assign VOUT = vc + vn + VSS;

endmodule
```

## シミュレーション結果

* Simulator: [ModelSim ALTERA STARTER EDITION vsim 10.5b Simulator 2016.10 Oct  5 2016](https://www.intel.com/content/www/jp/ja/programmable/downloads/download-center.html)
* Viewer: [GTKWave Analyzer v3.3.71](http://gtkwave.sourceforge.net/)

![add_noise_wave.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/d15715e0-f61d-e8dd-09c3-392c3fad04f1.png)

## おわりに

今日もおつかれさまなのじゃ
くたくたいふくに、やつれたおかお
よしよし、わらわにおまかせじゃ、そうじにせんたく、あったかごはん[^senkosan]

[^senkosan]: [もっふもふ DE よいのじゃよ - 世話やきキツネの仙狐さん](https://mora.jp/package/43000033/A72849/)