マルチパス伝搬の理論 — 遅延拡がり・ドップラー拡がりとチャネルモデル

スマートフォンで通話しながら歩いていると、突然音声が途切れたり歪んだりした経験はないでしょうか。屋内でWi-Fiの速度が場所によって大きく変わるのも、同じ現象が原因です。送信アンテナから放射された電磁波は、建物や地面、車両などに反射・回折・散乱しながら、複数の異なる経路を通って受信アンテナに到達します。これがマルチパス伝搬（multipath propagation）です。

各経路を通る信号は、それぞれ異なる遅延・振幅・位相を持って受信点に到着するため、互いに強め合ったり打ち消し合ったりします。この干渉パターンは受信点の位置や周波数によって激しく変動し、通信品質を大きく左右します。マルチパス伝搬を定量的に理解することは、以下のような技術の設計基盤となります。

• OFDM（直交周波数分割多重） — 周波数選択性フェージングに対抗するためにサブキャリア幅を設計する際、コヒーレンス帯域幅の知識が不可欠です
• 等化器（イコライザ）設計 — マルチパスによる符号間干渉を除去するために、チャネルのインパルス応答モデルが必要です
• MIMOアンテナ — マルチパスの空間的広がりを積極的に利用して通信容量を高めます
• 衛星測位（GPS/GNSS） — マルチパスは測距誤差の主要因であり、その特性を把握して軽減する必要があります

本記事の内容

• マルチパス伝搬の物理メカニズムと直感的理解
• 遅延プロファイルと RMS 遅延拡がりの定義・導出
• コヒーレンス帯域幅 $B_c \approx 1/(5\sigma_\tau)$ の導出
• 周波数選択性フェージングとフラットフェージングの判定基準
• ドップラー効果と最大ドップラー周波数の導出
• ドップラースペクトル（ジェイクスモデル）の理論
• コヒーレンス時間 $T_c \approx 1/(4f_d)$
• タップ遅延線チャネルモデルの定式化
• 広帯域チャネル vs 狭帯域チャネルの分類
• Python によるマルチパスチャネルのシミュレーションと可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 自由空間損失とフリスの伝達公式を導出する — 自由空間での電波伝搬の基本
• フェージングとダイバーシティ技術 — フェージングの分類とその対策

マルチパス伝搬とは

自由空間との違い

自由空間では送信アンテナから受信アンテナへ直接波（Line of Sight, LOS）だけが到達します。フリスの伝達公式が示すように、受信電力は距離の2乗に反比例して減衰するだけで、信号の波形そのものは変化しません。

しかし現実の環境では、電磁波は建物の壁、地面、山、車、人体など、さまざまな物体と相互作用します。この相互作用は大きく3つのメカニズムに分類できます。

1. 反射（Reflection） — 波長に比べて十分大きい平面に電磁波が当たると、入射角と等しい反射角で反射波が生じます。建物の壁面や地面での反射が典型例です
2. 回折（Diffraction） — 建物の角や山の稜線など、鋭いエッジ付近で電磁波が回り込む現象です。直接波が遮られた影の領域にも信号が到達する理由はこれです
3. 散乱（Scattering） — 波長と同程度かそれ以下のサイズの物体（街灯、木の葉、電線など）に電磁波が当たると、あらゆる方向に散乱波が生じます

これらのメカニズムにより、同じ送信信号が複数の経路（パス）を通って受信点に到達します。各パスは経路長が異なるため、到着時刻（遅延）が異なります。また、反射回数や散乱の態様によって振幅と位相も異なります。

受信信号の数学的表現

送信信号を $s(t)$ とすると、$N$ 本のパスを通って到達する受信信号は次のように書けます。

$$ \begin{equation} r(t) = \sum_{i=0}^{N-1} a_i(t) \, s(t – \tau_i(t)) \, e^{j\phi_i(t)} \end{equation} $$

ここで各パス $i$ について、

• $a_i(t)$: 振幅（時間とともにゆっくり変動する）
• $\tau_i(t)$: 遅延時間（経路長 $/$ 光速）
• $\phi_i(t)$: 位相シフト（反射による位相変化＋伝搬による位相回転）

です。時間 $t$ に依存しているのは、送受信機や散乱体が移動すると各パスの幾何学的関係が刻々と変化するためです。

この式の物理的意味を噛み砕くと、受信点には「同じ信号の、時間的にずれた・強さの違う・位相の異なるコピー」が大量に重ね合わされて届く、ということです。これらのコピーが互いに干渉することで、受信電力や信号波形が複雑に変動します。

マルチパス伝搬の影響を体系的に理解するには、時間領域と周波数領域の両面から特性を定量化する必要があります。まず時間領域の特性量である遅延プロファイルと遅延拡がりから見ていきましょう。

遅延プロファイルと遅延拡がり

パワー遅延プロファイル

マルチパスチャネルの時間領域特性を表す最も基本的な量がパワー遅延プロファイル（PDP: Power Delay Profile）です。

イメージとしては、受信点に非常に短いパルス（インパルス）を送った場合に、各パスを通ったパルスがどのタイミングでどれだけの強さで到着するかを示す「到着パターン」です。温泉の山から石を落としたとき、大小さまざまな石がバラバラのタイミングで谷底に到達するような状況を想像してください。

数学的には、チャネルのインパルス応答 $h(\tau)$ の電力の期待値として定義されます。

$$ \begin{equation} P(\tau) = \mathbb{E}\!\left[|h(\tau)|^2\right] \end{equation} $$

ここで $\tau$ は基準パス（通常は最初に到着するパス）からの超過遅延（excess delay）です。$P(\tau)$ は一般に $\tau$ が大きくなるにつれて減衰します。これは、長い経路を通った信号ほど多くの反射・散乱を経験し、エネルギーが失われるためです。

実測では、ある場所で多数のスナップショットを取得し、その平均として PDP を得ます。

遅延拡がりの定義

パワー遅延プロファイルから、マルチパスの「時間的な広がり度合い」を単一の数値で表す指標が遅延拡がり（delay spread）です。

まず、PDP の全電力（0次モーメント）を定義します。

$$ \begin{equation} P_{\text{total}} = \int_0^{\infty} P(\tau) \, d\tau \end{equation} $$

次に、PDP を確率密度関数のように正規化して、遅延の平均遅延（mean excess delay）を求めます。これは PDP の1次モーメントです。

$$ \begin{equation} \bar{\tau} = \frac{\int_0^{\infty} \tau \, P(\tau) \, d\tau}{\int_0^{\infty} P(\tau) \, d\tau} \end{equation} $$

さらに、遅延の2次モーメントは次の通りです。

$$ \begin{equation} \overline{\tau^2} = \frac{\int_0^{\infty} \tau^2 \, P(\tau) \, d\tau}{\int_0^{\infty} P(\tau) \, d\tau} \end{equation} $$

ここから、RMS 遅延拡がり（RMS delay spread） $\sigma_\tau$ を定義します。これは PDP の標準偏差に相当します。

