search
LoginSignup
1
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

SystemVerilog Advent Calendar 2020 Day 5

posted at

updated at

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

はじめに

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

  1. デジタル信号をなまらせよう
  2. デジタル信号をもっとなまらせよう
  3. デジタル信号をもっともっとなまらせよう
  4. デジタル信号にノイズをのせよう ← 今ここ

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

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

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

差分方程式を立てる

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

\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}$ の差分を表す式が得られます。

\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}

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

\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)$ に従う乱数です。これは確率微分方程式2の数値解法の中でも最も素朴な方法であり、オイラー丸山法3と呼ばれています。

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

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

ここでは、次式で表されるボックスミューラー法を使って、正規分布に従う real 型の乱数を生成します。

\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}

上式をSystemVerilogで実装すると、次のようになります5

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

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

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

SystemVerilogで記述する

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

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

シミュレーション結果

おわりに

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


  1. デートDV110番 - デートDVについて 

  2. 「確率微分方程式」という単語に聞き覚えがない場合は、AIcia Solid Project - 確率微分方程式 がおすすめ。厳密な話をおいといて、ノリと感覚で説明してくれるので、何となく分かった気になれます。たぶん。 

  3. mutsumunemitsu - 確率微分方程式を数値計算してみよう
    Mustafa Bayram - Numerical methods for simulation of stochastic differential equations
    Numerical SDE Simulation - Euler vs Milstein Methods
    Civil Engineering - Random Vibration 

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

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
1
Help us understand the problem. What are the problem?