走行中の車でラジオを聞いていると、信号待ちのほんの数メートルの移動で音が急にざらついたり、逆にクリアに戻ったりした経験はないでしょうか。あるいはスマートフォンで動画を見ながら歩いていると、立ち止まっているときより歩いているときの方が、通信速度の変動が激しく感じられることがあります。これらはどちらも、受信機が動くことで電波の受かり方が時々刻々と変わる「フェージング」という現象です。とりわけ厄介なのは、移動の速さに応じて変動の「速さ」そのものが変わる点です。歩く程度ならゆっくりとした変動ですが、新幹線の速度になると変動はずっと速くなります。
この「移動速度に応じて統計的性質が決まる時間変動するチャネル」を、計算機の中で忠実に再現したい——これがJakesモデルの動機です。実機で電波暗室や走行試験を繰り返すのは費用も時間もかかります。代わりに、現実のフェージングと同じ統計（包絡線の分布、変動の速さ、時間相関）を持つ複素利得の時系列をソフトウェアで生成できれば、受信機の設計やビット誤り率の評価が机上で完結します。Jakesモデルは、この要求に対して「少数の正弦波を足し合わせる」という驚くほど単純な仕掛けで応えた、移動体通信の古典にして今なお現役の手法です。
この図がJakesモデルの全体像です。移動する受信機の周囲には全方位に散乱体があり、各到来波は到来角 $\theta$ に応じて $f=f_d\cos\theta$ だけ周波数がずれます。前方（$\theta\approx0$）からの波は $+f_d$、後方（$\theta\approx\pi$）からは $-f_d$、真横からはシフトゼロです。到来角が一様に分布するという仮定から、後で見るU字型スペクトルと $J_0$ 相関がすべて導かれます。
Jakesモデルを理解すると、以下のような場面で役立ちます。
- 無線受信機の性能評価: LTE/5Gの端末やWi-Fi機器のビット誤り率（BER）を、走行速度を変えながらシミュレーションで検証できます。フェージングに強い等化器や符号化方式の比較が机上でできます。
- ダイバーシティ・MIMOの設計: アンテナ間隔やインターリーバ長を決める際、フェージングの時間相関（どれくらいの時間で利得が無相関になるか）が設計指針になります。Jakesモデルが与える相関関数が直接効いてきます。
本記事の内容
- 散乱体の到来角が一様分布することからClassical（U字）ドップラースペクトルを導出する
- ドップラースペクトルの逆フーリエ変換が、自己相関 $J_0(2\pi f_d \tau)$ になることを示す
- Jakesの正弦波和（sum-of-sinusoids）で複素利得の時系列を生成する手順
- 生成した利得の包絡線がレイリー分布、位相が一様分布になることをPythonで確認する
- 自己相関がベッセル関数 $J_0$ に一致することを数値的に実証する
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- マルチパス伝搬とは — 多重波がもたらす干渉とフェージング — 複数の反射波が干渉してフェージングが生じる仕組み
- 仲上・ライスフェージングチャネルの理論 — 見通し波がある場合のフェージング。レイリーはその特別な場合です
- フェージングとダイバーシティ — 変動の影響と対策 — フェージングへの対策技術と、相関の役割
フェージングはなぜ起きるのか — 多重波の干渉という出発点
まず、なぜ受信機が少し動くだけで信号強度が大きく変わるのか、その物理的なイメージを固めておきましょう。基地局から放たれた電波は、建物・地面・車などさまざまな物体で反射・散乱され、たくさんの経路を通って受信機に届きます。これがマルチパス伝搬です。それぞれの経路は長さが違うので、受信機に届くときの位相がバラバラです。同じ周波数の波が位相を変えて重なると、ある場所では山と山が重なって強め合い、別の場所では山と谷が重なって打ち消し合います。
ここで受信機が動くと、各経路の長さがミリ単位で変わり、位相関係が刻一刻と変化します。その結果、強め合っていた点が打ち消し合う点に変わったり、その逆が起きたりします。これが、移動に伴って受信レベルが激しく上下する正体です。電波の波長は、たとえば2 GHzで約15 cmです。受信機が波長の半分（7.5 cm）ほど動けば、位相は大きく変わり、信号レベルは「深いフェージング（ディープフェード）」と呼ばれる急激な落ち込みを起こすことがあります。
重要なのは、受信機の「移動」が時間軸での変動に直結する点です。空間的に7.5 cmごとに変動の節があるとして、その上を速度 $v$ で移動すれば、時間軸では一定のリズムで変動が現れます。このリズムの速さを定量化するのが、次に説明するドップラー周波数です。受信機の動きが各到来波の周波数をわずかにずらし、その「ずれの広がり」がフェージングの時間スケールを決めます。では、その広がりの形——ドップラースペクトル——がどう決まるのかを見ていきましょう。
ドップラー効果と到来角 — スペクトルが広がる仕組み
1本の到来波のドップラーシフト
救急車のサイレンが近づくと高く、遠ざかると低く聞こえる——これがドップラー効果です。電波でも同じことが起こります。受信機が速度 $v$ で動いているとき、ある方向から到来する電波の見かけの周波数は、受信機がその波源に近づくか遠ざかるかでずれます。
受信機の進行方向と、ある到来波の到来方向（電波が来る向き）のなす角を $\theta$ とします。このとき、その波が受ける周波数のずれ（ドップラーシフト）は次式で与えられます。
$$ f = f_d \cos\theta, \qquad f_d = \frac{v}{\lambda} = \frac{v f_c}{c} $$
ここで $f_d$ は最大ドップラー周波数で、$v$ は移動速度、$\lambda$ は波長、$f_c$ は搬送波周波数、$c$ は光速です。$\theta = 0$（受信機がまっすぐ波源に向かって進む）なら $f = +f_d$ と最大、$\theta = \pi$（まっすぐ遠ざかる）なら $f = -f_d$ と最小、$\theta = \pi/2$（真横から来る波）ならシフトはゼロです。
具体的な数値感覚をつかみましょう。$f_c = 2\,\mathrm{GHz}$、時速72 km（$v = 20\,\mathrm{m/s}$）なら、$f_d = 20 \times 2\times 10^9 / (3\times 10^8) \approx 133\,\mathrm{Hz}$ です。各到来波の周波数は搬送波の周りに $\pm 133\,\mathrm{Hz}$ の範囲でばらつくことになります。この「ばらつきの範囲」がドップラー拡がりです。
到来角が一様分布する理由
問題は、たくさんの到来波の周波数シフトが、$-f_d$ から $+f_d$ の間でどのように分布するか、です。これを決めるのが、到来角 $\theta$ の確率分布です。
Jakesモデル（より正確にはClarkeのモデル）では、受信機の周りの散乱体が水平面内であらゆる方向に均等に存在すると仮定します。市街地のように建物や車が四方八方にある環境では、電波がどの方向からも同じくらいの強さで到来する、という近似はかなり妥当です。この仮定のもとでは、到来角 $\theta$ は $[-\pi, \pi)$ 上で一様分布します。
$$ p_\theta(\theta) = \frac{1}{2\pi}, \qquad \theta \in [-\pi, \pi) $$
加えて、アンテナは水平面内で無指向性（どの方向にも同じ利得 $G$）と仮定します。これらの「あらゆる方向から等しく到来する」という単純な仮定こそが、後で出てくる特徴的なU字型スペクトルの源です。次節では、この一様な到来角分布を周波数領域に「変換」して、ドップラースペクトルの形を導きます。
Classicalドップラースペクトルの導出
ゴール：到来角の分布を周波数の分布に変換する
これから示したいのは、到来角 $\theta$ が一様分布するとき、ドップラー周波数 $f = f_d\cos\theta$ がどんな確率密度（=パワースペクトル密度の形）を持つか、です。結論を先に言うと、$f = \pm f_d$ の両端で発散するU字型になります。これがClarke/Jakesの「Classicalドップラースペクトル」です。
導出の方針は、確率変数の変数変換です。$\theta$ の密度がわかっているとき、$f = f_d\cos\theta$ の密度 $S(f)$ を、ヤコビアン（変数変換の倍率）を使って求めます。
変数変換による導出
受信パワーの方向分布（角度あたりのパワー）を $p(\theta)$ とし、アンテナ利得を $G(\theta)$ とします。周波数 $f$ と $f+df$ の間に入るパワーは、対応する角度区間 $\theta$ と $\theta + d\theta$（および対称な $-\theta$ 側）から来るパワーに等しい、という保存則を立てます。
$$ S(f)\,|df| = G(\theta)\,p(\theta)\,|d\theta| + G(-\theta)\,p(-\theta)\,|d\theta| $$
右辺が2項あるのは、$\cos\theta = \cos(-\theta)$ なので、同じ周波数 $f$ に $+\theta$ と $-\theta$ の2つの角度が寄与するためです。両辺を $|df|$ で割ると、
$$ S(f) = \big[G(\theta)p(\theta) + G(-\theta)p(-\theta)\big]\left|\frac{d\theta}{df}\right| $$
となります。ここで肝心の $\left|d\theta/df\right|$ を計算します。$f = f_d\cos\theta$ を $\theta$ で微分すると、
$$ \frac{df}{d\theta} = -f_d\sin\theta $$
です。$\sin\theta$ を $f$ で表すために、$\cos\theta = f/f_d$ と三角関数の関係 $\sin^2\theta + \cos^2\theta = 1$ を使います。$\sin\theta = \sqrt{1 – \cos^2\theta} = \sqrt{1 – (f/f_d)^2}$ を代入すると、
$$ \left|\frac{df}{d\theta}\right| = f_d\sqrt{1 – \left(\frac{f}{f_d}\right)^2} $$
逆数を取ると、求めたいヤコビアンが得られます。
$$ \left|\frac{d\theta}{df}\right| = \frac{1}{f_d\sqrt{1 – (f/f_d)^2}} $$
ここで、$\theta$ が一様分布（$p(\theta) = 1/2\pi$）で、アンテナが無指向性（$G(\theta) = G$ 一定）という仮定を入れます。すると $G(\theta)p(\theta) + G(-\theta)p(-\theta) = 2 \cdot G/(2\pi) = G/\pi$ となり、これは $f$ によらない定数です。これらをまとめると、
$$ S(f) = \frac{G/\pi}{f_d\sqrt{1 – (f/f_d)^2}} = \frac{G}{\pi f_d}\cdot\frac{1}{\sqrt{1 – (f/f_d)^2}}, \qquad |f| < f_d $$
を得ます。比例係数をまとめて、本質的な形だけ書けば次のとおりです。
$$ \boxed{\;S(f) \propto \frac{1}{\sqrt{1 – (f/f_d)^2}}, \qquad |f| \le f_d\;} $$
なぜU字型になるのか
導出した式 $S(f) \propto 1/\sqrt{1-(f/f_d)^2}$ を眺めましょう。$f \to \pm f_d$ では分母の平方根がゼロに近づくので、$S(f)$ は無限大に発散します。一方、$f = 0$ では $S(0) \propto 1$ で有限です。つまりスペクトルは中央でへこみ、両端で跳ね上がるU字型になります。
この発散には明快な物理的意味があります。$f = +f_d$ に対応する到来角は $\theta = 0$、$f = -f_d$ は $\theta = \pi$ です。ここで $\cos\theta$ は極値を取るため、$\theta$ が少し動いても $\cos\theta = f/f_d$ はほとんど変わりません。言い換えると、$\theta = 0$ や $\pi$ の近傍にある「広い角度範囲」の到来波が、$f = \pm f_d$ という「狭い周波数範囲」に集中して詰め込まれます。だから両端でパワー密度が高くなるのです。逆に $\theta = \pi/2$（真横）の近傍は、$\cos\theta$ の変化が急なので、角度のパワーが広い周波数範囲にばらけ、中央の密度は低くなります。
左図は写像 $f/f_d=\cos\theta$ で、$\theta=0,\pi$ の両端で曲線が平坦（傾きゼロ）になっており、広い角度がわずかな周波数幅に押し込まれることが見て取れます。右図は一様な到来角サンプル40万点を $\cos$ で写像した実測ヒストグラムで、両端 $f/f_d=\pm1$ で密度が跳ね上がり、理論曲線 $1/(\pi\sqrt{1-x^2})$ にぴたりと重なります。U字スペクトルが「幾何」から生まれることが定量的に確認できます。
このU字スペクトルは、Jakesモデルの「指紋」とも言える特徴です。実測のドップラースペクトルがこの形に近いことは、市街地伝搬の多くの場面で確認されています。次は、この周波数領域の特徴を時間領域に翻訳します。スペクトルの逆フーリエ変換が、フェージングの「時間相関」を与えるのです。
自己相関関数 — ベッセル関数 J0 の登場
スペクトルと自己相関はフーリエ対
ここで思い出したいのが、ウィーナー・ヒンチンの定理です。定常確率過程のパワースペクトル密度 $S(f)$ と自己相関関数 $R(\tau)$ は、フーリエ変換のペアになっています。
$$ R(\tau) = \int_{-\infty}^{\infty} S(f)\, e^{\,j 2\pi f \tau}\, df $$
フェージング利得を時間 $t$ の複素確率過程 $g(t)$ とみなすと、その自己相関 $R(\tau) = \mathbb{E}[g(t)g^*(t-\tau)]$ は、ドップラースペクトル $S(f)$ を逆フーリエ変換すれば得られます。「スペクトルの形が変動の速さ（相関の減衰）を決める」という直感が、この式に凝縮されています。
J0 を導く積分
Classicalスペクトル $S(f) = \dfrac{1}{\pi f_d\sqrt{1-(f/f_d)^2}}$（全体のパワーが1になるよう規格化したもの）を上の式に代入します。
$$ R(\tau) = \int_{-f_d}^{f_d} \frac{1}{\pi f_d\sqrt{1-(f/f_d)^2}}\, e^{\,j 2\pi f \tau}\, df $$
積分範囲が $[-f_d, f_d]$ なのは、$|f|>f_d$ では到来波が存在せず $S(f)=0$ だからです。ここで $f = f_d\cos\theta$ と置換します。すると $df = -f_d\sin\theta\, d\theta$ で、$\sqrt{1-(f/f_d)^2} = \sin\theta$（$\theta\in[0,\pi]$ で正）なので、分母とヤコビアンがきれいに打ち消し合います。積分範囲は $f: -f_d\to f_d$ が $\theta: \pi\to 0$ に対応します。
$$ R(\tau) = \int_{\pi}^{0} \frac{1}{\pi f_d \sin\theta}\, e^{\,j 2\pi f_d\tau\cos\theta}\,(-f_d\sin\theta)\,d\theta = \frac{1}{\pi}\int_{0}^{\pi} e^{\,j 2\pi f_d\tau\cos\theta}\, d\theta $$
最後の積分は、ベッセル関数の積分表示そのものです。第1種0次ベッセル関数 $J_0$ には次の表現があります。
$$ J_0(z) = \frac{1}{\pi}\int_{0}^{\pi} e^{\,j z\cos\theta}\, d\theta $$
$z = 2\pi f_d\tau$ と読み替えれば、ただちに次の結論を得ます。
$$ \boxed{\;R(\tau) = J_0(2\pi f_d \tau)\;} $$
この結果が意味すること
$R(\tau) = J_0(2\pi f_d\tau)$ は、フェージング利得の「記憶の長さ」を教えてくれます。$\tau=0$ では $J_0(0)=1$（自分自身とは完全相関）、$\tau$ が増えるにつれて振動しながら減衰します。$J_0$ の最初のゼロ点は引数が約2.405なので、$2\pi f_d\tau = 2.405$、すなわち $\tau \approx 0.38/f_d$ で相関がいったんゼロになります。つまり、おおよそ $1/f_d$ 程度の時間スケールで利得は無相関になります。
ここから「コヒーレンス時間」という重要な概念が出ます。利得がほぼ一定とみなせる時間は $T_c \sim 1/f_d$ のオーダーです。$f_d$ が大きい（=速く動く）ほどコヒーレンス時間は短く、変動が速い。これは冒頭の「速く動くと変動が激しい」という体感と完全に一致します。インターリーバ長やダイバーシティの設計は、この $T_c$ を基準に決めます。
理論はこれで出そろいました。一様到来角 → U字スペクトル → $J_0$ 相関、という連鎖です。問題は「では、この統計を満たす時系列をどうやって計算機で作るか」です。スペクトルを直接フィルタリングで作る方法もありますが、Jakesが提案したのは、もっと直接的で美しい正弦波の和でした。
Jakesの正弦波和モデル — 時系列の作り方
アイデア：少数の正弦波で散乱波を近似する
理論上は無限個の到来波の重ね合わせですが、計算機で無限和は扱えません。Jakesのアイデアは、到来角を $[0, 2\pi)$ 上に等間隔（あるいはほぼ等間隔）に配置した有限個 $N$ の正弦波で、散乱波の集合を近似することです。各正弦波は、対応する到来角のドップラーシフト $f_d\cos\theta_n$ で振動します。これらを足し合わせれば、U字スペクトルと $J_0$ 相関を近似的に再現できる、というわけです。
中心極限定理が効くのがポイントです。独立な正弦波を多数足すと、和は近似的にガウス分布になります。同相成分 $g_I$ と直交成分 $g_Q$ がそれぞれ平均0のガウス過程になれば、複素利得 $g = g_I + j g_Q$ の包絡線 $|g|$ は自動的にレイリー分布、位相は一様分布になります（これは後で確認します）。
同相・直交成分の構成
複素フェージング利得を次の形で生成します。本記事では、古典的Jakesモデルの相関再現性の弱点を改良したZheng–Xiao（鄭・肖）モデルの形を採用します。これは到来角に微小なオフセットを入れることで、$g_I$ と $g_Q$ を独立にし、相関を $J_0$ に精密に合わせたものです。
$$ g_I(t) = \sqrt{\frac{2}{N}}\sum_{n=1}^{N}\cos\!\big(2\pi f_d t\cos\alpha_n + \phi_n\big) $$
$$ g_Q(t) = \sqrt{\frac{2}{N}}\sum_{n=1}^{N}\cos\!\big(2\pi f_d t\sin\alpha_n + \psi_n\big) $$
各記号の意味は次のとおりです。
- $N$：正弦波（散乱波）の本数。$N$ を大きくするほどガウス近似が良くなります。実用上 $N = 8 \sim 16$ で十分です。
- $\alpha_n$：第 $n$ 波の到来角。$\displaystyle \alpha_n = \frac{2\pi n – \pi}{4N}$ と等間隔に配置します。これにより $\cos\alpha_n$ が $[-f_d, f_d]$ のドップラーシフトを満遍なくカバーします。
- $\phi_n, \psi_n$：各正弦波の初期位相。$[-\pi, \pi)$ 上の一様乱数で、試行ごとに独立に選びます。これがフェージングの「ランダムさ」の源です。
係数 $\sqrt{2/N}$ は、$g_I, g_Q$ の分散がそれぞれ1になるよう規格化するためのものです。実際、$\cos$ の二乗平均は $1/2$ なので、$N$ 本の和の分散は $(2/N)\cdot N \cdot (1/2) = 1$ になります。
包絡線がレイリー、位相が一様になる理由
$g_I, g_Q$ がそれぞれ独立に $\mathcal{N}(0, \sigma^2)$（ここでは $\sigma^2=1$）に従うとき、複素利得 $g = g_I + jg_Q$ の包絡線 $r = |g| = \sqrt{g_I^2 + g_Q^2}$ の分布を求めます。2次元ガウスの同時密度を極座標 $(r, \varphi)$ に変換します。$g_I = r\cos\varphi$、$g_Q = r\sin\varphi$ とおき、面積素片の関係 $dg_I\, dg_Q = r\, dr\, d\varphi$ を使うと、
$$ p(g_I, g_Q)\,dg_I\,dg_Q = \frac{1}{2\pi\sigma^2}\exp\!\left(-\frac{g_I^2+g_Q^2}{2\sigma^2}\right)dg_I\,dg_Q = \frac{r}{2\pi\sigma^2}\exp\!\left(-\frac{r^2}{2\sigma^2}\right)dr\,d\varphi $$
となります。ここで $g_I^2+g_Q^2 = r^2$ を代入しました。位相 $\varphi$ について $0$ から $2\pi$ まで積分すると、$r$ の周辺密度が得られます。
$$ p_r(r) = \frac{r}{\sigma^2}\exp\!\left(-\frac{r^2}{2\sigma^2}\right), \qquad r \ge 0 $$
これがまさにレイリー分布です。一方、$r$ について積分すれば位相の周辺密度 $p_\varphi(\varphi) = 1/(2\pi)$ となり、一様分布です。さらに $r$ と $\varphi$ の同時密度が積の形に分かれている（$p(r,\varphi) = p_r(r)p_\varphi(\varphi)$）ので、包絡線と位相は独立です。
つまり、$g_I, g_Q$ をガウス過程として作りさえすれば、レイリーフェージングの統計が自動的に従ってくるのです。Jakesモデルが正弦波和でガウス過程を作りにいくのは、この性質を狙ってのことです。理論的な道具立てが揃ったので、いよいよPythonで時系列を生成し、3つの統計（包絡線・位相・自己相関）を確かめましょう。
Pythonでの実装
ステップ1：複素フェージング利得の時系列生成
まず、Zheng–Xiao型の正弦波和でフェージング利得を生成する関数を作ります。最大ドップラー周波数 $f_d$、サンプリング周波数 $f_s$、サンプル数を引数に取り、複素利得の時系列を返します。
import numpy as np
import matplotlib.pyplot as plt
# 日本語フォント設定（豆腐文字を防ぐ）
import matplotlib.font_manager as fm
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in fm.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
def jakes_fading(fd, fs, num_samples, N=16, seed=None):
"""Zheng-Xiao型の正弦波和でレイリーフェージング利得を生成する
fd: 最大ドップラー周波数[Hz], fs: サンプリング周波数[Hz]
N: 正弦波の本数, 戻り値: 複素利得の時系列(分散1)"""
rng = np.random.default_rng(seed)
t = np.arange(num_samples) / fs
n = np.arange(1, N + 1)
# 到来角を等間隔に配置（Zheng-Xiaoの平均オフセット）
alpha = (2 * np.pi * n - np.pi) / (4 * N)
# 各正弦波の初期位相は一様乱数（ランダムさの源）
phi = rng.uniform(-np.pi, np.pi, N)
psi = rng.uniform(-np.pi, np.pi, N)
# 同相成分と直交成分（外積で時間×正弦波の行列を作って和を取る）
gI = np.sqrt(2.0 / N) * np.cos(2 * np.pi * fd * np.outer(t, np.cos(alpha)) + phi).sum(axis=1)
gQ = np.sqrt(2.0 / N) * np.cos(2 * np.pi * fd * np.outer(t, np.sin(alpha)) + psi).sum(axis=1)
return gI + 1j * gQ
# 時速72km, 2GHz相当のドップラー周波数を仮定
fd = 100.0 # 最大ドップラー周波数[Hz]
fs = 2000.0 # サンプリング周波数[Hz]
g = jakes_fading(fd, fs, num_samples=200000, N=16, seed=1)
print(f"同相成分の分散: {np.var(g.real):.4f}")
print(f"直交成分の分散: {np.var(g.imag):.4f}")
print(f"瞬時電力の平均: {np.mean(np.abs(g) ** 2):.4f}")
実行すると、同相成分・直交成分の分散がともに約1.00（実測で $\mathrm{var}(g_I)=1.0001$、$\mathrm{var}(g_Q)=1.0002$）、瞬時電力 $|g|^2$ の平均が約2.00（実測2.0003）になります。これは $g_I, g_Q$ がそれぞれ分散1のガウス過程として正しく規格化されている証拠です。瞬時電力が2になるのは $\mathbb{E}[|g|^2]=\mathbb{E}[g_I^2]+\mathbb{E}[g_Q^2]=1+1=2$ だからで、平均電力1に揃えたい場合は $g/\sqrt{2}$ とすればよいことがわかります。
生成した同相成分 $g_I(t)$ と直交成分 $g_Q(t)$ の時間波形を見てみましょう。
わずか16本の正弦波の和でありながら、2つの成分はともに平均0の周りで滑らかに揺らぐ、いかにもガウス過程らしい波形になっています。$g_I$ と $g_Q$ は互いに独立に揺らぎ、特定の相関を持たないことも見て取れます。この2つを実部・虚部とする複素利得が、次に見るレイリーフェージングを生み出します。
ステップ2：包絡線の時間波形を眺める
生成した利得の包絡線 $|g(t)|$ をdB表示で時間プロットし、フェージングの「深い落ち込み」を観察します。
t = np.arange(len(g)) / fs
env_db = 20 * np.log10(np.abs(g) + 1e-12)
plt.figure(figsize=(10, 4))
plt.plot(t[:2000], env_db[:2000], color="C0", lw=0.8)
plt.axhline(0, color="gray", ls="--", lw=1, label="平均レベル(0 dB付近)")
plt.xlabel("時間 [s]")
plt.ylabel("包絡線レベル [dB]")
plt.title(f"レイリーフェージングの包絡線波形 ($f_d$={fd:.0f} Hz)")
plt.legend(loc="lower right")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
この波形からは、フェージングの典型的な振る舞いが読み取れます。包絡線は平均レベル（0 dB付近）の周りで激しく上下し、ときおり $-20\,\mathrm{dB}$（赤い点線）を超える深い谷（ディープフェード）が現れます。この谷こそが通信が一瞬途切れる原因です。谷と谷の間隔がおおよそ $1/f_d \approx 10\,\mathrm{ms}$ 程度であることも見て取れ、これがコヒーレンス時間のスケールと一致します。
ステップ3：包絡線がレイリー分布、位相が一様分布になることの確認
統計の核心です。包絡線のヒストグラムをレイリー分布の理論曲線と、位相のヒストグラムを一様分布と重ねて確認します。
from scipy.stats import rayleigh
env = np.abs(g)
phase = np.angle(g)
sigma = np.sqrt(np.mean(env ** 2) / 2) # レイリー分布のスケール母数
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
# 左：包絡線のヒストグラム vs レイリー分布
axes[0].hist(env, bins=80, density=True, alpha=0.6, color="C0", label="シミュレーション")
r = np.linspace(0, env.max(), 400)
axes[0].plot(r, rayleigh.pdf(r, scale=sigma), "r-", lw=2, label="レイリー分布(理論)")
axes[0].set_xlabel("包絡線 $r=|g|$")
axes[0].set_ylabel("確率密度")
axes[0].set_title("包絡線の分布")
axes[0].legend()
axes[0].grid(alpha=0.3)
# 右：位相のヒストグラム vs 一様分布
axes[1].hist(phase, bins=80, density=True, alpha=0.6, color="C2", label="シミュレーション")
axes[1].axhline(1 / (2 * np.pi), color="r", lw=2, label="一様分布(理論) $1/2\\pi$")
axes[1].set_xlabel("位相 $\\varphi$ [rad]")
axes[1].set_ylabel("確率密度")
axes[1].set_title("位相の分布")
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
2つのヒストグラムは、それぞれ理論曲線にぴたりと重なります。左の包絡線は、原点で0から立ち上がり、$r=\sigma$ 付近にピークを持ち、その後なだらかに減衰するレイリー分布の形を正確に再現しています。右の位相は $[-\pi, \pi)$ 上でほぼ平坦で、$1/(2\pi)\approx 0.159$ の水平線に沿っています。前節で導いた「ガウス過程の包絡線はレイリー、位相は一様」が、わずか16本の正弦波の和から見事に成り立っていることが確認できました。
ステップ4：自己相関がベッセル関数 J0 に一致することの実証
最後に、この記事のクライマックスです。生成した利得の自己相関を数値的に計算し、理論値 $J_0(2\pi f_d\tau)$ と重ねます。
from scipy.special import jv
# 自己相関 R(tau) = E[g(t) g*(t-tau)] をラグごとに推定
max_lag = 150
g = jakes_fading(fd, fs, num_samples=400000, N=16, seed=1)
N_samp = len(g)
ac = np.array([np.mean(g[:N_samp - l] * np.conj(g[l:])) for l in range(max_lag)])
ac = ac.real / ac[0].real # tau=0で1になるよう規格化
taus = np.arange(max_lag) / fs
theory = jv(0, 2 * np.pi * fd * taus) # 理論値 J0(2 pi fd tau)
plt.figure(figsize=(10, 4.5))
plt.plot(fd * taus, ac, "o", ms=4, color="C0", label="シミュレーション")
plt.plot(fd * taus, theory, "r-", lw=2, label="$J_0(2\\pi f_d \\tau)$ (理論)")
plt.axhline(0, color="gray", lw=0.8)
plt.xlabel("規格化ラグ $f_d\\,\\tau$")
plt.ylabel("自己相関 $R(\\tau)$")
plt.title("フェージング利得の自己相関とベッセル関数 $J_0$")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print(f"最大誤差: {np.max(np.abs(ac - theory)):.5f}")
プロットを見ると、シミュレーションの点（青丸）が理論曲線 $J_0(2\pi f_d\tau)$（赤線）の上にぴったり乗ります。最大誤差は実測で $5.9\times10^{-5}$ と極めて小さく、Jakesの正弦波和が狙いどおりにClassicalドップラースペクトルに対応する時間相関を再現していることが定量的に確認できます。$f_d\tau \approx 0.38$（緑の点線）で相関が最初にゼロを横切る点も、$J_0$ の第1ゼロ点（引数2.405）から予測される $2.405/2\pi \approx 0.383$ と一致しています。これがコヒーレンス時間 $T_c\sim 1/f_d$ の正体です。
ステップ5：ドップラースペクトルがU字型になることの確認
時間波形のスペクトルを推定して、導出したU字型スペクトルが実際に現れることを見ます。
# 周期グラム（パワースペクトル密度）でドップラースペクトルを推定
g = jakes_fading(fd, fs, num_samples=400000, N=64, seed=3)
G = np.fft.fftshift(np.fft.fft(g))
psd = np.abs(G) ** 2
freqs = np.fft.fftshift(np.fft.fftfreq(len(g), d=1 / fs))
# 観測しやすいよう、近傍を移動平均で平滑化
def smooth(x, w=201):
k = np.ones(w) / w
return np.convolve(x, k, mode="same")
psd_s = smooth(psd)
mask = np.abs(freqs) < 1.3 * fd
plt.figure(figsize=(9, 4.5))
plt.plot(freqs[mask], psd_s[mask] / psd_s[mask].max(), color="C0", label="推定スペクトル(平滑化)")
# 理論のU字型曲線（|f|<fd）
ff = np.linspace(-0.999 * fd, 0.999 * fd, 800)
classic = 1.0 / (np.pi * fd * np.sqrt(1 - (ff / fd) ** 2))
plt.plot(ff, classic / classic[len(ff) // 2] * (psd_s[mask] / psd_s[mask].max())[np.argmin(np.abs(freqs[mask]))],
"r--", lw=2, label="Classicalスペクトル(理論)")
plt.axvline(fd, color="gray", ls=":", lw=1)
plt.axvline(-fd, color="gray", ls=":", lw=1, label="$\\pm f_d$")
plt.xlabel("周波数 [Hz]")
plt.ylabel("正規化パワー")
plt.title("ドップラースペクトル（U字型 / Classical）")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
推定したスペクトルは、中央でへこみ $f=\pm f_d$ の近くで跳ね上がる、特徴的なU字型を示します。理論のClassicalスペクトル（赤破線）と形状がよく一致し、$\pm f_d$ の外側（灰色の点線より外）ではパワーがほとんどゼロになっています。これは導出した $S(f)\propto 1/\sqrt{1-(f/f_d)^2}$（$|f|\le f_d$）そのものです。スペクトルが $\pm f_d$ で帯域制限されている事実は、フェージングの変動が $f_d$ で決まる速さを超えないことを意味しており、時間波形の滑らかさとも整合します。
ステップ6：移動速度（ドップラー周波数）による変動の速さの違い
最後に、$f_d$ を変えると変動の速さがどう変わるかを見て、移動速度との関係を体感します。
fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
for ax, fd_i, vname in zip(axes, [10, 50, 200], ["低速(歩行)", "中速(自動車)", "高速(高速鉄道)"]):
gi = jakes_fading(fd_i, fs, num_samples=4000, N=16, seed=7)
t_i = np.arange(len(gi)) / fs
ax.plot(t_i, 20 * np.log10(np.abs(gi) + 1e-12), color="C0", lw=0.9)
ax.set_ylabel("レベル[dB]")
ax.set_title(f"{vname}: $f_d$={fd_i} Hz（コヒーレンス時間 ≈ {1000/fd_i:.0f} ms）")
ax.grid(alpha=0.3)
ax.set_ylim(-35, 10)
axes[-1].set_xlabel("時間 [s]")
plt.tight_layout()
plt.show()
3つの波形を比べると、$f_d$ が大きいほど包絡線の変動が速く、谷の出現頻度も高いことが一目でわかります。$f_d=10\,\mathrm{Hz}$（歩行程度）ではゆったりと数百ミリ秒スケールで変動するのに対し、$f_d=200\,\mathrm{Hz}$（高速鉄道相当）では数ミリ秒の間に何度も深い谷が現れます。これは導出したコヒーレンス時間 $T_c\sim 1/f_d$ がそのまま現れた結果であり、「速く動くほど通信路の変動が速い」という移動体通信の根本原理を可視化しています。受信機の追従制御や符号化の設計が速度に依存する理由が、この図から直観できます。
ステップ7：レベルクロス率と平均フェード長の検証
フェージング対策の設計では「包絡線がどれくらいの頻度であるレベルを下回り、一度下回ると平均どれだけの時間そこに留まるか」が重要です。前者をレベルクロス率（LCR）、後者を平均フェード長（AFD）と呼びます。正規化しきい値 $\rho=R/R_{\mathrm{rms}}$ に対し、理論値は $\mathrm{LCR}=\sqrt{2\pi}\,f_d\rho\,e^{-\rho^2}$、$\mathrm{AFD}=(e^{\rho^2}-1)/(\sqrt{2\pi}\,f_d\rho)$ で与えられます。生成した利得から実測し、理論と比べます。
g = jakes_fading(fd, fs, num_samples=2000000, N=16, seed=11)
env = np.abs(g)
rrms = np.sqrt(np.mean(env ** 2))
total_t = len(g) / fs
rho_db = np.linspace(-30, 5, 28)
rho = 10 ** (rho_db / 20)
lcr_sim, afd_sim = [], []
for rl in rho:
above = env > rl * rrms
upcross = np.sum((~above[:-1]) & (above[1:])) # 上向きクロス回数
lcr_sim.append(upcross / total_t)
afd_sim.append((np.mean(~above) * total_t) / upcross)
lcr_th = np.sqrt(2 * np.pi) * fd * rho * np.exp(-rho ** 2)
afd_th = (np.exp(rho ** 2) - 1) / (np.sqrt(2 * np.pi) * fd * rho)
実測のLCRはしきい値 $\rho=-5$ dB付近でピークを持ち、理論曲線とよく一致します（$\rho=-10$ dBで実測65.6回/秒、理論68.1回/秒）。平均フェード長も $\rho=-10$ dBで実測1.24 ms・理論1.24 msとよく合います。深い谷ほどクロス頻度は下がるが一度落ちると長く留まる、という設計上重要な傾向が定量的に確認できます（極端に低いしきい値では稀なイベントのサンプル不足で理論からずれます）。
ステップ8：正弦波本数 N とガウス近似の収束
最後に、正弦波の本数 $N$ を変えると包絡線分布がどのようにレイリーへ近づくかを見ます。中心極限定理が効くには十分な本数が必要なはずです。
import numpy as np
r_grid = np.linspace(0, 3.5, 300)
for N in [2, 4, 8, 16]:
g = jakes_fading(fd, fs, num_samples=300000, N=N, seed=21)
env = np.abs(g)
hist, edges = np.histogram(env, bins=120, range=(0, 3.5), density=True)
# 各Nについて hist をプロット（中心 0.5*(edges[:-1]+edges[1:])）
$N=2$ では包絡線分布がレイリー曲線から大きく外れ、ガウス近似が破綻していることがわかります。$N$ を増やすにつれて分布は急速にレイリーへ近づき、$N=8$ で既にほぼ一致、$N=16$ では理論曲線と区別がつきません。これが本文で「実用上 $N=8\sim16$ で十分」とした根拠で、少数の正弦波でも中心極限定理が良く効くことを実測が裏付けています。
まとめ
本記事では、移動体フェージングを再現するJakesモデルを、理論の導出から実装・検証まで通して解説しました。
- 一様到来角からU字スペクトルへ: 散乱体が全方位に均等に分布し到来角が一様分布すると仮定すると、変数変換のヤコビアン $1/(f_d\sqrt{1-(f/f_d)^2})$ を通じて、Classical（U字型）ドップラースペクトル $S(f)\propto 1/\sqrt{1-(f/f_d)^2}$ が導かれます。両端の発散は「広い角度範囲が狭い周波数に詰め込まれる」幾何に由来します。
- スペクトルから J0 相関へ: ウィーナー・ヒンチンの定理でスペクトルを逆フーリエ変換し、$f=f_d\cos\theta$ の置換でベッセル関数の積分表示に帰着させると、自己相関 $R(\tau)=J_0(2\pi f_d\tau)$ が得られます。これがコヒーレンス時間 $T_c\sim 1/f_d$ を決めます。
- 正弦波和による時系列生成: 少数（$N=8\sim16$本）の正弦波の和でガウス過程を作ると、包絡線は自動的にレイリー分布、位相は一様分布になります。Zheng–Xiao型の角度オフセットを使えば、自己相関を $J_0$ に最大誤差 $6\times10^{-5}$ の精度で合わせられます。
- 数値的検証: Pythonで生成した利得が、包絡線レイリー・位相一様・自己相関 $J_0$・U字スペクトルという4つの統計をすべて満たすことを確認しました。移動速度を変えると変動の速さが $f_d$ に比例して変わることも可視化しました。
Jakesモデルは、わずかな正弦波の和という最小限の仕掛けで、移動体無線チャネルの本質的な統計を捉える優れたツールです。受信機のBER評価、ダイバーシティやインターリーバの設計、等化器の追従性能の検証など、無線システム設計のあらゆる場面で土台になります。
次のステップとして、以下の記事も参考にしてください。
- 仲上・ライスフェージングチャネルの理論 — 見通し波がある場合への拡張。レイリーは $K$ ファクタが0の極限です
- フェージングとダイバーシティ — 変動の影響と対策 — 本記事で求めた時間相関が、ダイバーシティ設計にどう効くか
- マルチパス伝搬とは — 多重波がもたらす干渉とフェージング — フェージングの物理的な出発点