$$ \begin{equation} \sigma_\tau = \sqrt{\overline{\tau^2} – \bar{\tau}^2} \end{equation} $$

$\sigma_\tau$ は、マルチパスによる信号の時間的な「散らばり具合」を表す最も重要なパラメータです。値が大きいほど、遅延の広がりが大きい（＝遠回りして来るパスの影響が相対的に大きい）ことを意味します。

離散パスモデルの場合

実際のチャネル測定やシミュレーションでは、$N$ 本の離散パスとして扱うことが多いです。各パスの電力を $P_k = |a_k|^2$、遅延を $\tau_k$ とすると、積分が和に置き換わります。

$$ \begin{equation} \bar{\tau} = \frac{\sum_{k=0}^{N-1} P_k \, \tau_k}{\sum_{k=0}^{N-1} P_k} \end{equation} $$

同様に、遅延の2次モーメントも電力重み付き平均で求まります:

$$ \begin{equation} \overline{\tau^2} = \frac{\sum_{k=0}^{N-1} P_k \, \tau_k^2}{\sum_{k=0}^{N-1} P_k} \end{equation} $$

これら2つのモーメントを用いると、連続の場合と同じ公式で RMS 遅延拡がりが計算できます:

$$ \begin{equation} \sigma_\tau = \sqrt{\overline{\tau^2} – \bar{\tau}^2} \end{equation} $$

典型的な値

環境によって $\sigma_\tau$ は桁違いに異なります。

屋内環境では反射体までの距離が短いため遅延拡がりが小さく、山岳地帯のように遠方からの反射波がある環境では非常に大きくなります。光速 $c = 3 \times 10^8$ m/s で考えると、$\sigma_\tau = 1\,\mu$s は約 300 m の経路長差に対応します。

遅延拡がりの大きさは、信号の帯域幅との相対関係によって通信への影響が決まります。次に、この関係を定量化するコヒーレンス帯域幅を導出します。

コヒーレンス帯域幅

周波数相関関数

マルチパスチャネルを周波数領域で見ると、各周波数成分が異なる振幅・位相変化を受けます。直感的には、2つの周波数が十分に近ければチャネルの応答はほぼ同じですが、離れすぎると互いに無相関になります。この「チャネルが同じ応答を示すとみなせる周波数の範囲」がコヒーレンス帯域幅（coherence bandwidth）です。

これを定式化するために、チャネルの周波数相関関数を導入します。チャネルの周波数応答 $H(f)$ はインパルス応答 $h(\tau)$ のフーリエ変換です。

$$ \begin{equation} H(f) = \int_{-\infty}^{\infty} h(\tau) \, e^{-j2\pi f\tau} \, d\tau \end{equation} $$

2つの周波数 $f_1$ と $f_2$ におけるチャネル応答の相関をとります。周波数間隔を $\Delta f = f_2 – f_1$ として、周波数相関関数を次のように定義します。

$$ \begin{equation} R_H(\Delta f) = \mathbb{E}\!\left[H^*(f) \, H(f + \Delta f)\right] \end{equation} $$

ウィーナー–ヒンチンの定理による導出

ここで、WSSUS（Wide-Sense Stationary Uncorrelated Scattering）チャネルモデルを仮定します。これは、異なる遅延のパスが統計的に無相関であるという仮定で、実用上広く用いられる標準的な近似です。

WSSUS仮定のもとでは、次の重要な関係が成り立ちます。

$$ \begin{equation} \mathbb{E}\!\left[h^*(\tau_1) \, h(\tau_2)\right] = P(\tau_1) \, \delta(\tau_1 – \tau_2) \end{equation} $$

異なる遅延のパス間の相関がゼロ（デルタ関数で表現）であるという意味です。

この関係を用いて、$R_H(\Delta f)$ を導出します。$H(f)$ の定義を代入すると、

$$ R_H(\Delta f) = \mathbb{E}\!\left[\int h^*(\tau_1) e^{j2\pi f\tau_1} d\tau_1 \int h(\tau_2) e^{-j2\pi(f+\Delta f)\tau_2} d\tau_2\right] $$

期待値を積分の中に入れ、WSSUS仮定を適用すると、

$$ R_H(\Delta f) = \int\!\int P(\tau_1)\,\delta(\tau_1 – \tau_2)\,e^{j2\pi f\tau_1}\,e^{-j2\pi(f+\Delta f)\tau_2}\,d\tau_1\,d\tau_2 $$

$\delta(\tau_1 – \tau_2)$ により $\tau_2 = \tau_1$ として積分が潰れます。

$$ R_H(\Delta f) = \int P(\tau)\,e^{j2\pi f\tau}\,e^{-j2\pi(f+\Delta f)\tau}\,d\tau $$

指数関数をまとめると $e^{j2\pi f\tau} \cdot e^{-j2\pi(f+\Delta f)\tau} = e^{-j2\pi \Delta f \cdot \tau}$ なので、

$$ \begin{equation} R_H(\Delta f) = \int_0^{\infty} P(\tau) \, e^{-j2\pi \Delta f \, \tau} \, d\tau \end{equation} $$

これは極めて重要な結果です。周波数相関関数はパワー遅延プロファイルのフーリエ変換に等しいのです。これはウィーナー–ヒンチンの定理のチャネル版といえます。

コヒーレンス帯域幅の導出

コヒーレンス帯域幅 $B_c$ は、周波数相関関数 $R_H(\Delta f)$ がある閾値以上を保つ周波数間隔として定義されます。閾値の取り方によって異なる定義がありますが、最も広く用いられるのは以下の2つです。

相関 0.9 以上の場合:

$$ \begin{equation} B_c \approx \frac{1}{50\,\sigma_\tau} \end{equation} $$

相関 0.5 以上の場合（よく使われる）:

$$ \begin{equation} B_c \approx \frac{1}{5\,\sigma_\tau} \end{equation} $$

後者の導出を概略的に示します。$P(\tau)$ がガウス型や指数減衰型など、RMS幅 $\sigma_\tau$ を持つ分布であるとき、そのフーリエ変換のスペクトル幅はフーリエ変換の不確定性関係により $\sigma_\tau$ の逆数に比例します。

具体的に、$P(\tau) = \frac{1}{\sigma_\tau}\exp\!\left(-\frac{\tau}{\sigma_\tau}\right)$ （指数減衰型）を仮定してフーリエ変換を計算すると、

$$ R_H(\Delta f) = \int_0^{\infty} \frac{1}{\sigma_\tau} e^{-\tau/\sigma_\tau} \, e^{-j2\pi \Delta f \, \tau} \, d\tau $$

これは $a = 1/\sigma_\tau + j2\pi\Delta f$ とおいてラプラス変換の基本公式を適用すれば、

$$ R_H(\Delta f) = \frac{1}{\sigma_\tau} \cdot \frac{1}{1/\sigma_\tau + j2\pi\Delta f} = \frac{1}{1 + j2\pi\sigma_\tau \Delta f} $$

