デジタル信号をもっとなまらせよう

## はじめに

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では満足できないんだ！」
そんなときは、2次LPFを使いましょう。
<img width="400" alt="2nd_lpf_gain.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/b409ebe1-c489-1b75-a23e-4f888c04e362.png">

## 微分方程式を立てる

ここでは、波形をなまらせるために、下図のような2次LPFを使います[^5]。
<img width="550" alt="2nd_lpf.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/54af0630-ca76-c66f-240d-a832f84407e6.png">
上図から、2次LPFの出力電圧が満たすべき式として、次の連立微分方程式が得られます[^1]。

[^5]: [Electronics Tutorials - Passive Low Pass Filter](https://www.electronics-tutorials.ws/filter/filter_2.html)
[^1]: 2次LPFの出力電圧は2階微分方程式で書き表せます。2階微分方程式を解くためには、初期条件が2つ必要になります。2次LPFにおいて初期条件としてすぐに思いつくのは、コンデンサの初期電圧がゼロという条件、すなわち $V_1(0)=0$ と $V_\mathrm{out}(0)=0$ の2つです。シミュレーション実行時に、これらの初期条件を使って計算を進めたいので、$V_1$ と $V_\mathrm{out}$ を状態変数にしています。


```math
\begin{cases}
\displaystyle \frac{dV_1}{dt}= -\frac{2}{RC}V_1 + \frac{1}{RC}V_\mathrm{out} + \frac{1}{RC}V_\mathrm{in}\\
\displaystyle \frac{dV_\mathrm{out}}{dt}= \frac{1}{RC}V_1 - \frac{1}{RC}V_\mathrm{out}
\end{cases}
```

状態方程式[^2]っぽく書くと、次のようになります。

[^2]: [Tutorialspoint: Control Systems - State Space Model](https://www.tutorialspoint.com/control_systems/control_systems_state_space_model.htm)

```math
\begin{align}
\begin{bmatrix}
\dot{V}_1\\
\dot{V}_\mathrm{out}\\
\end{bmatrix}
=
\frac{1}{RC}
\begin{bmatrix}
-2 & 1\\
1 & -1
\end{bmatrix}
\begin{bmatrix}
V_1\\
V_\mathrm{out}
\end{bmatrix}
+
\frac{1}{RC}
\begin{bmatrix}
1\\
0
\end{bmatrix}
V_\mathrm{in}
\end{align}
```

## 微分方程式の解析解を求める

先ほどの微分方程式は、大学の教科書に書かれているような解法を使えば、手計算で解くことができます……が、時間に追われる現代人にそんな余裕はないでしょう。そのため、ここでは、[Wolfram Cloud](https://www.wolframcloud.com/)を使って解を求めることにします[^3]。

[^3]: [Wolfram Cloud](https://www.wolfram.com/cloud/index.php.ja?source=footer)では、無料でWolfram言語が使えます。複雑な計算をおこないたい場合は、[WolframAlpha](https://www.wolframalpha.com/)よりも[Wolfram Cloud](https://www.wolframcloud.com/)のほうが便利です。

```
eqns ={
V1'[t]==(-2*V1[t]+Vout[t]+Vin)/tau,
Vout'[t]==(V1[t]-Vout[t])/tau,
V1[0]==V10,Vout[0]==Vout0
}
sol = DSolve[eqns, {V1,Vout},t]
Collect[(V1/.sol[[1]])[t],{Vin,V10,Vout0},Simplify]
Collect[(Vout/.sol[[1]])[t],{Vin,V10,Vout0},Simplify]
```

上記のコマンドを[Wolfram Cloud](https://www.wolframcloud.com/)上で実行すると、次の解析解が得られます。ただし、テキスト表示の都合上、$V_1(0)$ を`V10`、$V_\mathrm{out}(0)$ を`Vout0`と書いています。

2次LPFの内部電圧 $V_1(t)$
<img width="650" alt="2nd_lpf_v1_large.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/df08526d-3e9c-ec2a-c51c-c5615ec33823.png">
2次LPFの出力電圧 $V_\mathrm{out}(t)$
<img width="650" alt="2nd_lpf_vout_large.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/b26ea942-b1d2-91ca-b14e-79cdd318943b.png">
1次LPFの出力電圧の解析解はシンプルでしたが、2次LPFの場合は、ちょっと複雑です。この解析解をそのままモデルに書くのはしんどいので、考え方を変えましょう。今回は、「方程式の解ではなく、方程式そのものをモデルに組み込んで、シミュレーション実行時に逐次的に出力電圧を求めればいいや」と考えることにします。

## 微分方程式を離散時間で近似する

微小時間 $dt$ の代わりに、有限の短さの時間ステップ $\Delta t$ を使って、微分方程式を次のように近似します。

```math
\begin{cases}
\displaystyle \frac{\Delta V_1}{\Delta t}= -\frac{2}{RC}V_1 + \frac{1}{RC}V_\mathrm{out} + \frac{1}{RC}V_\mathrm{in}\\
\displaystyle \frac{\Delta V_\mathrm{out}}{\Delta t}= \frac{1}{RC}V_1 - \frac{1}{RC}V_\mathrm{out}
\end{cases}
```

上式の両辺に $\Delta t$ を掛けると、時間ステップごとの電圧の変化量 $\Delta V_1$ と $\Delta V_\mathrm{out}$ が得られます。これを利用して、$k$ 番目の時間ステップにおける内部電圧 $V_1[k]$ と出力電圧 $V_\mathrm{out}[k]$ を次のように表します[^4]。

[^4]: これは「前進オイラー法」と呼ばれる、微分方程式の数値解法のなかで最も素朴な方法です。精度が必要な場合は、少し式がややこしくなりますが、ルンゲクッタ法を使ったほうがいいと思います。なお、前進オイラー法、ルンゲクッタ法は、どちらも陽解法に分類されます。SPICEなどの回路シミュレータでは安定性の理由から、陽解法ではなく陰解法が使われます。しかしながら、「数値解の安定性が気になるようなきわどい回路解析をVerilog simulator上でやろう」なんてことを考える人はそうそういないはずなので、実装の簡単さを優先して陽解法を採用するのがいいんじゃないかな、と何となく思っています。

```math
\begin{cases}
\displaystyle V_1[k+1] = V_1[k] + \left(-\frac{2}{RC}V_1[k] + \frac{1}{RC}V_\mathrm{out}[k] + \frac{1}{RC}V_\mathrm{in}[k]\right) \Delta t\\
\displaystyle V_\mathrm{out}[k+1] = V_\mathrm{out}[k] + \left(\frac{1}{RC}V_1[k] - \frac{1}{RC}V_\mathrm{out}[k]\right) \Delta t
\end{cases}
```

## SystemVerilogで記述する

先ほど得られた差分方程式をもとに、SystemVerilogで2次LPFの動作モデルを記述します。

```verilog:LPF2.sv
module LPF2 #(
  parameter realtime timestep = 100ps,
  parameter real initial_vc1 = 0,
  parameter real initial_vc2 = 0,
  parameter real R = 1e3,
  parameter real C = 500e-15
)(
  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

  localparam realtime tau = R*C;
  realtime t;
  realtime last_t;
  realtime delta_t;
  real last_vin;
  real vc1;
  real vc2;
  real delta_vc1;
  real delta_vc2;

  always @(initial_step, next_step) begin
    if (initial_step.triggered()) begin
      vc1 = initial_vc1;
      vc2 = initial_vc2;
      last_vin = _vin;
      last_t = $realtime();
      -> start_timer;
    end else begin
      -> stop_timer;
      t = $realtime();
      delta_t = t - last_t;
      delta_vc1 = (-2*vc1 + vc2 + last_vin)*delta_t/tau;
      delta_vc2 = (vc1 - vc2)*delta_t/tau;
      vc1 += delta_vc1;
      vc2 += delta_vc2;
      last_vin = _vin;
      last_t = t;
      -> start_timer;
    end
  end

  assign VOUT = vc2 + 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/)

![2nd_lpf_wave.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/863735/f7110f37-a98b-2731-8277-de7e73472984.png)


## おわりに

最終回じゃないぞよ
もうちっとだけ続くんじゃ