スマートフォンで通話しながら部屋の中を歩くと、ある場所では声がクリアに聞こえ、別の場所では急にノイズが増えることがあります。これは、基地局から届く電波が壁や家具で何度も反射し、複数の経路を通って受信機に到達した結果、互いに強め合ったり弱め合ったりするためです。この現象を「マルチパスフェージング」と呼びます。とりわけ、基地局が見える窓際のような環境では、反射波に加えて遮蔽物のない直接波（見通し波, Line-of-Sight）が強く届きます。このような「強い直接波 + たくさんの弱い散乱波」という状況を統計的に記述するのが、本記事で扱うライスフェージング（Rician fading）です。
ライスフェージングを理解すると、無線通信の設計や評価において次のような場面で見通しが格段によくなります。
- 衛星通信・GNSS: 衛星からの直接波が支配的だが、地表や建物による反射波が混入する。直接波と散乱波の比（Kファクタ）が受信品質を左右する
- マイクロセル・屋内無線（Wi-Fi, 5G）: 基地局やアクセスポイントが見える環境ではライス分布、見えない環境ではレイリー分布に近づく。チャネルモデルの選択が誤り率予測の精度を決める
- ドローン・航空機の通信リンク: 上空からは地表が見通せるため強い直接波が存在し、ライスフェージングのKファクタが大きくなる
本記事の内容
- 直接波 + 複素ガウス散乱波という物理モデルの設定
- 包絡線が従うライス分布の確率密度関数の導出（省略なし）
- Kファクタ$K$の定義と、$K \to 0$でレイリー分布、$K \to \infty$で無フェージングに漸近する連続性
- ライス分布のモーメントとモーメント法によるKファクタ推定
- Pythonによる分布生成・ヒストグラム照合・包絡線時系列シミュレーション・K推定・BER比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ライスフェージングとは
電波が空間を伝わって受信アンテナに届くとき、その経路は1本とは限りません。送信機から一直線に届く直接波のほかに、ビルの壁・地面・天井・家具などで反射した無数の散乱波が、それぞれ少しずつ異なる時間差（位相差）をもって受信機に重なります。受信機が観測するのは、これらすべての波のベクトル和です。
ここで重要なのは、各波がもつ「位相」がランダムに散らばっているという点です。反射経路の長さは波長（数センチ〜数十センチ）のスケールで見ると簡単に1波長以上変動するため、受信機が少し動いただけで各散乱波の位相は $0$ から $2\pi$ までほぼ一様にランダムに変わります。強め合う位相関係になれば信号は大きくなり、打ち消し合えば信号はほとんど消えてしまう。これがフェージング（信号強度のランダムな変動）の正体です。
散乱波だけが存在し、突出して強い波が1本もない環境では、受信信号の包絡線（振幅）はレイリー分布に従うことが知られています。市街地で基地局が建物の陰に隠れているような「見通しなし（NLOS）」環境がこれにあたります。
一方、基地局が直接見える窓際や、衛星のように遮蔽されない送信源がある場合は、ランダムな散乱波の山の上に、振幅が一定で位相も安定した直接波が乗ります。このとき包絡線が従うのがライス分布です。直感的には「ランダムなゆらぎ（散乱波）+ どっしりした土台（直接波）」という構図で、土台が高いほど信号は安定し、土台が消えるとレイリー分布に戻ります。この「土台の強さ」を定量化するのがKファクタです。
以下では、この物理的な描像を数式モデルに落とし込み、包絡線の確率密度関数を厳密に導出していきます。
受信信号の物理モデル
同相成分・直交成分への分解
受信信号を複素ベースバンド表現（同相 $I$ 成分と直交 $Q$ 成分の複素数）で表します。受信複素振幅を $g$ とすると、直接波と散乱波の和は次のように書けます。
$$ g = \underbrace{A e^{j\theta_0}}_{\text{直接波}} + \underbrace{\sum_{i=1}^{N} a_i e^{j\phi_i}}_{\text{散乱波}} $$
ここで $A$ は直接波の振幅（定数）、$\theta_0$ は直接波の位相、$a_i$ と $\phi_i$ は $i$ 番目の散乱波の振幅と位相です。散乱波の本数 $N$ は十分大きいと仮定します。
散乱波の和を実部と虚部に分けます。
$$ \sum_{i=1}^{N} a_i e^{j\phi_i} = \underbrace{\sum_{i=1}^{N} a_i \cos\phi_i}_{X’} + j\underbrace{\sum_{i=1}^{N} a_i \sin\phi_i}_{Y’} $$
各 $\phi_i$ は $[0, 2\pi)$ で独立にほぼ一様分布し、$a_i$ は互いに独立とします。すると $X’$ と $Y’$ は、独立で同分布な確率変数を多数足し合わせたものになります。ここで中心極限定理が効いてきます。$N$ が大きいとき、$X’$ と $Y’$ はそれぞれ平均ゼロのガウス分布に従い、しかも互いに独立になります。この事実は、散乱波の細かい振幅・位相分布の詳細によらず成り立つ、非常に強い結論です。
複素ガウス散乱波としての定式化
中心極限定理の帰結を整理しましょう。直接波の位相 $\theta_0$ を基準（$\theta_0 = 0$）に取り直すと、受信複素振幅 $g$ の実部・虚部は次のように書けます。
$$ g = X + jY, \qquad X = A + X’, \quad Y = Y’ $$
ここで $X’$, $Y’$ は平均 $0$、分散 $\sigma^2$ の独立な正規分布に従う確率変数です。
$$ X’ \sim \mathcal{N}(0, \sigma^2), \qquad Y’ \sim \mathcal{N}(0, \sigma^2) $$
したがって、同相成分 $X$ は平均 $A$、直交成分 $Y$ は平均 $0$ のガウス分布に従います。
$$ X \sim \mathcal{N}(A, \sigma^2), \qquad Y \sim \mathcal{N}(0, \sigma^2) $$
直接波は同相成分の平均を $0$ から $A$ へ「押し上げる」役割を果たし、直交成分には影響しません（位相基準を直接波に合わせたため）。$\sigma^2$ は散乱波の作るゆらぎの大きさを表し、各方向（$I$, $Q$）の散乱電力を意味します。
私たちが最終的に知りたいのは、受信信号の包絡線（振幅）
$$ R = |g| = \sqrt{X^2 + Y^2} $$
がどんな確率分布に従うか、という点です。包絡線が大きいときは信号が強く（通信品質が良く）、小さいときは深いフェージングに落ち込んで通信が途絶えやすくなります。
物理モデルが整いました。次は、$X$ と $Y$ の同時分布から、包絡線 $R$ の確率密度関数を導出します。
ライス分布の確率密度関数の導出
同時確率密度から極座標へ
$X$ と $Y$ は独立なガウス分布なので、その同時確率密度関数は各々の密度の積になります。
$$ p_{X,Y}(x, y) = \frac{1}{2\pi\sigma^2} \exp\!\left(-\frac{(x-A)^2 + y^2}{2\sigma^2}\right) $$
私たちが欲しいのは包絡線 $R = \sqrt{X^2 + Y^2}$ の分布なので、直交座標 $(x, y)$ から極座標 $(r, \psi)$ へ変数変換します。
$$ x = r\cos\psi, \qquad y = r\sin\psi $$
ここで $r \geq 0$ は包絡線、$\psi \in [0, 2\pi)$ は受信信号の位相です。変数変換のヤコビアンは $|\partial(x,y)/\partial(r,\psi)| = r$ なので、$(r, \psi)$ の同時密度は次のようになります。
$$ p_{R,\Psi}(r, \psi) = r \cdot p_{X,Y}(r\cos\psi, r\sin\psi) $$
これを具体的に書き下します。指数の中身を展開する準備として、$(x-A)^2 + y^2$ に $x = r\cos\psi$, $y = r\sin\psi$ を代入します。
$$ (r\cos\psi – A)^2 + (r\sin\psi)^2 = r^2\cos^2\psi – 2Ar\cos\psi + A^2 + r^2\sin^2\psi $$
$\cos^2\psi + \sin^2\psi = 1$ を使うと、$r^2\cos^2\psi + r^2\sin^2\psi = r^2$ にまとまります。したがって、
$$ (x-A)^2 + y^2 = r^2 – 2Ar\cos\psi + A^2 $$
これを同時密度に代入すると、
$$ p_{R,\Psi}(r, \psi) = \frac{r}{2\pi\sigma^2} \exp\!\left(-\frac{r^2 + A^2 – 2Ar\cos\psi}{2\sigma^2}\right) $$
ここまでで、包絡線 $r$ と位相 $\psi$ の同時分布が得られました。次に、位相 $\psi$ について積分して消去し、$r$ だけの分布（周辺分布）を求めます。
位相についての周辺化とベッセル関数
包絡線 $R$ の周辺密度は、$\psi$ を $0$ から $2\pi$ まで積分すれば得られます。
$$ p_R(r) = \int_0^{2\pi} p_{R,\Psi}(r, \psi)\, d\psi $$
指数の中身を $\psi$ に依存する項とそうでない項に分けます。$\psi$ に依存するのは $2Ar\cos\psi$ の項だけなので、
$$ p_R(r) = \frac{r}{2\pi\sigma^2} \exp\!\left(-\frac{r^2 + A^2}{2\sigma^2}\right) \int_0^{2\pi} \exp\!\left(\frac{Ar\cos\psi}{\sigma^2}\right) d\psi $$
と、$\psi$ に依存しない因子を積分の外へ出せます。残った積分は、第1種0次変形ベッセル関数 $I_0$ の積分表示そのものです。変形ベッセル関数 $I_0$ は次のように定義されます。
$$ I_0(z) = \frac{1}{2\pi}\int_0^{2\pi} e^{z\cos\psi}\, d\psi $$
この定義式と見比べると、$z = Ar/\sigma^2$ と置けば積分は $2\pi I_0(Ar/\sigma^2)$ に等しいことがわかります。
$$ \int_0^{2\pi} \exp\!\left(\frac{Ar\cos\psi}{\sigma^2}\right) d\psi = 2\pi I_0\!\left(\frac{Ar}{\sigma^2}\right) $$
これを代入すると、$2\pi$ が約分されて、ついにライス分布の確率密度関数が得られます。
$$ \boxed{\; p_R(r) = \frac{r}{\sigma^2} \exp\!\left(-\frac{r^2 + A^2}{2\sigma^2}\right) I_0\!\left(\frac{Ar}{\sigma^2}\right), \quad r \geq 0 \;} $$
これがライス分布（Rician distribution）です。$A$ が直接波の振幅、$\sigma^2$ が散乱波の各成分の分散を表します。$I_0$ は引数が大きいほど指数関数的に増大する関数で、直接波と散乱波の相互作用（コヒーレントな強め合い）を担っています。
なお、変形ベッセル関数 $I_0$ は級数で次のように表せます。
$$ I_0(z) = \sum_{k=0}^{\infty} \frac{1}{(k!)^2}\left(\frac{z}{2}\right)^{2k} = 1 + \frac{z^2}{4} + \frac{z^4}{64} + \cdots $$
$z$ が小さい（$A$ が小さい）ときは $I_0(z) \approx 1$ となり、後で見るようにレイリー分布に帰着します。
確率密度関数が導けたので、次はこの分布を特徴づける最も重要なパラメータ、Kファクタを導入します。
Kファクタの定義と極限挙動
Kファクタの定義
ライス分布には $A$（直接波振幅）と $\sigma^2$（散乱分散）という2つのパラメータがあります。しかし、フェージングの「質」を決めるのは、これらの絶対値ではなく直接波と散乱波の電力比です。これを表すのがKファクタ（Rician K-factor）です。
$$ K = \frac{\text{直接波の電力}}{\text{散乱波の電力}} = \frac{A^2}{2\sigma^2} $$
分子の $A^2$ は直接波の電力（振幅の2乗）、分母の $2\sigma^2$ は散乱波の電力です。散乱波の電力が $2\sigma^2$ になるのは、$I$ 成分と $Q$ 成分がそれぞれ分散 $\sigma^2$ を持ち、合計の散乱電力が $\sigma^2 + \sigma^2 = 2\sigma^2$ になるためです。
Kファクタはしばしばデシベルで表されます。
$$ K_{\text{dB}} = 10\log_{10} K $$
直感的には、$K$ が大きいほど「どっしりした直接波の土台」が高く、フェージングが浅くなります。$K$ が小さいほど散乱波のランダムなゆらぎが支配的になり、深いフェージングが頻繁に起こります。典型的な値として、見通しの良い衛星リンクでは $K \approx 10$〜$20$ dB、市街地のマイクロセルでは $K \approx 0$〜$10$ dB、見通しなしの環境では $K \to -\infty$ dB（$K=0$）です。
総受信電力とパラメータの関係
総受信電力 $\Omega$ は、直接波電力と散乱波電力の和です。
$$ \Omega = \mathbb{E}[R^2] = A^2 + 2\sigma^2 $$
Kファクタの定義 $K = A^2/(2\sigma^2)$ と合わせると、$A^2$ と $\sigma^2$ を $K$ と $\Omega$ で表せます。$A^2 = K \cdot 2\sigma^2$ を $\Omega = A^2 + 2\sigma^2$ に代入すると $\Omega = 2\sigma^2(K+1)$ となり、
$$ 2\sigma^2 = \frac{\Omega}{K+1}, \qquad A^2 = \frac{K\,\Omega}{K+1} $$
が得られます。この関係は、シミュレーションで「総電力を $\Omega=1$ に固定したまま $K$ だけ変えたい」というときに便利です。$K$ を決めれば $A$ と $\sigma$ が一意に定まります。
$K \to 0$ の極限：レイリー分布
直接波が消える（$A \to 0$、すなわち $K \to 0$）と、ライス分布はどうなるでしょうか。確率密度関数で $A \to 0$ とすると、$I_0(Ar/\sigma^2) \to I_0(0) = 1$ となり、ベッセル関数の項が消えます。
$$ p_R(r) \xrightarrow{A \to 0} \frac{r}{\sigma^2}\exp\!\left(-\frac{r^2}{2\sigma^2}\right) $$
これはまさにレイリー分布の確率密度関数です。物理的にも自然で、直接波がなくなれば「ランダムな散乱波だけ」の状況、すなわちレイリーフェージングに帰着します。つまりレイリー分布はライス分布の特別な場合（$K=0$）なのです。
$K \to \infty$ の極限：無フェージング（AWGN）
逆に、散乱波が直接波に比べて無視できるほど小さくなる（$\sigma \to 0$、すなわち $K \to \infty$）と何が起きるでしょうか。このとき包絡線 $R$ はほぼ $A$ に集中します。
これを定量的に見るために、$K$ が大きいときのライス分布の挙動を調べます。$z = Ar/\sigma^2$ が大きいとき、変形ベッセル関数は次の漸近形を持ちます。
$$ I_0(z) \approx \frac{e^z}{\sqrt{2\pi z}} \quad (z \gg 1) $$
これをライス分布に代入し、$r$ が $A$ の近傍にあるとして整理すると、指数部は
$$ -\frac{r^2 + A^2}{2\sigma^2} + \frac{Ar}{\sigma^2} = -\frac{(r-A)^2}{2\sigma^2} $$
と平方完成できます。したがって $K \to \infty$ では、
$$ p_R(r) \approx \frac{1}{\sqrt{2\pi}\,\sigma}\exp\!\left(-\frac{(r-A)^2}{2\sigma^2}\right) $$
となり、包絡線は平均 $A$、標準偏差 $\sigma$ のガウス分布に近づきます。$\sigma \to 0$ の極限では分布の幅がゼロに縮み、$R = A$ の一定値（デルタ関数）になります。これは振幅が揺らがない無フェージングの状況、すなわち加法性白色ガウス雑音（AWGN）だけが残るチャネルに対応します。
このように、Kファクタという1つのパラメータを $0$ から $\infty$ まで動かすだけで、レイリーフェージング（最悪のケース）から無フェージング（理想的なケース）まで連続的につながるのが、ライス分布の美しさです。
極限挙動を押さえたので、次はKファクタを実データから推定するための土台として、ライス分布のモーメントを計算します。
ライス分布のモーメントとKファクタ推定
モーメントの計算
ライス分布の $n$ 次モーメント $\mathbb{E}[R^n]$ は、一般には合流型超幾何関数 $_1F_1$（クンマー関数）で表されます。
$$ \mathbb{E}[R^n] = (2\sigma^2)^{n/2}\,\Gamma\!\left(1+\frac{n}{2}\right)\, {}_1F_1\!\left(-\frac{n}{2};\, 1;\, -\frac{A^2}{2\sigma^2}\right) $$
ここで $\Gamma$ はガンマ関数です。一般形は複雑ですが、推定に使う2次モーメント（電力）は前節で見たように非常にシンプルです。
$$ \mathbb{E}[R^2] = A^2 + 2\sigma^2 = \Omega $$
これは超幾何関数を経由せずとも、$R^2 = X^2 + Y^2$ の期待値が $\mathbb{E}[X^2] + \mathbb{E}[Y^2] = (A^2 + \sigma^2) + \sigma^2 = A^2 + 2\sigma^2$ となることから直接確かめられます。
4次モーメントも比較的簡単な形になります。$X, Y$ がガウスであることを使ってモーメントを展開すると、
$$ \mathbb{E}[R^4] = \mathbb{E}[(X^2+Y^2)^2] = A^4 + 8A^2\sigma^2 + 8\sigma^4 $$
が得られます。これらの2次・4次モーメントを使うと、Kファクタを観測データから推定するモーメント法が構成できます。
モーメント法によるK推定
包絡線の電力 $R^2$ の分散を考えます。$\mathrm{Var}[R^2] = \mathbb{E}[R^4] – (\mathbb{E}[R^2])^2$ を計算すると、
$$ \mathrm{Var}[R^2] = (A^4 + 8A^2\sigma^2 + 8\sigma^4) – (A^2 + 2\sigma^2)^2 $$
右辺を展開します。$(A^2 + 2\sigma^2)^2 = A^4 + 4A^2\sigma^2 + 4\sigma^4$ なので、
$$ \mathrm{Var}[R^2] = (8A^2\sigma^2 + 8\sigma^4) – (4A^2\sigma^2 + 4\sigma^4) = 4A^2\sigma^2 + 4\sigma^4 = 4\sigma^2(A^2 + \sigma^2) $$
ここで、よく使われる推定量として、電力の2乗平均と分散の比に着目します。具体的には次の量を定義します。
$$ \gamma = \frac{(\mathbb{E}[R^2])^2}{\mathrm{Var}[R^2]} = \frac{(A^2 + 2\sigma^2)^2}{4\sigma^2(A^2 + \sigma^2)} $$
この $\gamma$ を $K = A^2/(2\sigma^2)$ で書き直してみましょう。分子・分母を $\sigma^4$ で割り、$A^2/\sigma^2 = 2K$ を代入すると、$A^2 + 2\sigma^2 = \sigma^2(2K+2)$、$A^2 + \sigma^2 = \sigma^2(2K+1)$ なので、
$$ \gamma = \frac{\sigma^4(2K+2)^2}{4\sigma^4(2K+1)} = \frac{(2K+2)^2}{4(2K+1)} = \frac{(K+1)^2}{2K+1} $$
このように、$\gamma$ は $K$ だけの関数になります。逆に解けば、観測した $\gamma$ から $K$ が求まります。$\gamma(2K+1) = (K+1)^2$ を $K$ について整理すると、$K^2 + (2 – 2\gamma)K + (1 – \gamma) = 0$ という2次方程式になり、物理的に意味のある正の根を取ると、
$$ \boxed{\; \hat{K} = \frac{(\gamma – 1) + \sqrt{\gamma^2 – \gamma}}{1} = (\gamma – 1) + \sqrt{\gamma(\gamma-1)} \;} $$
が得られます（$\gamma \geq 1$ で実数解を持ち、$\gamma=1$ のとき $K=0$ すなわちレイリーに対応します）。
実用上は、観測された包絡線サンプル $\{r_i\}$ から $\mathbb{E}[R^2]$ と $\mathrm{Var}[R^2]$ を標本平均・標本分散で置き換え、$\gamma$ を計算して上式に代入すれば、Kファクタが推定できます。この方法は「モーメント法（method of moments）」と呼ばれ、計算が軽くロバストなため、実機での測定によく使われます。
モーメント法のほかに、2次・4次モーメントを直接使うGreenstein法として知られる推定式
$$ \hat{K} = \frac{\sqrt{1 – 1/\gamma}}{1 – \sqrt{1 – 1/\gamma}} $$
もよく用いられます。これは上で導いた式と数学的に等価です（$\gamma = (K+1)^2/(2K+1)$ を逆に解いた別表現）。後ほどPythonで両者が一致することを確認します。
推定式が整いました。理論はここまでとして、次はPythonでライス分布を実際に生成し、これらの結果が正しいことを数値的に検証していきましょう。
Pythonによる実装と可視化
ライス分布の確率密度関数とサンプル生成
まず、導出したライス分布の確率密度関数を関数として実装し、異なるKファクタでの形状を比較します。総電力を $\Omega = 1$ に固定し、$K$ を変えたときの分布の変化を見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import i0 # 第1種0次変形ベッセル関数 I0
def rician_pdf(r, K, Omega=1.0):
"""ライス分布の確率密度関数 (総電力 Omega を固定)"""
# K, Omega から A^2 と sigma^2 を計算
two_sigma2 = Omega / (K + 1.0) # 2*sigma^2
sigma2 = two_sigma2 / 2.0
A2 = K * Omega / (K + 1.0) # A^2
A = np.sqrt(A2)
# p(r) = (r/sigma^2) exp(-(r^2 + A^2)/(2 sigma^2)) I0(A r / sigma^2)
return (r / sigma2) * np.exp(-(r**2 + A2) / (2 * sigma2)) * i0(A * r / sigma2)
r = np.linspace(0, 3, 500)
K_list = [0.0, 1.0, 4.0, 16.0] # K=0 はレイリー
plt.figure(figsize=(8, 5))
for K in K_list:
label = "K=0 (Rayleigh)" if K == 0 else f"K={K:.0f} ({10*np.log10(K):.1f} dB)"
plt.plot(r, rician_pdf(r, K), label=label, lw=2)
plt.xlabel("envelope r")
plt.ylabel("probability density p(r)")
plt.title("Rician PDF for various K-factors (Omega=1)")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("rician_pdf.png", dpi=150)
plt.show()
上のグラフから、Kファクタが分布の形状を劇的に変えることが読み取れます。$K=0$（レイリー）では分布が原点近くに偏り、低い包絡線（深いフェージング）が頻繁に起こります。$K$ を大きくしていくと、ピークが包絡線 $r=1$ 付近（直接波の振幅に対応）へ移動し、分布の幅が狭くなって左右対称なガウス分布に近づいていきます。これは前節で導いた「$K\to\infty$ でガウス分布に漸近する」という結論と完全に一致しています。直接波が強いほど、受信信号は安定するのです。
ヒストグラムと理論曲線の照合
次に、物理モデル（同相・直交成分のガウス和）から直接サンプルを生成し、その包絡線のヒストグラムが理論曲線と一致することを確認します。これは導出が正しいことの数値的な裏付けになります。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import i0
def generate_rician_envelope(K, N, Omega=1.0, rng=None):
"""物理モデル: 直接波(定数) + 複素ガウス散乱波 の包絡線を生成"""
if rng is None:
rng = np.random.default_rng(0)
sigma2 = Omega / (2 * (K + 1.0)) # 各成分(I,Q)の分散
A = np.sqrt(K * Omega / (K + 1.0)) # 直接波振幅
X = A + np.sqrt(sigma2) * rng.standard_normal(N) # 同相: 平均A
Y = np.sqrt(sigma2) * rng.standard_normal(N) # 直交: 平均0
return np.sqrt(X**2 + Y**2)
def rician_pdf(r, K, Omega=1.0):
sigma2 = Omega / (2 * (K + 1.0))
A2 = K * Omega / (K + 1.0)
A = np.sqrt(A2)
return (r / sigma2) * np.exp(-(r**2 + A2) / (2 * sigma2)) * i0(A * r / sigma2)
rng = np.random.default_rng(42)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
for ax, K in zip(axes, [1.0, 16.0]):
samples = generate_rician_envelope(K, 200000, rng=rng)
ax.hist(samples, bins=120, density=True, alpha=0.5,
color="steelblue", label="histogram (simulated)")
r = np.linspace(0, samples.max(), 400)
ax.plot(r, rician_pdf(r, K), "r-", lw=2, label="theoretical PDF")
ax.set_title(f"K = {K:.0f} ({10*np.log10(K):.1f} dB)")
ax.set_xlabel("envelope r"); ax.set_ylabel("density")
ax.legend(); ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("rician_hist.png", dpi=150)
plt.show()
2つのヒストグラムは、いずれも赤い理論曲線にぴったり重なっています。これは、独立な2つのガウス成分（平均 $A$ の同相成分と平均 $0$ の直交成分）の和の包絡線がライス分布に従うという、本記事の導出が正しいことを裏付けています。左の $K=1$ では分布がやや非対称で裾が広い一方、右の $K=16$ では $r=1$ 付近に鋭くピークが立ち、ガウス分布に近づいているのが見て取れます。
包絡線時系列のシミュレーション
実際の無線環境では、受信機や反射物が動くことで包絡線が時間とともに変動します。ここでは、相関を持たせた散乱成分を時系列として生成し、レイリー（$K=0$）とライス（$K=10$）でフェージングの深さがどう違うかを比較します。
import numpy as np
import matplotlib.pyplot as plt
def fading_time_series(K, N, Omega=1.0, smooth=20, rng=None):
"""簡易な相関付きフェージング包絡線の時系列"""
if rng is None:
rng = np.random.default_rng(1)
sigma = np.sqrt(Omega / (2 * (K + 1.0)))
A = np.sqrt(K * Omega / (K + 1.0))
# 白色ガウスを移動平均で平滑化し、時間相関(ドップラー)を模擬
win = np.ones(smooth) / smooth
Xp = np.convolve(rng.standard_normal(N + smooth), win, mode="valid")[:N]
Yp = np.convolve(rng.standard_normal(N + smooth), win, mode="valid")[:N]
# 平滑化で分散が縮むので元の sigma に正規化
Xp *= sigma / Xp.std(); Yp *= sigma / Yp.std()
return np.sqrt((A + Xp)**2 + Yp**2)
rng = np.random.default_rng(7)
N = 2000
env_rayleigh = fading_time_series(0.0, N, rng=rng)
env_rician = fading_time_series(10.0, N, rng=rng)
plt.figure(figsize=(10, 5))
plt.plot(20*np.log10(env_rayleigh), label="K=0 (Rayleigh)", alpha=0.8)
plt.plot(20*np.log10(env_rician), label="K=10 (Rician)", alpha=0.8)
plt.axhline(0, color="k", lw=0.8, ls="--")
plt.xlabel("sample index (time)")
plt.ylabel("envelope level [dB]")
plt.title("Fading envelope time series: Rayleigh vs Rician")
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("fading_timeseries.png", dpi=150)
plt.show()
このグラフから、レイリー（$K=0$）の包絡線はしばしば $-20$ dB を下回る深いフェードに落ち込んでいるのに対し、ライス（$K=10$）の包絡線は $0$ dB 付近に安定して張り付き、深いフェードがほとんど起きていないことがわかります。直接波の土台がある分、瞬間的に信号が消える「ディープフェード」のリスクが大幅に減るのです。これが、見通しのある環境で通信が安定する理由を直感的に示しています。
モーメント法によるKファクタ推定
次に、生成した包絡線サンプルから真のKファクタを復元できるか、モーメント法を実装して検証します。本記事で導いた推定式と、等価なGreenstein法の両方を実装し、結果が一致することも確認します。
import numpy as np
def generate_rician_envelope(K, N, Omega=1.0, rng=None):
if rng is None:
rng = np.random.default_rng(0)
sigma2 = Omega / (2 * (K + 1.0))
A = np.sqrt(K * Omega / (K + 1.0))
X = A + np.sqrt(sigma2) * rng.standard_normal(N)
Y = np.sqrt(sigma2) * rng.standard_normal(N)
return np.sqrt(X**2 + Y**2)
def estimate_K_moment(envelope):
"""包絡線サンプルからモーメント法でKを推定"""
P = envelope**2 # 瞬時電力 R^2
gamma = np.mean(P)**2 / np.var(P) # gamma = (E[R^2])^2 / Var[R^2]
# 本記事の式: K = (gamma-1) + sqrt(gamma(gamma-1))
K_ours = (gamma - 1) + np.sqrt(max(gamma * (gamma - 1), 0.0))
# Greenstein法: K = sqrt(1-1/gamma) / (1 - sqrt(1-1/gamma))
s = np.sqrt(max(1 - 1/gamma, 0.0))
K_green = s / (1 - s) if (1 - s) > 0 else np.inf
return K_ours, K_green
rng = np.random.default_rng(123)
print(f"{'true K':>8} | {'K_dB':>6} | {'est (ours)':>10} | {'est (Green)':>11}")
for K_true in [0.0, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0]:
env = generate_rician_envelope(K_true, 500000, rng=rng)
Ko, Kg = estimate_K_moment(env)
kdb = (10*np.log10(K_true)) if K_true > 0 else -np.inf
print(f"{K_true:>8.1f} | {kdb:>6.1f} | {Ko:>10.3f} | {Kg:>11.3f}")
この出力表を見ると、推定されたKファクタ（est (ours) と est (Green)）が、真の値 true K にきわめて近い値を返していることが確認できます。さらに、本記事で導いた式とGreenstein法の推定値が小数点以下まで一致しており、両者が数学的に等価であるという主張も裏付けられました。$K=0$（レイリー）では推定値もほぼ $0$ になり、モーメント法が極限ケースでも破綻しないことがわかります。サンプル数を増やすほど推定精度は上がります。
BERへの影響：レイリーとライスの比較
最後に、フェージングが通信のビット誤り率（BER）に与える影響を見ます。BPSK変調を仮定し、平均SNRを横軸に、レイリーとライス（複数のK）でBERがどう違うかをモンテカルロシミュレーションで比較します。瞬時SNRは包絡線の2乗に比例するため、フェージングが浅い（Kが大きい）ほど誤り率が下がるはずです。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc
def ber_bpsk_fading(K, avg_snr_dB, N=2_000_000, rng=None):
"""フェージング下のBPSK BER をモンテカルロで計算"""
if rng is None:
rng = np.random.default_rng(0)
avg_snr = 10**(avg_snr_dB / 10)
Omega = 1.0
sigma2 = Omega / (2 * (K + 1.0))
A = np.sqrt(K * Omega / (K + 1.0))
# 包絡線の2乗 |h|^2 (平均 Omega=1)
X = A + np.sqrt(sigma2) * rng.standard_normal(N)
Y = np.sqrt(sigma2) * rng.standard_normal(N)
h2 = X**2 + Y**2
# 瞬時SNR = 平均SNR * |h|^2 / E[|h|^2]
inst_snr = avg_snr * h2 / np.mean(h2)
# BPSK の瞬時BER を平均
return np.mean(0.5 * erfc(np.sqrt(inst_snr)))
snr_dB = np.arange(0, 31, 2)
rng = np.random.default_rng(2024)
plt.figure(figsize=(8, 6))
for K, style in [(0.0, "o-"), (3.0, "s-"), (10.0, "^-")]:
ber = [ber_bpsk_fading(K, s, rng=rng) for s in snr_dB]
lbl = "K=0 (Rayleigh)" if K == 0 else f"K={K:.0f}"
plt.semilogy(snr_dB, ber, style, label=lbl)
# AWGN (フェージングなし) の理論値
ber_awgn = 0.5 * erfc(np.sqrt(10**(snr_dB/10)))
plt.semilogy(snr_dB, ber_awgn, "k--", label="AWGN (no fading)")
plt.xlabel("average SNR [dB]"); plt.ylabel("BER")
plt.title("BPSK BER under Rician/Rayleigh fading")
plt.ylim(1e-6, 1); plt.legend(); plt.grid(alpha=0.3, which="both")
plt.tight_layout()
plt.savefig("ber_fading.png", dpi=150)
plt.show()
このBER曲線から、フェージングが通信品質に与える深刻な影響が読み取れます。レイリー（$K=0$）のBERはSNRに対してゆるやかにしか下がらず、低いBERを得るには非常に高いSNRが必要です。一方、$K$ が大きくなるにつれて曲線は急峻になり、AWGN（フェージングなし、黒破線）の理論曲線に近づいていきます。これは、直接波が強いほどディープフェードによる誤りが減り、無フェージングの理想状態に漸近するためです。前節で導いた「$K\to\infty$ でAWGNに帰着する」という極限挙動が、BER性能の面からも確認できました。実システムでは、ダイバーシティ技術を使ってレイリー環境を実効的にライス的（フェージングの浅い状態）に近づけることが、性能改善の鍵になります。
まとめ
本記事では、見通し成分のあるマルチパス環境を記述するライスフェージングチャネルについて、物理モデルから確率分布、推定、性能評価までを一貫して解説しました。
- 物理モデル: 受信信号は「振幅一定の直接波」と「中心極限定理により複素ガウスとなる散乱波」の和であり、同相成分は平均 $A$、直交成分は平均 $0$ のガウス分布に従う
- ライス分布: 包絡線 $R=\sqrt{X^2+Y^2}$ は $p_R(r)=\frac{r}{\sigma^2}\exp(-\frac{r^2+A^2}{2\sigma^2})I_0(\frac{Ar}{\sigma^2})$ に従い、これを極座標変換と変形ベッセル関数 $I_0$ の積分表示から省略なく導いた
- Kファクタ: $K=A^2/(2\sigma^2)$ は直接波と散乱波の電力比。$K\to 0$ でレイリー分布、$K\to\infty$ でガウス（無フェージング・AWGN）へ連続的に漸近する
- K推定: 電力の2乗平均と分散の比 $\gamma=(\mathbb{E}[R^2])^2/\mathrm{Var}[R^2]$ から $\hat{K}=(\gamma-1)+\sqrt{\gamma(\gamma-1)}$ で推定でき、Greenstein法と等価
- 性能への影響: Kが大きいほどディープフェードが減り、BER曲線がAWGNに近づく
ライスフェージングは、レイリーフェージングとAWGNという2つの極端なチャネルモデルを1つのパラメータでつなぐ、無線伝搬の中核的なモデルです。Kファクタを実測・推定することで、リンクの品質を定量的に評価し、変調方式や誤り訂正符号、ダイバーシティの設計に活かせます。
次のステップとして、以下の記事も参考にしてください。