この絶対値は次の通りです。

$$ |R_H(\Delta f)| = \frac{1}{\sqrt{1 + (2\pi\sigma_\tau \Delta f)^2}} \end{equation} $$

$|R_H(\Delta f)| = 0.5$ となる $\Delta f$ を求めると、

$$ \frac{1}{\sqrt{1 + (2\pi\sigma_\tau \Delta f)^2}} = 0.5 $$

両辺を2乗して整理すると、

$$ 1 + (2\pi\sigma_\tau \Delta f)^2 = 4 $$

$(2\pi\sigma_\tau \Delta f)^2 = 3$ を $\Delta f$ について解くと、コヒーレンス帯域幅の目安が得られます。

$$ \Delta f = \frac{\sqrt{3}}{2\pi\sigma_\tau} \approx \frac{0.276}{\sigma_\tau} $$

厳密には $1/(5\sigma_\tau) = 0.2/\sigma_\tau$ とは一致しませんが、$P(\tau)$ の形状に依存する近似であり、工学的な目安として $B_c \approx 1/(5\sigma_\tau)$ が広く採用されています。この近似は多くの実測チャネルモデルに対して十分な精度を持つことが経験的に確認されています。

数値例

先ほどの典型的な遅延拡がりの値を使って、コヒーレンス帯域幅を計算してみましょう。

屋内では 4 MHz もの帯域にわたってチャネルが「平坦」ですが、市街地や山岳地帯ではわずか数十 kHz の帯域幅でもチャネル応答が変化します。この違いが、次に説明する周波数選択性フェージングの判定に直結します。

周波数選択性フェージングとフラットフェージング

信号帯域幅との比較

コヒーレンス帯域幅 $B_c$ は、それ単独では「大きい」とも「小さい」とも言えません。重要なのは、送信信号の帯域幅 $W$ との比較です。

この比較は、郵便受けの幅（$B_c$）と届く荷物の幅（$W$）の関係に似ています。荷物が郵便受けに余裕をもって入るなら形は崩れませんが、郵便受けより大きい荷物を無理に押し込めば変形してしまいます。

フラットフェージング（flat fading）: $W \ll B_c$

信号帯域幅がコヒーレンス帯域幅より十分小さい場合、信号の全周波数成分が同じチャネル応答を受けます。信号の波形は変化せず、振幅だけがスケーリングされます。これがフラットフェージングです。

$$ \begin{equation} W \ll B_c \iff T_s \gg \sigma_\tau \end{equation} $$

ここで $T_s \approx 1/W$ はシンボル時間です。右辺は「シンボル時間がマルチパスの遅延拡がりよりはるかに長い」ことを意味し、各パスの信号が同一シンボル内に収まるため、隣接シンボルへの漏れ込み（符号間干渉, ISI）が起きないという直感と一致します。

周波数選択性フェージング（frequency-selective fading）: $W \gg B_c$

信号帯域幅がコヒーレンス帯域幅を超える場合、周波数によってチャネルの応答が異なります。ある周波数成分は強く増幅され、別の周波数成分は深く減衰するなど、信号スペクトルが歪みます。

$$ \begin{equation} W \gg B_c \iff T_s \ll \sigma_\tau \end{equation} $$

シンボル時間よりもマルチパスの広がりが長いため、前のシンボルの遅延波が次のシンボルに被って ISI が生じます。

具体的な判定例

携帯電話（LTE）を例に考えます。LTEのサブキャリア間隔は 15 kHz です。市街地の $B_c \approx 67$ kHz と比較すると、$15 \text{ kHz} \ll 67 \text{ kHz}$ なので各サブキャリアはフラットフェージングを受けます。しかし、全体の帯域幅が 20 MHz であれば $20 \text{ MHz} \gg 67 \text{ kHz}$ なので、帯域全体では周波数選択性フェージングです。これこそが OFDM の設計思想 — 広帯域信号を狭帯域のサブキャリアに分割して、各サブキャリアではフラットフェージングとして扱う — の根拠です。

ここまでの議論は、マルチパスチャネルの時間分散（遅延拡がり → コヒーレンス帯域幅 → 周波数選択性）に関するものでした。次に、もう一つの重要な側面である時間変動（ドップラー効果）を見ていきます。

ドップラー効果と最大ドップラー周波数

移動に伴う周波数シフト

送信機または受信機（あるいは周囲の散乱体）が移動すると、受信される電磁波の周波数が変化します。これがドップラー効果です。救急車のサイレンが近づくと高く、遠ざかると低く聞こえるのと同じ原理です。

受信機が速度 $v$ で移動しており、ある電磁波が受信機の移動方向と角度 $\theta$ をなして到来する場合、そのドップラーシフトは次のように表されます。

$$ \begin{equation} f_d(\theta) = \frac{v}{\lambda}\cos\theta = f_{\max}\cos\theta \end{equation} $$

ここで $\lambda$ は搬送波の波長、$f_{\max}$ は最大ドップラー周波数（maximum Doppler frequency）です。

$$ \begin{equation} f_{\max} = \frac{v}{\lambda} = \frac{v \, f_c}{c} \end{equation} $$

$f_c$ は搬送波周波数、$c$ は光速です。

この式の導出を示します。受信機が方向 $\hat{\bm{v}}$ に速度 $v$ で移動し、電磁波が方向 $\hat{\bm{k}}$ から到来する場合を考えます。時間 $\Delta t$ の間に受信機は $v\Delta t$ だけ移動し、電磁波の到来方向への射影成分は $v\Delta t \cos\theta$ です。この距離変化による位相変化は、

$$ \Delta\phi = \frac{2\pi}{\lambda} \cdot v\Delta t\cos\theta $$

周波数は位相の時間変化率を $2\pi$ で割ったものなので、ドップラーシフトは次のようになります。

$$ f_d = \frac{1}{2\pi}\frac{d\phi}{dt} = \frac{v\cos\theta}{\lambda} $$

$\cos\theta = 1$（受信機に正面から向かってくる波）のとき最大値 $f_{\max} = v/\lambda$ をとり、$\cos\theta = -1$（受信機から遠ざかる方向）のとき $-f_{\max}$ となります。

数値例

歩行者の Wi-Fi 利用では $f_{\max}$ はわずか 12 Hz ですが、新幹線で 5G を使うと約 1 kHz にも達します。高速移動環境ではドップラー効果への対策が不可欠であることがわかります。

マルチパス環境では、各パスの到来角 $\theta_i$ がそれぞれ異なるため、各パスが受けるドップラーシフトも異なります。この「ドップラーシフトのばらつき」を統計的に記述するのが、次に紹介するドップラースペクトルです。

ドップラースペクトルとジェイクスモデル

ドップラーパワースペクトル

マルチパス環境で受信される信号の周波数変動を統計的に記述するために、ドップラーパワースペクトル $S(f)$ を導入します。これは受信信号のパワースペクトル密度を搬送波周波数からの偏差 $f$ の関数として表したものです。

