この記事はD言語 Advent Calendar 2019の17日目の記事です （LDPC符号の記事を書く予定でしたが，私の理解不足により，BPSK/QPSKのお話に変わりました）.
（あと，まだ書きかけなので．．．許してください）
様々な分野で共通だと思いますが，なにかの技術の性能や特徴を調べるためには，実際に機器を使用して実験をしてみるのが一番です．
ですが，実機実験が不可能な場合というのはよく起こりえます． たとえば，高価な装置が必要な場合や，実機実験のために膨大な労力が必要になるなど，実験ができない状況はいろいろ考えられます． そのようなときにシミュレーションや理論解析が効力を発揮します．
これは主観的な意見になりますが，通信工学は他の工学分野と比較して，理論解析とシミュレーションがよく浸透している分野だと思います． それは，数学・物理学の一分野として発展した電気回路や電磁気学の影響を大きく受けつつ工学に昇華したものが通信工学であるため，だと思います．
通信工学の中でも，さらに実験が難しい分野として無線通信があります． なぜ難しいかというと，まず電波法に制限されています． たとえば，実験目的で無線の送受信機を作成したとしても，特定実験試験局という免許が必要になり，手軽に実験できるわけではありません． （東海地区の特定実験試験局の免許一覧はここから見れますが，うちの研究室は14局もありますね）
次に，実験系を構築するのが難しいというのがあります． たしかに，FPGAやDSP, CPUで信号処理を記述するだけで無線機が作れるソフトウェア無線があることで実機実験はやりやすくなったのだと思います． でも，ソフトウェア無線では無理なこと，たとえばハードウェアである高周波回路を変更する必要がある場合には高周波回路のエキスパートになる必要があります．
他にもいろいろ理由はあると思いますが，これらの理由から無線通信の分野ではよくシミュレーションや理論解析で評価をします．
そして，僕が面白いなあと思うことは，現在の無線機は理論やシミュレーションモデルとよく一致するように回路が作られているという点です． つまり，「無線機ではどのような回路が最も望まれるか？」という答えには「理論結果と一致する回路」になるということです． 後ほど説明しますが，無線通信では複素数を使って信号を数学的に記述しますが，実際の無線機の高周波回路では複素数の信号を実数の高周波信号に変換する回路が入っています！
さて，話が長くなりましたので，そろそろ本題に入りますね．
等価低域系（ベースバンド）信号
たとえば広く一般に利用されている無線通信であるWi-Fiでは2.4GHzとか5.2GHzの周波数を利用して情報をやり取りします． このときに電波として飛んでいる高周波の信号$\hat{x}(t)$は，キャリア周波数$f_c(=2.4\times 10^9, 5.2\times 10^9)$と振幅変化$A(t)$や位相変化$\theta(t)$を用いて次のような形で記述されます．
\[\begin{align} \hat{x}(t) &= A(t) \cos(2\pi f_c t + \theta(t)) \end{align}\]この周波数帯の信号のシミュレーションをするには，サンプリング定理より5GHz～10GHzのサンプリングレートが必要になります． つまり，1秒間の信号をシミュレーションするのに最低でも$5 \times 10^9$サンプルも処理する必要があるわけです． たった1秒のシミュレーションでこんなに処理していては，いつまで経ってもシミュレーションが終わりません．
でも安心してください． 通信工学では（もっというと交流回路理論も），このような高周波帯域（パスバンド）の信号を扱わずに信号をモデル化する方法が古くから利用されています． それが等価低域系信号（等価ベースバンド信号）です．
さきほどのパスバンド信号に対応するベースバンド信号は次のように記述されます． ここで，$j$は虚数単位です．
\[\begin{align} x(t) &= A(t) \cos(\theta(t)) + j A(t) \sin(\theta(t)) \\ &= A(t) e^{j \theta(t)} \end{align}\]パスバンド信号に比べてベースバンド信号は複素数になりましたが，非常に簡潔に記述されていることがわかると思います． さらに，周波数に着目すると，ベースバンド信号には$2\pi f_c$が無いのです． つまり，中心周波数は0Hzです！ もちろん，$A(t)$も$\theta(t)$も時間変化するため，サンプリング周波数を0Hzにすることはできませんが，サンプリング周波数をパスバンド信号の10分の1から100分の1程度まで減らすことができます！ これは，シミュレーション時間が10分の1や100分の1になることを意味します．
このような，実際の信号は実数なのに数式モデル上では複素数を用いるというテクニック（？）は，交流回路理論や制御工学などでもよく見ると思います．
簡潔さやサンプリングレート以外にも，パスバンド信号よりも扱いやすい特徴がベースバンド信号にはまだあります． たとえば簡単のために$A(t)=A$，$\theta(t)=\theta$と時間固定したときに，遅延した信号をパスバンドで考えると，
\[\begin{align} \hat{y}(t) &= \hat{x}(t - \tau) = A \cos(2 \pi f_c (t - \tau) + \theta) \\ &= \cos(2 \pi f_c \tau) \hat{x}(t) + A \sin(2\pi f_c \tau) \sin(2 \pi f_c t + \theta) \end{align}\]というように，元の信号$\hat{x}(t)$をある係数$\cos(2 \pi f_c \tau)$倍した信号と，元の信号$\hat{x}(t)$には含まれていない$\sin$成分が出現します（この$\sin$成分は元の$\hat{x}(t)$のヒルベルト変換になっています）．
一方で，ベースバンド信号ですと，以下のように単に係数$e^{j2\pi f_c \tau}$を乗じるだけで遅延が表現できます．
\[\begin{align} y(t) = x(t - \tau) = x(t) e^{j2\pi f_c \tau} \end{align}\]このような性質を持つため，理論解析やモデルの記述には，実際に回路や空間を飛んでいる波形であるパスバンド信号ではなく，等価低域系信号が使われており，「解析信号」とも呼ばれています．
では，ベースバンド信号からパスバンド信号を生成するにはどうすればいいのでしょうか． 実はそれほど難しくありません． パスバンド信号とベースバンド信号には以下の関係があります．
\[\hat{x}(t) = \mathrm{Re}[x(t)]\cos(2\pi f_c t) - \mathrm{Im}[x(t)]\sin(2\pi f_c t)\]そして，皆さんが持っているスマートフォンでは，この数学的な変換がI/Qミキサという高周波回路によって達成されています！ つまり，現実上の無線機でもディジタル信号処理部分ではベースバンド信号を使うようになっていますが，実際に無線信号として送受信する場合にはI/Qミキサによってパスバンド信号に変換して伝送するのです．
このように，ベースバンド信号は理論解析やシミュレーションのみならず，実際の無線機でも使われています． そして，ベースバンド信号からパスバンド信号への変換は高周波回路がやってくれるので，私達はベースバンド信号のみに注目すればよいのです．
通信路の等価低域系表現
全く歪みのなく，雑音も無しに，送った信号がそのまま受信されるような通信路では，送信信号$x(t)$と受信信号$y(t)$は次のような関係になります．
\[y(t) = x(t)\]$y(t)$と$x(t)$はそれぞれ複素数であることを忘れないでくださいね？
さて，通信路で雑音を受ける状況を考えます． 通信ではよく加法性白色ガウス雑音（Additive White Gaussian Noise: AWGN）通信路という通信路を考えます． このAWGN通信路における送信信号$x(t)$と受信信号$y(t)$の関係は次のようになります．
\[y(t) = x(t) + n(t)\]ここで，$n(t)$は平均0で電力がある値（ここでは$N$)の複素ガウス分布$\mathcal{CN}(0,N)$に従う雑音です．
信号帯雑音比（Signal-to-Noise ratio: SNR）という言葉はいろいろな分野で出現しますが，この通信路のSN比は次のようになります．
\[\mathrm{SNR} = \frac{\mathbb{E}[|x(t)|^2]}{N}\]分子は信号の平均電力，分母は雑音の平均電力です．
では，電力$N$の雑音を付加するAWGN通信路をD言語で作ってみましょう．
まず，複素ガウス分布に従う雑音を生成しますが，ここではよく知られたBox-Muller法を利用します． Box-Muller法では，$(0, 1)$の一様分布に従う確率変数$X, Y$に対して次の変換を加えることで正規分布$\mathcal{N}(0, 1/2)$に従う独立な確率変数$U, V$を生成します．
\[\begin{align} U &= \sqrt{-\log X} \cos(2 \pi Y) \\ V &= \sqrt{-\log X} \sin(2 \pi Y) \end{align}\]$Z = U + jV$とすれば$Z$は標準複素正規分布$\mathcal{CN}(0,1)$に従う確率変数ですから，
\[Z = \sqrt{-\log X} e^{j2\pi Y}\]という変換によって複素数バージョンのBox-Muller法になります（振幅方向にレイリー分布，位相方向に一様分布なことがよくわかりますね）． これをプログラムにすると次のようになります． 単体テストには，プログラムが正しいかどうかを確認するために，実部と虚部の相関行列を計算して，それぞれの電力が0.5であることと，実部と虚部が無相関であることを確認しています．
module complex_gaussian;
import std.complex;
import std.math;
import std.random;
/// 平均0, 電力1の複素ガウス分布する確率変数の乱数
Complex!F complexGaussian01(F, UniformRNG)(ref UniformRNG rnd)
{
F x = uniform01!F(rnd);
F y = uniform01!F(rnd);
return cast(Complex!F)(sqrt(-log(x)) * std.complex.expi(2 * PI * y));
}
unittest
{
import std;
immutable N = 1024 * 1024;
auto rnd = Xorshift(0);
auto rcoef =
generate!(() => complexGaussian01!double(rnd)).take(N)
.map!("a.re * a.re", "a.re * a.im", "a.im * a.im")
.fold!("a + b[0]", "a + b[1]", "a + b[2]")(0.0, 0.0, 0.0)
.tupleof.only.map!(a => a / N);
assert(rcoef[0].approxEqual(0.5, 0, 1E-2)); // 実部の電力
assert(rcoef[1].approxEqual(0.0, 0, 1E-2)); // 実部と虚部の相関
assert(rcoef[2].approxEqual(0.5, 0, 1E-2)); // 虚部の電力
}