$S(f)$ の形状は、電磁波の到来角分布に依存します。到来角分布が一様であるか、特定の方向に偏っているかによって、ドップラースペクトルの形状が変わります。

ジェイクスモデル（クラークモデル）

都市部の移動通信環境でよく用いられるのがジェイクスモデル（Jakes model）、正式にはクラークの散乱モデルです。このモデルでは以下を仮定します。

1. 多数の散乱波が水平面内のあらゆる方向から一様に到来する（2D等方散乱）
2. 各散乱波の振幅は等しく、位相は独立一様分布

この仮定のもとで、到来角 $\theta$ は $[0, 2\pi)$ 上の一様分布に従います。$\theta$ の確率密度関数は $p(\theta) = 1/(2\pi)$ です。

ドップラーシフトは $f_d = f_{\max}\cos\theta$ なので、$f_d$ と $\theta$ の関係から変数変換を行ってドップラースペクトルを導きます。

$f_d = f_{\max}\cos\theta$ を $\theta$ について解くと、

$$ \theta = \arccos\!\left(\frac{f_d}{f_{\max}}\right) $$

$\cos\theta$ は $[0, \pi]$ と $[\pi, 2\pi]$ でそれぞれ単調なので、同じ $f_d$ に対応する $\theta$ が2つあります。変数変換の公式 $p(f_d) = p(\theta) \cdot |d\theta/df_d|$ の寄与を2つ分足し合わせます。

$$ \frac{d\theta}{df_d} = \frac{-1}{f_{\max}\sin\theta} $$

$\sin\theta = \sqrt{1 – \cos^2\theta} = \sqrt{1 – (f_d/f_{\max})^2}$ を代入すると、

$$ \left|\frac{d\theta}{df_d}\right| = \frac{1}{f_{\max}\sqrt{1 – (f_d/f_{\max})^2}} $$

$p(\theta) = 1/(2\pi)$ と2つの $\theta$ からの寄与を合わせると、

$$ S(f_d) = 2 \times \frac{1}{2\pi} \times \frac{1}{f_{\max}\sqrt{1 – (f_d/f_{\max})^2}} $$

整理して次のジェイクスドップラースペクトルを得ます。

$$ \begin{equation} S(f_d) = \frac{1}{\pi f_{\max}\sqrt{1 – \left(\frac{f_d}{f_{\max}}\right)^2}}, \quad |f_d| < f_{\max} \end{equation} $$

$|f_d| \geq f_{\max}$ では $S(f_d) = 0$ です。

ジェイクススペクトルの特徴

このスペクトルにはいくつかの重要な特徴があります。

1. バスタブ型 — $f_d = \pm f_{\max}$ 近傍でスペクトル密度が発散します（特異点）。これは、$\cos\theta \approx \pm 1$ の方向（受信機の移動方向とその逆方向）からの波が $\theta$ の微小変化に対して $f_d$ がほとんど変化しないため、エネルギーが集中するからです
2. 有限帯域 — $|f_d| > f_{\max}$ ではスペクトルがゼロです。ドップラーシフトは $f_{\max}$ を超えません
3. 対称性 — $S(f_d)$ は $f_d = 0$ に関して対称です。前方と後方から到来する波の統計的寄与が等しいためです

実際の環境では散乱が完全に一様ではないため、特異点は有限値にとどまりますが、ジェイクスモデルは2D散乱環境の良い近似として広く使われています。

ドップラースペクトルの幅はチャネルの時間変動の速さを決定します。次に、この幅と関連するコヒーレンス時間を導出します。

コヒーレンス時間

時間相関関数

コヒーレンス帯域幅がチャネルの周波数領域での「安定幅」を表すのと双対的に、コヒーレンス時間（coherence time） $T_c$ はチャネルの時間領域での「安定幅」を表します。コヒーレンス時間内であれば、チャネルの応答はほぼ一定とみなせます。

チャネルの時間相関関数は次のように定義されます。

$$ \begin{equation} R_h(\Delta t) = \mathbb{E}\!\left[h^*(t) \, h(t + \Delta t)\right] \end{equation} $$

遅延拡がりと周波数相関関数がフーリエ変換対であったのと同様に、ドップラーパワースペクトルと時間相関関数もフーリエ変換対です。

$$ \begin{equation} R_h(\Delta t) = \int_{-\infty}^{\infty} S(f_d) \, e^{j2\pi f_d \Delta t} \, df_d \end{equation} $$

ジェイクスモデルのスペクトルを代入して積分を評価すると、結果は第一種ベッセル関数 $J_0$ になります。

$$ \begin{equation} R_h(\Delta t) = J_0(2\pi f_{\max} \Delta t) \end{equation} $$

$J_0$ は $x = 0$ で値1をとり、振動しながら減衰していきます。$J_0$ が最初にゼロを横切るのは $x \approx 2.4048$ のときなので、$R_h(\Delta t) = 0$ となる最初の時間差は、

$$ 2\pi f_{\max} \Delta t \approx 2.4048 \implies \Delta t \approx \frac{0.383}{f_{\max}} $$

コヒーレンス時間の近似

実用上のコヒーレンス時間には複数の定義が使われます。

相関が0.5を下回る時間差の定義:

$$ \begin{equation} T_c \approx \frac{9}{16\pi f_{\max}} \approx \frac{0.179}{f_{\max}} \end{equation} $$

これは $J_0(2\pi f_{\max} T_c) = 0.5$ から得られますが、やや保守的な値です。

幾何平均的な定義（広く使われる）:

ドップラー拡がりと遅延拡がりの双対関係から、コヒーレンス帯域幅の定義と同様の考え方で、

$$ \begin{equation} T_c \approx \frac{1}{4 f_{\max}} \end{equation} $$

とする定義が工学的に広く用いられます。

高速フェージングと低速フェージング

コヒーレンス帯域幅がフラットフェージングと周波数選択性フェージングの判定に使われたのと同様に、コヒーレンス時間は高速フェージングと低速フェージングの判定に使われます。

低速フェージング（slow fading）: $T_s \ll T_c$

シンボル時間がコヒーレンス時間より十分短い場合、1シンボルの間にチャネルはほぼ変化しません。

高速フェージング（fast fading）: $T_s \gg T_c$

シンボル時間がコヒーレンス時間を超える場合、1シンボルの間にチャネルが変動し、信号が歪みます。

数値例

先ほどのドップラー周波数の値を用いてコヒーレンス時間を計算します。

LTEのスロット長は 0.5 ms、サブフレーム長は 1 ms なので、自動車走行中の $T_c \approx 1$ ms はちょうどサブフレーム長と同程度です。新幹線では $T_c$ がスロット長を下回るため、チャネル推定の頻度を上げるなどの対策が必要になります。

ここまでで、マルチパスチャネルの特性を記述する4つの基本量 — 遅延拡がり $\sigma_\tau$、コヒーレンス帯域幅 $B_c$、最大ドップラー周波数 $f_{\max}$、コヒーレンス時間 $T_c$ — が出揃いました。次に、これらを統一的に扱うチャネルモデルを定式化します。

タップ遅延線チャネルモデル

モデルの動機

通信システムのシミュレーションや等化器の設計では、マルチパスチャネルを離散時間のフィルタとして表現できると便利です。連続時間のインパルス応答をサンプリング定理に基づいて離散化したものがタップ遅延線（TDL: Tapped Delay Line）チャネルモデルです。

イメージとしては、水道管を流れる水が複数に分岐し、それぞれ異なる長さ・太さの管を通って再び合流する様子を想像してください。各分岐が「タップ」に対応し、管の長さが遅延、太さが振幅に対応します。

定式化

信号のサンプリング間隔を $T_s = 1/W$（$W$ は信号帯域幅）とし、チャネルのインパルス応答を $T_s$ 間隔でサンプリングします。タップ数 $L$ は遅延拡がりと信号帯域幅から決まり、次のように見積もられます。

$$ \begin{equation} L \approx \left\lfloor \frac{\tau_{\max}}{T_s} \right\rfloor + 1 = \left\lfloor \tau_{\max} \cdot W \right\rfloor + 1 \end{equation} $$

ここで $\tau_{\max}$ は最大超過遅延（PDP が有意な値を持つ最大の遅延）です。

離散時間チャネルモデルは、時刻 $n$ における入力サンプル $x[n]$ に対して、出力を次のように表します。

$$ \begin{equation} y[n] = \sum_{l=0}^{L-1} h_l[n] \, x[n – l] + w[n] \end{equation} $$

ここで、

• $h_l[n]$: 第 $l$ タップの複素係数（時刻 $n$ で変動する）
• $x[n-l]$: $l$ サンプル遅延した入力信号
• $w[n]$: 加法性白色ガウス雑音（AWGN）

各タップ係数 $h_l[n]$ は、遅延 $\tau_l = l \cdot T_s$ 近傍に到着するパス群の合成を表します。多数の散乱パスの重ね合わせであるため、中心極限定理により $h_l[n]$ は複素ガウス分布に従います。

各タップの統計的性質

直接波が存在しない場合（NLOS: Non-Line of Sight）、各タップの実部と虚部はそれぞれ独立な平均ゼロのガウス分布に従います。したがって振幅 $|h_l[n]|$ はレイリー分布に従います。

$$ p(r) = \frac{r}{\sigma_l^2} \exp\!\left(-\frac{r^2}{2\sigma_l^2}\right), \quad r \geq 0 $$

ここで $\sigma_l^2 = P(\tau_l)/2$ は第 $l$ タップの PDP から決まる分散です。

直接波が存在する場合（LOS: Line of Sight）には、振幅はライス分布に従います。

$$ p(r) = \frac{r}{\sigma_l^2} \exp\!\left(-\frac{r^2 + A^2}{2\sigma_l^2}\right) I_0\!\left(\frac{rA}{\sigma_l^2}\right) $$

$A$ は直接波の振幅、$I_0$ は変形ベッセル関数です。ライスファクター $K = A^2/(2\sigma_l^2)$ が直接波と散乱波の電力比を表します。

タップ間の相関

WSSUS仮定のもとでは、異なるタップの係数は統計的に無相関です。

$$ \mathbb{E}\!\left[h_l^*[n] \, h_m[n]\right] = P_l \, \delta_{lm} $$

ここで $P_l = \mathbb{E}[|h_l[n]|^2]$ は第 $l$ タップの平均電力、$\delta_{lm}$ はクロネッカーのデルタです。

一方、同一タップの時間方向の相関は、ドップラースペクトルから決まるベッセル関数で記述されます。

$$ \mathbb{E}\!\left[h_l^*[n] \, h_l[n+m]\right] = P_l \, J_0(2\pi f_{\max} m T_s) $$

代表的なTDLモデル

国際電気通信連合（ITU）や 3GPP が標準化した代表的な TDL モデルがあります。

ITU Pedestrian A（歩行者環境）:

ITU Vehicular A（車両環境）:

Vehicular A は Pedestrian A に比べて最大遅延が大きく、かつ各タップの電力差も小さいことがわかります。これは車両環境のほうが遠方の散乱体が多く、遅延拡がりが大きいことを反映しています。

タップ遅延線モデルは、信号帯域幅に応じて「広帯域チャネル」と「狭帯域チャネル」に大別されます。次にこの分類について整理します。

広帯域チャネル vs 狭帯域チャネル

分類基準

マルチパスチャネルを通る信号が周波数選択性を経験するかどうかは、信号帯域幅 $W$ とコヒーレンス帯域幅 $B_c$ の関係（等価的に、シンボル時間 $T_s$ と遅延拡がり $\sigma_\tau$ の関係）で決まることはすでに述べました。この関係に基づいて、チャネルを「狭帯域」と「広帯域」に分類します。

狭帯域チャネル（narrowband channel）: $W \ll B_c$

TDL モデルでタップ数 $L = 1$ に対応します。チャネルは単一の複素乗算係数 $h[n]$ で表されます。

$$ y[n] = h[n] \, x[n] + w[n] $$

信号の波形は変化せず、受信電力が時間変動するフラットフェージングチャネルです。$h[n]$ の振幅がレイリー分布やライス分布に従い、時間相関がベッセル関数で与えられます。

広帯域チャネル（wideband channel）: $W \gg B_c$

TDL モデルのタップ数 $L > 1$ であり、チャネルは FIR（有限インパルス応答）フィルタとして振る舞います。信号の波形が歪み、ISI が発生するため、等化器やOFDMなどの対策が必要です。

二重分散チャネル

現実のマルチパスチャネルは、時間分散（遅延拡がり）と周波数分散（ドップラー拡がり）を同時に持ちます。これを二重分散チャネル（doubly dispersive channel）と呼びます。

チャネルの4つの特性パラメータと、それに対する判定基準をまとめると、以下のようになります。

これら2つの軸の組み合わせで、チャネルは4つのカテゴリに分類されます。

1. フラット・低速 ($W \ll B_c$, $T_s \ll T_c$) — 最も扱いやすいチャネル。単純な振幅変動のみ
2. フラット・高速 ($W \ll B_c$, $T_s \gg T_c$) — 波形は保たれるが1シンボル内でチャネルが変動
3. 周波数選択性・低速 ($W \gg B_c$, $T_s \ll T_c$) — ISI が生じるがチャネルは安定。等化器やOFDMで対処可能
4. 周波数選択性・高速 ($W \gg B_c$, $T_s \gg T_c$) — 最も過酷。ISI とチャネル変動の両方に対処が必要

実用上、多くの移動通信システムは「周波数選択性・低速」または「フラット・低速」のカテゴリに入るよう設計されています。5Gの高速移動シナリオでは「周波数選択性・高速」が問題になり、高速なチャネル推定アルゴリズムが研究されています。

ここまでの理論をPythonで実装し、マルチパスチャネルの振る舞いを可視化してみましょう。

Pythonによるシミュレーション

パワー遅延プロファイルとRMS遅延拡がり

まず、ITU Vehicular Aモデルのパワー遅延プロファイルを描画し、RMS遅延拡がりを計算します。

import numpy as np
import matplotlib.pyplot as plt

# ITU Vehicular A チャネルモデル
delays_ns = np.array([0, 310, 710, 1090, 1730, 2510])  # [ns]
powers_dB = np.array([0.0, -1.0, -9.0, -10.0, -15.0, -20.0])  # [dB]

# 線形スケールに変換
delays_us = delays_ns * 1e-3  # [μs]
powers_linear = 10 ** (powers_dB / 10)

# RMS 遅延拡がりの計算
P_total = np.sum(powers_linear)
mean_delay = np.sum(powers_linear * delays_us) / P_total
mean_delay_sq = np.sum(powers_linear * delays_us**2) / P_total
rms_delay_spread = np.sqrt(mean_delay_sq - mean_delay**2)

# コヒーレンス帯域幅
Bc = 1 / (5 * rms_delay_spread)  # [MHz]

print(f"平均遅延: {mean_delay:.4f} μs")
print(f"RMS 遅延拡がり: {rms_delay_spread:.4f} μs")
print(f"コヒーレンス帯域幅 (0.5相関): {Bc*1e3:.1f} kHz")

# 描画
fig, ax = plt.subplots(figsize=(9, 5))
markerline, stemlines, baseline = ax.stem(delays_us, powers_dB, basefmt=" ")
plt.setp(stemlines, linewidth=2, color='#00bcd4')
plt.setp(markerline, markersize=8, color='#00bcd4')
ax.set_xlabel("Excess Delay [μs]", fontsize=12)
ax.set_ylabel("Relative Power [dB]", fontsize=12)
ax.set_title("Power Delay Profile — ITU Vehicular A", fontsize=14)
ax.set_ylim(-25, 5)
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='gray', linewidth=0.5)

# RMS delay spread を矢印で表示
ax.annotate(f"σ_τ = {rms_delay_spread:.3f} μs\n"
            f"B_c ≈ {Bc*1e3:.0f} kHz",
            xy=(mean_delay, -12), fontsize=11,
            bbox=dict(boxstyle='round,pad=0.3', facecolor='#263238', edgecolor='#00bcd4', alpha=0.8),
            color='white')

plt.tight_layout()
plt.savefig("pdp_vehicular_a.png", dpi=150, bbox_inches='tight')
plt.show()

このコードを実行すると、6本の棒グラフとして ITU Vehicular A のパワー遅延プロファイルが描画されます。最初のタップ（0 ns）が最も強く、遅延が大きいタップほど電力が低下していることが確認できます。2番目のタップ（310 ns, -1.0 dB）は1番目とほぼ同じ電力を持っており、近距離の強い反射体が存在する環境を模擬しています。計算される RMS 遅延拡がりは約 0.37 $\mu$s、対応するコヒーレンス帯域幅は約 540 kHz となります。

マルチパスチャネルのインパルス応答シミュレーション

次に、TDLモデルに基づいてレイリーフェージングチャネルのインパルス応答を時間変動とともに生成します。

import numpy as np
import matplotlib.pyplot as plt

# チャネルパラメータ（ITU Vehicular A）
delays_ns = np.array([0, 310, 710, 1090, 1730, 2510])
powers_dB = np.array([0.0, -1.0, -9.0, -10.0, -15.0, -20.0])
powers_linear = 10 ** (powers_dB / 10)

# 正規化（全電力を1に）
powers_linear /= np.sum(powers_linear)

num_taps = len(delays_ns)
num_snapshots = 500  # 時間スナップショット数
f_max = 100  # 最大ドップラー周波数 [Hz]
Ts_channel = 1e-3  # チャネルサンプリング間隔 [s]

# 各タップのレイリーフェージング係数を生成（Jakes相関を近似）
np.random.seed(42)
N_sinusoids = 32  # ジェイクスモデルの正弦波数

tap_coefficients = np.zeros((num_taps, num_snapshots), dtype=complex)
time_axis = np.arange(num_snapshots) * Ts_channel

for l in range(num_taps):
    sigma = np.sqrt(powers_linear[l] / 2)
    # Sum-of-sinusoids (Jakes) 法
    h_l = np.zeros(num_snapshots, dtype=complex)
    for n in range(N_sinusoids):
        alpha_n = (2 * np.pi * n + np.random.uniform(0, 2*np.pi)) / N_sinusoids
        phi_n = np.random.uniform(0, 2 * np.pi)
        psi_n = np.random.uniform(0, 2 * np.pi)
        f_n = f_max * np.cos(alpha_n)
        h_l += np.exp(1j * (2 * np.pi * f_n * time_axis + phi_n))
    h_l *= sigma * np.sqrt(2.0 / N_sinusoids)
    tap_coefficients[l, :] = h_l

ここでは ITU Vehicular A の6タップに対して、Sum-of-Sinusoids 法を用いて時間変動するレイリーフェージング係数を生成しています。各タップの平均電力は PDP に従って設定されており、ドップラー効果に起因するランダムな時間変動を持ちます。次に、生成されたタップ係数の振幅変動と瞬時 PDP を可視化します。

# インパルス応答の3Dプロット
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左: タップごとの振幅の時間変動
for l in range(num_taps):
    amp_dB = 20 * np.log10(np.abs(tap_coefficients[l, :]) + 1e-10)
    axes[0].plot(time_axis * 1e3, amp_dB, alpha=0.7,
                 label=f"Tap {l} ({delays_ns[l]} ns)")

axes[0].set_xlabel("Time [ms]", fontsize=12)
axes[0].set_ylabel("Amplitude [dB]", fontsize=12)
axes[0].set_title("Tap Amplitude vs Time", fontsize=13)
axes[0].legend(fontsize=8, loc='lower left')
axes[0].set_ylim(-50, 5)
axes[0].grid(True, alpha=0.3)

# 右: ある時刻のインパルス応答（PDP のスナップショット）
snapshot_idx = [0, 100, 250]
colors = ['#00bcd4', '#ff9800', '#e91e63']
for i, idx in enumerate(snapshot_idx):
    amps = np.abs(tap_coefficients[:, idx]) ** 2
    amps_dB = 10 * np.log10(amps + 1e-30)
    axes[1].stem(delays_ns * 1e-3, amps_dB, linefmt=f'{colors[i]}-',
                 markerfmt=f'o', basefmt=' ',
                 label=f"t = {time_axis[idx]*1e3:.0f} ms")

axes[1].set_xlabel("Delay [μs]", fontsize=12)
axes[1].set_ylabel("Power [dB]", fontsize=12)
axes[1].set_title("Instantaneous PDP (snapshots)", fontsize=13)
axes[1].legend(fontsize=10)
axes[1].set_ylim(-55, 5)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("tdl_impulse_response.png", dpi=150, bbox_inches='tight')
plt.show()

左のグラフでは、6つのタップの振幅が時間とともにランダムに変動する様子が確認できます。Tap 0（直接波近傍）と Tap 1 が最も強く、Tap 5 は平均的に -20 dB 以上低いですが、フェージングにより瞬間的に他のタップと同程度まで強くなることもあります。深いフェージング（振幅が-30 dB以下に落ち込む「ディップ」）がランダムに発生していることも読み取れます。

右のグラフでは、3つの異なる時刻における瞬時のパワー遅延プロファイルを重ねて表示しています。同じチャネルモデルでも、時刻によってタップごとの電力配分が大きく異なることがわかります。これがフェージングの本質 — チャネルの時間変動性 — を視覚的に示しています。

ジェイクスドップラースペクトルの可視化

理論的なジェイクスドップラースペクトルと、シミュレーションで得たタップ係数のパワースペクトルを比較します。

import numpy as np
import matplotlib.pyplot as plt

f_max = 100  # Hz

# 理論的なジェイクススペクトル
f_theory = np.linspace(-f_max * 0.999, f_max * 0.999, 2000)
S_jakes = 1.0 / (np.pi * f_max * np.sqrt(1 - (f_theory / f_max)**2))

# シミュレーション（sum-of-sinusoidsで生成した時系列のPSD）
np.random.seed(42)
N_samples = 10000
Ts_ch = 1e-3
N_sin = 64
time_vec = np.arange(N_samples) * Ts_ch

h_sim = np.zeros(N_samples, dtype=complex)
for n in range(N_sin):
    alpha_n = (2 * np.pi * n + np.random.uniform(0, 2*np.pi)) / N_sin
    phi_n = np.random.uniform(0, 2 * np.pi)
    f_n = f_max * np.cos(alpha_n)
    h_sim += np.exp(1j * (2 * np.pi * f_n * time_vec + phi_n))
h_sim /= np.sqrt(N_sin)

# パワースペクトル密度の推定（ウェルチ法）
from scipy.signal import welch
fs_ch = 1.0 / Ts_ch  # 1000 Hz
freqs, psd = welch(h_sim, fs=fs_ch, nperseg=1024, noverlap=512,
                   return_onesided=False)

# 周波数を正しくソート
sort_idx = np.argsort(freqs)
freqs = freqs[sort_idx]
psd = psd[sort_idx]

# 正規化（ピーク値を合わせる）
psd_norm = psd / np.max(psd) * np.max(S_jakes) * 0.6

理論的なジェイクススペクトルは解析式から直接計算し、シミュレーション側はウェルチ法で PSD を推定しています。次にこれらを重ねてプロットし、スペクトル形状の一致を確認します。

fig, ax = plt.subplots(figsize=(9, 5))

ax.plot(f_theory, 10 * np.log10(S_jakes / np.max(S_jakes)),
        'w-', linewidth=2, label='Jakes theory', alpha=0.9)

mask = np.abs(freqs) < f_max * 1.5
ax.plot(freqs[mask], 10 * np.log10(psd_norm[mask] / np.max(psd_norm[mask]) + 1e-10),
        color='#00bcd4', linewidth=1.2, alpha=0.8, label='Simulation (Welch PSD)')

ax.axvline(x=f_max, color='#ff9800', linestyle='--', linewidth=1.5,
           label=f'f_max = {f_max} Hz')
ax.axvline(x=-f_max, color='#ff9800', linestyle='--', linewidth=1.5)

ax.set_xlabel("Doppler Frequency [Hz]", fontsize=12)
ax.set_ylabel("Normalized PSD [dB]", fontsize=12)
ax.set_title("Jakes Doppler Spectrum", fontsize=14)
ax.set_xlim(-200, 200)
ax.set_ylim(-30, 5)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("jakes_doppler_spectrum.png", dpi=150, bbox_inches='tight')
plt.show()

理論曲線（白線）は $f_d = \pm f_{\max}$ で鋭いピーク（特異点）を持つ典型的なバスタブ型を示しています。シミュレーション結果（シアン線）は有限長のデータからウェルチ法で推定しているため、ピーク付近が滑らかになっていますが、全体的な形状は理論曲線とよく一致しています。$|f_d| > f_{\max}$ の領域ではスペクトルが急激に減衰し、ドップラーシフトが $f_{\max}$ を超えないことが確認できます。

周波数応答と時間変動の可視化

最後に、マルチパスチャネルの周波数応答の時間変動を2Dカラーマップで可視化します。これにより、チャネルが周波数軸と時間軸の両方で変動する「二重分散」の様子を直感的に捉えることができます。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
delays_ns = np.array([0, 310, 710, 1090, 1730, 2510])
powers_dB = np.array([0.0, -1.0, -9.0, -10.0, -15.0, -20.0])
powers_linear = 10 ** (powers_dB / 10)
powers_linear /= np.sum(powers_linear)

num_taps = len(delays_ns)
num_time = 200
f_max = 100
Ts_ch = 1e-3
N_sin = 32
num_freq = 256

np.random.seed(123)
time_vec = np.arange(num_time) * Ts_ch

# タップ係数の生成
tap_coeff = np.zeros((num_taps, num_time), dtype=complex)
for l in range(num_taps):
    sigma = np.sqrt(powers_linear[l] / 2)
    h_l = np.zeros(num_time, dtype=complex)
    for n in range(N_sin):
        alpha_n = (2 * np.pi * n + np.random.uniform(0, 2*np.pi)) / N_sin
        phi_n = np.random.uniform(0, 2 * np.pi)
        f_n = f_max * np.cos(alpha_n)
        h_l += np.exp(1j * (2 * np.pi * f_n * time_vec + phi_n))
    h_l *= sigma * np.sqrt(2.0 / N_sin)
    tap_coeff[l, :] = h_l

# 周波数応答の計算
freq_axis = np.linspace(-2, 2, num_freq)  # [MHz]
H = np.zeros((num_freq, num_time), dtype=complex)

for t_idx in range(num_time):
    for l in range(num_taps):
        tau_us = delays_ns[l] * 1e-3  # μs → MHz の逆数
        H[:, t_idx] += tap_coeff[l, t_idx] * np.exp(-1j * 2 * np.pi * freq_axis * tau_us)

各時刻について、タップ係数と遅延に基づいてチャネルの周波数応答 $H(f, t)$ を計算しました。次にこの $|H(f, t)|$ を 2D カラーマップとして描画し、時間・周波数の両方向での変動（二重分散性）を視覚的に確認します。

# 2Dカラーマップ
fig, ax = plt.subplots(figsize=(10, 6))

H_dB = 20 * np.log10(np.abs(H) + 1e-10)
im = ax.pcolormesh(time_vec * 1e3, freq_axis, H_dB,
                   shading='auto', cmap='inferno', vmin=-25, vmax=5)

cbar = fig.colorbar(im, ax=ax, label='|H(f,t)| [dB]')
ax.set_xlabel("Time [ms]", fontsize=12)
ax.set_ylabel("Frequency [MHz]", fontsize=12)
ax.set_title("Time-Varying Channel Frequency Response", fontsize=14)

# コヒーレンス帯域幅のスケールバー
sigma_tau = 0.37  # μs (Vehicular A の近似値)
Bc_MHz = 1 / (5 * sigma_tau)
ax.annotate('', xy=(5, -1.5), xytext=(5, -1.5 + Bc_MHz),
            arrowprops=dict(arrowstyle='<->', color='#00bcd4', lw=2))
ax.text(8, -1.5 + Bc_MHz/2, f'B_c ≈ {Bc_MHz:.2f} MHz',
        color='#00bcd4', fontsize=10, va='center')

plt.tight_layout()
plt.savefig("channel_frequency_response_2d.png", dpi=150, bbox_inches='tight')
plt.show()

この2Dカラーマップでは、横軸が時間、縦軸が周波数を表し、色がチャネルの振幅を表しています。以下の特徴が読み取れます。

1. 周波数軸方向の変動 — 同一時刻で見ると、周波数によって振幅が大きく異なります。明るい領域（高い振幅）と暗い領域（深いフェード）が周波数方向に交互に現れています。この変動の「幅」がコヒーレンス帯域幅 $B_c$ に対応します
2. 時間軸方向の変動 — 同一周波数で見ると、時間とともに振幅が変動しています。明暗のパターンが時間とともにゆるやかにシフトする様子が確認できます。この変動の「速さ」がコヒーレンス時間 $T_c$ で特徴づけられます
3. 二重分散 — 周波数と時間の両方でチャネルが変動しており、「明暗の島」が2次元的に分布しています。これが二重分散チャネルの典型的な振る舞いです

時間相関関数とベッセル関数の比較

理論で導いた $R_h(\Delta t) = J_0(2\pi f_{\max} \Delta t)$ をシミュレーション結果と比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j0

f_max = 100  # Hz
Ts_ch = 1e-3
num_time = 5000
N_sin = 64

np.random.seed(77)
time_vec = np.arange(num_time) * Ts_ch

# 単一タップのレイリーフェージング系列を生成
h_sim = np.zeros(num_time, dtype=complex)
for n in range(N_sin):
    alpha_n = (2 * np.pi * n + np.random.uniform(0, 2*np.pi)) / N_sin
    phi_n = np.random.uniform(0, 2 * np.pi)
    f_n = f_max * np.cos(alpha_n)
    h_sim += np.exp(1j * (2 * np.pi * f_n * time_vec + phi_n))
h_sim /= np.sqrt(N_sin)

# 自己相関の計算
max_lag = 200
lags = np.arange(max_lag)
R_sim = np.zeros(max_lag, dtype=complex)
for m in range(max_lag):
    R_sim[m] = np.mean(np.conj(h_sim[:num_time-max_lag]) * h_sim[m:num_time-max_lag+m])
R_sim /= R_sim[0]  # 正規化

# 理論値
delta_t = lags * Ts_ch
R_theory = j0(2 * np.pi * f_max * delta_t)

# コヒーレンス時間
Tc = 1 / (4 * f_max)  # [s]

シミュレーション結果を $R_h(0)$ で正規化しているため、$\Delta t = 0$ で値が1となる相関係数として理論値と比較できます。次に理論曲線のベッセル関数 $J_0(2\pi f_{\max} \Delta t)$ と重ねてプロットします。

fig, ax = plt.subplots(figsize=(9, 5))

ax.plot(delta_t * 1e3, np.real(R_sim), color='#00bcd4', linewidth=1.5,
        alpha=0.8, label='Simulation')
ax.plot(delta_t * 1e3, R_theory, 'w--', linewidth=2, alpha=0.9,
        label=r'Theory: $J_0(2\pi f_{max} \Delta t)$')

ax.axhline(y=0.5, color='#ff9800', linestyle=':', linewidth=1.5, alpha=0.7)
ax.axvline(x=Tc * 1e3, color='#e91e63', linestyle='--', linewidth=1.5,
           label=f'T_c = 1/(4f_max) = {Tc*1e3:.1f} ms')

ax.set_xlabel(r"$\Delta t$ [ms]", fontsize=12)
ax.set_ylabel(r"Normalized $R_h(\Delta t)$", fontsize=12)
ax.set_title("Time Autocorrelation of Rayleigh Fading Channel", fontsize=14)
ax.legend(fontsize=11)
ax.set_xlim(0, delta_t[-1] * 1e3)
ax.set_ylim(-0.5, 1.1)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("time_autocorrelation.png", dpi=150, bbox_inches='tight')
plt.show()

シミュレーション結果（シアン実線）は理論曲線のベッセル関数 $J_0$（白破線）と極めてよく一致しています。相関関数は $\Delta t = 0$ で1（完全相関）から出発し、振動しながら減衰していきます。コヒーレンス時間 $T_c = 1/(4f_{\max}) = 2.5$ ms（ピンク破線）の付近で相関が十分に低下していることが確認でき、この定義が実用的に妥当であることが視覚的にもわかります。黄色の点線は相関 0.5 のラインを示しており、$J_0$ がこのラインを横切る時刻は $T_c$ よりやや小さい値であることも読み取れます。

まとめ

本記事では、マルチパス伝搬の物理メカニズムから定量的なチャネルモデルまでを体系的に解説しました。

• マルチパス伝搬とは、反射・回折・散乱により同じ信号が複数の経路で受信点に到達する現象であり、各経路の遅延・振幅・位相の違いが干渉を引き起こします
• RMS遅延拡がり $\sigma_\tau$ はパワー遅延プロファイルの標準偏差であり、マルチパスの時間的広がりを定量化します。これと信号帯域幅の関係で周波数選択性フェージングかフラットフェージングかが決まります
• コヒーレンス帯域幅 $B_c \approx 1/(5\sigma_\tau)$ は周波数相関関数がパワー遅延プロファイルのフーリエ変換であることから導出されます
• 最大ドップラー周波数 $f_{\max} = v/\lambda$ は移動速度と搬送波周波数で決まり、ジェイクスモデルはバスタブ型のドップラースペクトルを予測します
• コヒーレンス時間 $T_c \approx 1/(4f_{\max})$ はシンボル時間との比較で高速フェージングか低速フェージングかを判定します。時間相関関数はベッセル関数 $J_0$ で表されます
• タップ遅延線モデルはマルチパスチャネルを離散時間 FIR フィルタとして表現し、各タップの振幅がレイリー分布やライス分布に従います
• 遅延拡がり（周波数軸）とドップラー拡がり（時間軸）の2つの軸で二重分散チャネルを分類することで、通信システムの設計指針が得られます

マルチパスチャネルの理解は、そのまま通信システム設計に直結します。次のステップとして、以下の記事も参考にしてください。

• OFDMの原理と実装 — 周波数選択性フェージングに対抗する OFDM の仕組み
• 等化器の設計 — マルチパスによる ISI を除去する適応等化技術