ドップラースプレッドとコヒーレンス時間 — 時間選択性フェージングの尺度

歩きながらスマートフォンで通話していると、数十ミリ秒おきに音がわずかに途切れることがあります。歩みを止めると途切れが消え、逆に電車や車で移動すると途切れの頻度が上がる — この「移動すると通信品質が細かく揺れる」現象は、電波の強さが時間とともに激しく変動する時間選択性フェージングが原因です。同じ場所に立ち止まっていれば電界強度は安定しているのに、わずか数センチ動いただけで受信レベルが 20 dB 落ちることさえあります。では、その揺れは「どれくらいの速さ」で起きるのでしょうか。この速さを定量化する尺度が、本記事の主役であるドップラースプレッド $B_D$ とコヒーレンス時間 $T_c$ です。

この2つの量を理解すると、無線システム設計の具体的な判断が数式で下せるようになります。

• チャネル推定とパイロット設計: 受信機はパイロット信号（既知シンボル）でチャネルを測りますが、測った値が「いつまで有効か」を決めるのがコヒーレンス時間です。パイロット間隔が $T_c$ を超えると推定は一気に破綻します
• CSI フィードバックと適応変調（AMC）: 端末が測ったチャネル状態を基地局に返して変調方式を切り替える仕組みは、往復遅延が $T_c$ より短いことを前提にしています。高速移動時に AMC が機能しなくなるのはこのためです
• OFDM のサブキャリア間干渉（ICI）: 1 OFDM シンボルの間にチャネルが動くと、サブキャリアの直交性が崩れて干渉が生じます。その大きさは $f_d T$（正規化ドップラー）だけで決まります
• レーダー・ソナーの積分時間: 目標からの反射をコヒーレント積分できる時間の上限も、同じドップラー広がりで決まります

つまりコヒーレンス時間は「チャネルの賞味期限」です。この記事では、移動体まわりに散乱体がぐるりと分布する状況からドップラースペクトルの形（いわゆるバスタブ形）を変数変換で導出し、そのフーリエ変換が第 1 種 0 次ベッセル関数 $J_0$ になること、そこから $T_c \approx 0.423/f_d$ という有名な目安が出てくることを、途中式を省かずに追いかけます。最後に Python でフェージング波形を作り、スペクトル・相関・コヒーレンス時間・パイロット間隔と推定誤差の関係を実測して、理論式が本当に成り立っているかを確認します。

全体像を先に絵で掴んでおきましょう。次の図は、これから数式で追いかける話の骨格をそのまま描いたものです。

左の絵では、移動する受信機に全方位から波が到来しています。進行方向の前から来る波（赤）は周波数が最大 $+f_d$ だけ上がり、後ろから来る波（青）は $-f_d$ だけ下がり、真横から来る波（灰）はまったくシフトしません。右のグラフはその結果として現れるスペクトルで、真正面・真後ろ付近の広い角度範囲の電力が両端 $\pm f_d$ に押し込まれるため、真ん中がへこみ両端が跳ね上がる「バスタブ形」になります。この左から右への変換（角度の分布 → 周波数の分布）を数式で書き下すことが、本記事前半の目標です。

本記事の内容

• なぜ移動すると受信レベルが激しく揺れるのか（マルチパスとドップラーの直感）
• 最大ドップラー周波数 $f_d$ と到来角の関係 $f = f_d\cos\theta$ の導出
• 2 次元等方散乱を仮定したドップラースペクトル（Clarke/Jakes スペクトル）の導出（ヤコビアンつき）
• RMS ドップラースプレッド $\sigma_D = f_d/\sqrt{2}$ の計算
• スペクトルの逆フーリエ変換が $J_0(2\pi f_d \tau)$ になることの導出
• コヒーレンス時間の定義と $T_c \approx 9/(16\pi f_d) \approx 0.179/f_d$、幾何平均定義 $T_c \approx 0.423/f_d$
• 遅延スプレッド／コヒーレンス帯域幅との双対関係とフェージングの 4 分類
• レベルクロス率と平均フェード時間
• Python 実装: フェージング波形生成、スペクトル推定、自己相関と $J_0$ の一致、$T_c$ の実測、パイロット間隔と推定誤差

前提知識

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

• ドップラー効果を詳しく理解する
• マルチパス伝搬 — 反射・回折・散乱がつくる受信波の重ね合わせ
• Jakes モデルによるフェージングシミュレーション

なぜ移動すると受信レベルが揺れるのか

まず、時間選択性フェージングがどこから生まれるのかを、できるだけ単純な状況で掴んでおきましょう。

都市部で電波を受けるとき、受信アンテナに届くのは送信アンテナからの 1 本の波ではありません。建物の壁、地面、車体、街路樹などで反射・散乱された無数の波が、それぞれ別々の経路長を通ってやってきます。経路長が違えば位相も違うので、受信点での合成波は「たくさんのベクトルを足し合わせた結果」になります。ベクトルの向きがそろえば大きな振幅（ピーク）、打ち消し合えば深い落ち込み（フェードやヌル）です。

ここで受信機が動くと何が起こるでしょうか。受信点が $\Delta d$ だけ動くと、進行方向の前方から来る波は経路長が $\Delta d$ 短くなり、後方から来る波は $\Delta d$ 長くなります。位相にすると $2\pi \Delta d/\lambda$ の増減です。波長 $\lambda$ が 2 GHz で 15 cm しかないことを思い出してください。わずか 7.5 cm 動くだけで、前方波と後方波の位相差は $2\pi$ も変わります。合成ベクトルの向きが総取り替えになるので、受信レベルはガラリと変わります。これが「数センチ動くと 20 dB 落ちる」現象の正体です。

もっとも単純化して、前方（進行方向、$\theta = 0$）から来る波と後方（$\theta = \pi$）から来る波の 2 波だけを考えてみましょう。速度 $v$ で動くと、前方波は $+f_d = v/\lambda$ だけ周波数が上がり、後方波は $-f_d$ だけ下がります。搬送波を取り除いたベースバンドで見ると、2 つの複素正弦波 $e^{+j2\pi f_d t}$ と $e^{-j2\pi f_d t}$ の和です。その振幅は $|2\cos(2\pi f_d t)|$ となり、周期 $1/(2f_d)$ のうなり（ビート）になります。時速 60 km、2 GHz なら $f_d = 111$ Hz なので、4.5 ms ごとに深いフェードが訪れる計算です。

上段の黒い太線が合成波の振幅 $|2\cos(2\pi f_d t)|$ で、灰色の細い線は搬送波を含めた瞬時値（見やすさのため搬送波周波数を大幅に落とした模式）です。振幅がゼロを横切るたびに 2 波が完全に打ち消し合っており、その間隔がちょうど $1/(2f_d) = 4.50$ ms になっています。下段の dB 表示にすると、この打ち消しが「底なしに落ち込むヌル」として現れることがはっきり見え、実際の受信機が経験する深いフェードの正体がこれと同じ機構であることがわかります。

現実には到来方向は 2 つではなく連続的に分布しますから、うなりの周期はきれいには決まらず、レベル変動はランダムになります。それでも「変動の時間スケールが $1/f_d$ 程度」という骨格は変わりません。ドップラーによる周波数の広がりが大きいほど、時間変動が速い — この対応関係を厳密に定式化するのが、次以降で導くドップラースペクトルとその相関関数です。

まずは、到来角と周波数のずれを結ぶ基本式 $f = f_d\cos\theta$ をきちんと書き下しておきましょう。

最大ドップラー周波数と到来角

幾何から $f = f_d\cos\theta$ を導く

移動体が速度ベクトル $\bm{v}$（大きさ $v$）で動いているとします。ある平面波が、単位方向ベクトル $\hat{\bm{k}}$ の向きに進んで移動体に入射しているとしましょう。移動体から見ると、この波の波面は「移動体の速度の、波の進行方向成分だけ」余分に近づいたり遠ざかったりします。

この余分な接近速度は $-\bm{v}\cdot\hat{\bm{k}}$ です（波の進行方向と同じ向きに逃げていれば波面から遠ざかるのでマイナス）。到来方向、すなわち「波が来る方向」を進行方向から測った角度 $\theta$ で表すと、$\bm{v}\cdot\hat{\bm{k}} = -v\cos\theta$ となるので、接近速度は $v\cos\theta$ です。単位時間に $v\cos\theta$ だけ余分に波面を横切るので、波長 $\lambda$ あたり 1 周期と数えれば、観測される周波数のずれは次のようになります。

$$ \begin{equation} f_{\text{Doppler}} = \frac{v\cos\theta}{\lambda} = f_d\cos\theta, \qquad f_d \equiv \frac{v}{\lambda} = \frac{v f_c}{c} \end{equation} $$

ここで $f_c$ は搬送波周波数、$c$ は光速、$f_d$ を最大ドップラー周波数と呼びます。$\theta = 0$（真正面から来る波）で $+f_d$、$\theta = \pi$（真後ろから追いかけてくる波）で $-f_d$、$\theta = \pi/2$（真横）でゼロです。$v \ll c$ なので相対論補正は不要で、この 1 次の式で十分です。

大切なのは、ドップラーシフトが到来角にしか依存しないことです。散乱体がどこにあろうと、受信点で見た到来角が同じなら同じ周波数ずれを受けます。したがって「到来角の分布」が決まれば「周波数ずれの分布」、つまりスペクトルの形が決まります。この論理が次節の導出の核心です。

数値感覚をつける

$f_d = v f_c/c$ は、速度と搬送波周波数の積に比例します。移動体通信で「高い周波数帯ほど高速移動に弱い」と言われるのは、この式が理由です。代表的な値を並べてみましょう（$c = 3.0\times 10^8$ m/s とします）。

表の右端の $T_c$ はこの後で導く目安値です。歩行者が 800 MHz 帯を使うときの $T_c$ は 100 ms を超え、チャネルは「ほぼ静止」と見なせます。一方、ミリ波（28 GHz）を車の中で使うと、たかだか時速 60 km でも $T_c$ は 0.3 ms 弱しかありません。同じ速度でも、周波数を 14 倍にすればチャネルの寿命は 1/14 になるのです。60 GHz で歩行する場合の $f_d$ が、2 GHz で時速 150 km 走るのと同程度だという事実は、ミリ波屋内システムの設計者にとって重要な直感です。

なお、LEO 衛星のように毎秒 7.5 km で動く場合、2 GHz では $f_d$ が 50 kHz にも達します。ただしこのケースでは、地上のマルチパス散乱と違って視線方向のはっきりした 1 本の波が支配的なので、生じるのは「広がり」ではなく大きな「シフト」です。シフトは周波数プリコンペンセーションで打ち消せますが、広がりは打ち消せません。両者を区別することが実務では大事になります。

両対数で描くと、$f_d = v f_c/c$ も $T_c \approx 0.423/f_d$ もすべて傾き $\pm 1$ の平行な直線群になります。黒丸は上の表の代表例で、たとえば 60 GHz で歩行する点（左上の明るい緑）が、2 GHz で時速 150 km 走る点とほぼ同じ高さにあることが読み取れます。右のグラフで赤い破線（1 ms）を下回る領域は「チャネルが LTE のサブフレーム 1 個ぶんも持たない」危険地帯で、ミリ波帯では歩行速度でもすぐそこまで来ていることがわかります。

到来角と周波数の対応がついたので、いよいよ「到来角が一様分布するとき、周波数の分布はどうなるか」を計算しましょう。

ドップラースペクトルの導出（Clarke/Jakes スペクトル）

モデルの仮定

Clarke が 1968 年に提案し、Jakes が 1974 年の教科書で広めたモデルは、市街地の移動体受信を次のように理想化します。

1. 受信点のまわり 水平面内に一様に散乱体が分布し、多数の平面波が全方位から到来する（2 次元等方散乱）
2. 各到来波の振幅はほぼ等しく、位相は互いに独立で $[0, 2\pi)$ に一様
3. 受信アンテナは水平面内で無指向性（利得 $G$ が $\theta$ によらず一定）
4. 直接波（見通し成分）は存在しない

仮定 2 と多数波の重ね合わせから、中心極限定理によりベースバンドの複素チャネル利得 $g(t)$ は複素ガウス過程になり、その包絡線 $|g(t)|$ はレイリー分布に従います。ここまではレイリーフェージングの標準的な話です。本記事で追いかけたいのは「その複素ガウス過程が、時間方向にどんな色をもつか」、つまりパワースペクトル密度の形です。

仮定 1 を数式にすると、到来角 $\theta$ の確率密度関数は次のようになります。

$$ p_\Theta(\theta) = \frac{1}{2\pi}, \qquad -\pi \le \theta < \pi $$

受信総電力を $P$（以下では $P = 1$ に正規化）とすると、微小角 $d\theta$ から到来する電力は $P\,p_\Theta(\theta)\,d\theta$ です。この電力が、周波数軸上では $f = f_d\cos\theta$ の位置に置かれます。「角度の上に分布していた電力を、周波数の上に載せ替える」— これが変数変換です。

変数変換とヤコビアン

確率変数 $\Theta$ の関数 $F = f_d\cos\Theta$ の密度を求めます。単調でない変換なので、逆像が複数ある点に注意が必要です。$f \in (-f_d, f_d)$ を固定すると、$\cos\theta = f/f_d$ を満たす $\theta$ は区間 $[-\pi, \pi)$ に 2 つあります。

$$ \theta_{1,2} = \pm\arccos\!\left(\frac{f}{f_d}\right) $$

進行方向の右側から来る波と左側から来る波が、同じドップラーシフトを与えるということです。変数変換の公式は、逆像すべてにわたる和になります。

$$ S(f) = \sum_{i=1,2} \frac{p_\Theta(\theta_i)}{\left|\dfrac{df}{d\theta}\right|_{\theta=\theta_i}} $$

ヤコビアン（1 変数なので単なる導関数）を計算しましょう。$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 = 1 – (f/f_d)^2$ であり、$|\sin\theta| = \sqrt{1 – (f/f_d)^2}$ です。したがって、

$$ \left|\frac{df}{d\theta}\right| = f_d\sqrt{1 – \left(\frac{f}{f_d}\right)^2} $$

この値は $\theta_1$ と $\theta_2$ のどちらでも同じ（$|\sin\theta|$ は偶関数）です。$p_\Theta(\theta_i) = 1/(2\pi)$ も共通なので、和は単に 2 倍になります。以上を代入すると、

$$ S(f) = 2\cdot\frac{1/(2\pi)}{f_d\sqrt{1 – (f/f_d)^2}} $$

分子の 2 と $2\pi$ の 2 が約分されて、有名なClarke（Jakes）のドップラースペクトルが得られます。

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

$|f| \ge f_d$ では $S(f) = 0$ です。物理的に当然で、$|\cos\theta| \le 1$ なので $f_d$ を超えるドップラーシフトは起こりえません。

$f_d = 100$ Hz として描いたのがこの曲線です。中央での値は $1/(\pi f_d) = 0.00318$ ちょうどで、そこから両端に向かってじわじわ持ち上がり、$\pm 100$ Hz の直前で急激に立ち上がります。水色で塗った面積が全電力 1 に等しく、両端で発散しているのに面積が有限にとどまっている点が、次に確認する「可積分な特異点」の意味です。オレンジの点線で示した $\pm\sigma_D = \pm 70.7$ Hz は後述する RMS ドップラースプレッドの位置で、電力の「重心からの広がり」がスペクトル端の 7 割ほどのところにあることを表しています。

正規化の確認

導出が正しいかは、全電力が 1 に戻るかで検算できます。$u = f/f_d$ と置くと $df = f_d\,du$ なので、

$$ \int_{-f_d}^{f_d} S(f)\,df = \frac{1}{\pi f_d}\int_{-1}^{1}\frac{f_d\,du}{\sqrt{1-u^2}} = \frac{1}{\pi}\Big[\arcsin u\Big]_{-1}^{1} = \frac{1}{\pi}\left(\frac{\pi}{2} + \frac{\pi}{2}\right) = 1 $$

$\pm f_d$ で被積分関数が発散しているにもかかわらず、積分は有限値 1 に収束します。これは $1/\sqrt{1-u^2}$ という可積分な特異点（べき $-1/2$ の発散）だからです。数値計算で端の点をそのまま評価すると無限大になるので、実装では端を少しだけ内側に寄せる必要があります。

なぜ「バスタブ形」なのか

$S(f)$ は $f = 0$ で最小値 $1/(\pi f_d)$、$f \to \pm f_d$ で発散します。真ん中がへこみ両端が立ち上がる形なので、俗にバスタブ（U 字）スペクトルと呼ばれます。この形は直感的にも理解できます。

鍵はヤコビアン $|df/d\theta| = f_d|\sin\theta|$ です。$\theta = \pi/2$（真横）付近では $|\sin\theta| = 1$ で最大、つまり「角度を少し変えると周波数が大きく変わる」ため、電力は周波数軸上に広くばらまかれ、密度は薄くなります。逆に $\theta = 0$ や $\pi$（真正面・真後ろ）付近では $\cos\theta$ が停留点をもち、$|\sin\theta| \to 0$ になります。角度が大きく変わっても周波数がほとんど動かないので、広い角度範囲の電力が $\pm f_d$ のごく狭い周波数帯に押し込められ、密度が発散するのです。

同じ現象は光学や分光学でも見られます。等方的に運動する粒子の吸収線を 1 方向から見たときのプロファイル、あるいは回転する円板上の点の速度を 1 軸に射影したときの分布 — いずれも「射影の停留点で密度が発散する」構造で、$1/\sqrt{1-u^2}$ はアークサイン分布（第 1 種ベータ分布 $\mathrm{Beta}(1/2,1/2)$）そのものです。

(a) を見ると、周波数軸上で同じ幅 $0.1f_d$ を切り取っても、それに対応する到来角の範囲は端（赤帯）では合計 $52^\circ$ もあるのに、中央（緑帯）では $11^\circ$ しかありません。この差がそのまま密度の差になります。(b) はその倍率であるヤコビアン $|\sin\theta|$ で、$\theta = 0, \pi$ で完全にゼロへ落ちること（写像が潰れること）が発散の直接の原因だとわかります。(c) では円周上に等間隔に並べた到来方向を実際に周波数軸へ射影しており、点が両端にぎゅっと密集してヒストグラムが理論曲線に一致する様子が確認できます。

非等方の場合への一般化

実際のアンテナには指向性があり、到来角も一様とは限りません。アンテナの水平面利得を $G(\theta)$、到来角密度を $p_\Theta(\theta)$ とすると、同じ変数変換の手順で次の一般形が得られます。

$$ S(f) = \frac{G(\theta_1)p_\Theta(\theta_1) + G(\theta_2)p_\Theta(\theta_2)}{f_d\sqrt{1 – (f/f_d)^2}}, \qquad \theta_{1,2} = \pm\arccos(f/f_d) $$

分母のヤコビアン因子 $\sqrt{1-(f/f_d)^2}$ は変わらないので、両端が立ち上がる性質はアンテナや散乱分布によらず残ります。変わるのは分子の重みだけです。たとえばビームフォーミングで進行方向にビームを絞れば $\theta \approx 0$ の成分だけが残り、スペクトルは $+f_d$ 付近に集中した「片側だけのバスタブ」になります。ミリ波システムで狭ビームを使うと実効的なドップラー広がりが減るのは、この効果です。

左が仮定した到来角分布、右が対応するスペクトルです。等方散乱（青）では両端が同じ高さになりますが、進行方向に偏らせる（橙、赤）につれて電力が $+f_d$ 側だけに集まり、$-f_d$ 側はほとんど消えます。それでも $+f_d$ の直前で立ち上がる形は 3 つとも共通で、これが「ヤコビアン因子は分布によらない」という主張の可視化です。赤い鋭いビームの場合、電力の大半が $+f_d$ のごく近くに集中するため、実効的な広がり（スペクトルの実質的な幅）は等方散乱よりずっと小さくなります。

RMS ドップラースプレッド

スペクトルの「広がり」を 1 つの数字で表すには、遅延プロファイルの RMS 遅延広がりと同じ発想で、2 次モーメントを使います。平均ドップラーシフト $\bar{f}$ と RMS ドップラースプレッド $\sigma_D$ は次のように定義されます。

$$ \bar{f} = \frac{\displaystyle\int f S(f)\,df}{\displaystyle\int S(f)\,df}, \qquad \sigma_D = \sqrt{\frac{\displaystyle\int (f-\bar{f})^2 S(f)\,df}{\displaystyle\int S(f)\,df}} $$

Clarke スペクトルは偶関数なので $\bar{f} = 0$ です。2 次モーメントを計算しましょう。$f = f_d\cos\theta$ に戻して積分するのが簡単です。分母は 1 なので分子だけ計算すると、

$$ \int_{-f_d}^{f_d} f^2 S(f)\,df = \mathbb{E}\!\left[f_d^2\cos^2\Theta\right] = f_d^2\,\mathbb{E}[\cos^2\Theta] $$

$\Theta$ が一様分布なら $\mathbb{E}[\cos^2\Theta] = 1/2$ なので、

$$ \begin{equation} \sigma_D = \frac{f_d}{\sqrt{2}} \approx 0.707\,f_d \end{equation} $$

「ドップラースプレッド $B_D$」という言葉は文脈によって $f_d$ そのもの（スペクトルの片側幅）、$2f_d$（全幅）、$\sigma_D$（RMS 値）のどれを指すこともあります。論文や規格書を読むときは定義を確認する癖をつけてください。本記事では以降、混乱を避けるため $f_d$ を基準にして議論します。

スペクトルの形が決まりました。次は、このスペクトルが時間領域でどんな相関を意味するのかを、フーリエ変換で調べます。

時間自己相関はベッセル関数になる

ウィーナー・ヒンチンの定理から出発する

パワースペクトル密度と自己相関関数はフーリエ変換のペアです（ウィーナー・ヒンチンの定理）。複素チャネル利得 $g(t)$ の自己相関関数を

$$ R_g(\tau) = \mathbb{E}\!\left[g^*(t)\,g(t+\tau)\right] $$

と定義すると、$S(f)$ の逆フーリエ変換が $R_g(\tau)$ になります。

$$ R_g(\tau) = \int_{-\infty}^{\infty} S(f)\,e^{j2\pi f\tau}\,df = \int_{-f_d}^{f_d} \frac{e^{j2\pi f\tau}}{\pi f_d\sqrt{1-(f/f_d)^2}}\,df $$

この積分は、素直に $f = f_d\cos\theta$ という置換で計算できます。導出時と逆向きに戻すわけです。

置換して積分表示に持ち込む

$f = f_d\cos\theta$ とすると $df = -f_d\sin\theta\,d\theta$ であり、積分区間は $f: -f_d \to f_d$ が $\theta: \pi \to 0$ に対応します。また分母の平方根は $\sqrt{1-\cos^2\theta} = \sin\theta$（$0\le\theta\le\pi$ で $\sin\theta \ge 0$）です。これらを代入すると、

$$ R_g(\tau) = \int_{\pi}^{0} \frac{e^{j2\pi f_d\tau\cos\theta}}{\pi f_d \sin\theta}\,(-f_d\sin\theta)\,d\theta $$

マイナス符号で積分の上下限を入れ替え、$f_d\sin\theta$ が約分されることに注目すると、非常にすっきりした形になります。

$$ R_g(\tau) = \frac{1}{\pi}\int_{0}^{\pi} e^{j2\pi f_d\tau\cos\theta}\,d\theta $$

ここで、第 1 種 0 次ベッセル関数の積分表示（ベッセルの積分）

$$ J_0(z) = \frac{1}{\pi}\int_{0}^{\pi} e^{jz\cos\theta}\,d\theta = \frac{1}{\pi}\int_{0}^{\pi}\cos(z\cos\theta)\,d\theta $$

を使えば、$z = 2\pi f_d\tau$ として次の結論が得られます。

$$ \begin{equation} R_g(\tau) = J_0(2\pi f_d \tau) \end{equation} $$

きわめて美しい結果です。バスタブ形という一見扱いにくいスペクトルが、時間領域では 0 次ベッセル関数という 1 個の特殊関数にまとまってしまいます。ちなみに虚部が消えるのは、$e^{jz\cos\theta}$ の虚部 $\sin(z\cos\theta)$ が $\theta \to \pi – \theta$ の置換で符号反転するため、$0$ から $\pi$ の積分でちょうど打ち消し合うからです。つまり $R_g(\tau)$ は実関数で、しかも偶関数です。

$J_0$ の性質から読み取れること

$J_0(z)$ のふるまいを押さえておくと、フェージングの時間構造がそのまま見えてきます。

(a) 原点付近の 2 次的な落ち込み: $J_0$ の級数展開は

$$ J_0(z) = 1 – \frac{z^2}{4} + \frac{z^4}{64} – \cdots $$

なので、$\tau$ が小さいところでは

$$ R_g(\tau) \approx 1 – \frac{(2\pi f_d\tau)^2}{4} = 1 – (\pi f_d\tau)^2 $$

と放物線的に落ちます。相関が「1 のところから水平に出発して二次で落ちる」形は、チャネル予測（channel prediction）が短時間なら高精度でできることを意味します。1 次で落ちる（微分不可能な）過程だったら予測は絶望的です。

(b) 最初のゼロ点: $J_0(z) = 0$ の最小の正の解は $z \approx 2.4048$ です。したがって、

$$ \tau_0 = \frac{2.4048}{2\pi f_d} \approx \frac{0.3827}{f_d} $$

で相関がちょうどゼロになります。時速 60 km、2 GHz なら $f_d = 111$ Hz なので $\tau_0 \approx 3.4$ ms です。「3.4 ms 離れた 2 時点のチャネルは無相関」という、時間ダイバーシチ設計に直結する数字が出てきました。

(c) 減衰しながらの振動: 大きな $z$ では

$$ J_0(z) \approx \sqrt{\frac{2}{\pi z}}\cos\!\left(z – \frac{\pi}{4}\right) $$

と、$1/\sqrt{z}$ で減衰しながら振動します。相関がゼロを何度も横切り、負の値もとることに注意してください。これは「$\pm f_d$ に強いスペクトル成分がある」ことの時間領域での表れで、前節で 2 波モデルが $\cos(2\pi f_d t)$ のうなりを生んだのと同じ現象です。相関が完全にゼロへ落ちきらず振動し続けるのは、チャネル予測アルゴリズムにとってはむしろ好都合で、遠い過去の観測にも情報が残っていることを意味します。

(a) では、原点で水平に出発して二次で落ちる様子と、その近似 $1-(\pi f_d\tau)^2$（橙の破線）が $f_d\tau \lesssim 0.15$ までよく一致することが見て取れます。近似のほうが早く落ちるので、短時間予測の誤差を見積もる際は保守側の評価になります。(b) では $f_d\tau$ を 3 まで伸ばしており、相関がゼロを何度も横切って負値をとりながら、包絡線 $\pm\sqrt{2/(\pi z)}$（灰の点線）に沿ってゆっくり減衰していく様子と、漸近形（緑の破線）が $f_d\tau \gtrsim 0.4$ でほぼ完全に重なることが確認できます。

自己相関関数が手に入りました。次は、この関数から「チャネルはどれくらいの時間、同じと見なせるか」という数字を切り出しましょう。

コヒーレンス時間の定義

定義が一意でないことを先に認める

コヒーレンス時間 $T_c$ は「チャネルがほぼ変化しないと見なせる時間幅」です。しかし $R_g(\tau)$ は連続的に減衰する関数なので、「ほぼ変化しない」の境界をどこに引くかは本質的に恣意的です。そのため文献ごとに定義が違い、定数倍の食い違いが生じます。ここでは代表的な 3 つを整理し、それぞれがどこから来るのかを明示します。

定義 1: 複素利得の相関がしきい値を下回るまで

もっとも素直な定義は、$R_g(\tau) = J_0(2\pi f_d\tau)$ が 0.5 を下回る時刻です。$J_0(z) = 0.5$ の解は $z \approx 1.5211$ なので、

$$ T_c^{(1)} = \frac{1.5211}{2\pi f_d} \approx \frac{0.242}{f_d} $$

ただしこの定義には落とし穴があります。$g(t)$ の位相まで含めた相関を見ているので、受信機が搬送波位相を追尾している場合の実感とはずれます。実務でよく使われるのは、次の「包絡線・電力の相関」です。

定義 2: 電力（包絡線 2 乗）の相関がしきい値を下回るまで

$g(t)$ が複素ガウス過程のとき、その電力 $|g(t)|^2$ の自己相関には次の厳密な関係があります。平均ゼロの複素ガウス過程では、4 次モーメントが 2 次モーメントで書けるので（複素ガウスモーメント定理）、

$$ \mathbb{E}\!\left[|g(t)|^2|g(t+\tau)|^2\right] = \mathbb{E}[|g|^2]^2 + \left|R_g(\tau)\right|^2 $$

が成り立ちます。したがって、平均を引いた電力の正規化相関係数は

$$ \begin{equation} \rho_{|g|^2}(\tau) = \frac{\mathbb{E}\big[(|g(t)|^2-\bar{P})(|g(t+\tau)|^2-\bar{P})\big]}{\mathrm{Var}[|g|^2]} = \left|\frac{R_g(\tau)}{R_g(0)}\right|^2 = J_0^2(2\pi f_d\tau) \end{equation} $$

となります。$J_0$ が 2 乗されるのがポイントです。包絡線 $|g|$ 自体の相関も、厳密には超幾何関数を含みますが、実用上は $J_0^2$ で近似して差し支えありません。

そこで $J_0^2(2\pi f_d T_c) = 0.5$ を解きます。$J_0(z) = 1/\sqrt{2} \approx 0.7071$ の解は $z \approx 1.1264$ なので、

$$ T_c^{(2)} = \frac{1.1264}{2\pi f_d} \approx \frac{0.1793}{f_d} $$

ここで、$9/(16\pi) = 0.17905$ という値を計算してみてください。上で得た $0.1793$ とほぼ一致します。これが教科書でよく見る

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

の正体です。$9/(16\pi)$ という一見謎めいた定数は、「レイリーフェージングの電力相関が 0.5 に落ちる時刻」を簡潔な分数で書いたものにすぎません。誤差はわずか 0.1 % ほどです。

定義 3: 幾何平均（Rappaport の目安）

定義 2 の $0.179/f_d$ は「かなり厳しめ」の見積もりです。一方、もっとも粗い目安として「フェージングの典型的な時間スケールは $1/f_d$」という言い方もあります。この 2 つは 5 倍以上違うので、実務ではその中間として幾何平均をとった値が広く使われます。

$$ T_c = \sqrt{\frac{9}{16\pi f_d}\cdot\frac{1}{f_d}} = \frac{1}{f_d}\sqrt{\frac{9}{16\pi}} = \frac{3}{4\sqrt{\pi}\,f_d} $$

数値を入れると、

$$ \begin{equation} T_c \approx \frac{0.423}{f_d} \end{equation} $$

これが記事冒頭の表で使った定義で、Rappaport の教科書をはじめ多くの文献で「コヒーレンス時間の経験則」として引用されています。要するに、$T_c$ は $1/f_d$ の 0.18 倍から 1 倍のどこかで、実用的な代表値が 0.42 倍ということです。桁の議論をしているうちは、どの定義を使っても結論は変わりません。

ナイキスト間隔との美しい一致

$T_c \approx 0.423/f_d$ という数字には、もう 1 つ深い意味があります。チャネル利得 $g(t)$ は帯域が $[-f_d, f_d]$ に完全に制限された（帯域制限された）確率過程です。標本化定理より、$g(t)$ を欠落なく再構成するには

$$ T_{\text{Nyquist}} = \frac{1}{2f_d} = \frac{0.5}{f_d} $$

以下の間隔でサンプリングしなければなりません。ここで $T_c \approx 0.423/f_d$ と比べると、

$$ \frac{T_c}{T_{\text{Nyquist}}} = 2 \times 0.423 = 0.846 $$

つまり コヒーレンス時間はナイキスト間隔の約 85 % です。「チャネルを $T_c$ ごとに測る」という設計は、実は「ドップラー過程をぎりぎりナイキスト率で標本化する」ことにほぼ等しいのです。$T_c$ を超える間隔でチャネルを測ると、どんなに賢い補間器を使っても情報が原理的に失われる（エイリアシングする）ことになります。後の Python 実験で、パイロット間隔が $T_c$ を超えた瞬間に推定誤差が跳ね上がるのは偶然ではなく、この標本化定理の帰結です。

同じ 1 本の相関関数の上に、ここまでに出てきた 4 つの「境界」が並んでいます。もっとも厳しい定義 2（電力相関 0.5、$T_c f_d = 0.1793$）から、もっとも緩いナイキスト間隔（$T_c f_d = 0.5$）まで、その差はわずか 2.79 倍しかありません。定義 3 の幾何平均 $0.4231$ がちょうどナイキスト間隔の少し手前に位置していることも、この図で一目でわかります。定義の違いは係数を 2〜3 倍動かすだけで桁は変えない、という感覚を持てば、文献ごとの数字の食い違いに惑わされずに済みます。

コヒーレンス時間の意味が固まったので、次は同じ構造をもつ「遅延スプレッドとコヒーレンス帯域幅」と並べて、フェージングの全体像を整理しましょう。

遅延スプレッドとの双対性、そしてフェージングの 4 分類

2 つの独立した「広がり」

無線チャネルには 2 種類の広がりがあり、それぞれが別のフーリエ対を作ります。

この表の対称性は覚える価値があります。遅延の広がりは周波数方向の相関を壊し、ドップラーの広がりは時間方向の相関を壊す。どちらも「一方の領域での広がりが、他方の領域での相関距離の逆数を与える」という不確定性関係そのものです。

そして重要なのは、この 2 つが独立だということです。遅延スプレッドは散乱体の空間配置（環境の広さ）で決まり、ドップラースプレッドは移動速度と搬送波周波数で決まります。狭い部屋の中を速く動けば「フラットかつ高速フェージング」、広い山間部で止まっていれば「周波数選択性かつ低速フェージング」になります。

4 分類

送信信号の帯域幅 $B_s$（＝シンボルレート $1/T_s$ 程度）と比較して、チャネルは 4 通りに分類されます。

1. フラット・低速フェージング ($B_s \ll B_c$, $T_s \ll T_c$): もっとも扱いやすい。1 シンボルの間にチャネルは動かず、等化器も不要
2. フラット・高速フェージング ($B_s \ll B_c$, $T_s \gtrsim T_c$): シンボル内でチャネルが動く。低ビットレートかつ高速移動の場合に起こり、位相追尾が困難になる
3. 周波数選択性・低速フェージング ($B_s \gtrsim B_c$, $T_s \ll T_c$): 広帯域伝送の典型。等化器や OFDM が必要だが、チャネルはゆっくり変わるので推定しやすい
4. 周波数選択性・高速フェージング ($B_s \gtrsim B_c$, $T_s \gtrsim T_c$): 最悪ケース。高速鉄道での広帯域通信などが該当し、ドップラー領域まで含めた 2 次元の推定が必要になる

実務でよく使われる言い方に注意しましょう。「高速フェージング（fast fading）」という語は、シンボル長と比べて速いという意味であって、絶対的な速さではありません。同じ時速 100 km でも、シンボル長 66 μs の LTE では低速フェージング（$T_c/T_s \approx 34$）、シンボル長 10 ms の低速テレメトリでは高速フェージングになります。

正規化ドップラー $f_d T_s$ という指標

そこで、速さを表す無次元量として正規化ドップラー周波数 $f_d T_s$（$T_s$ はシンボル長、または OFDM シンボル長）がよく使われます。$T_c/T_s \approx 0.423/(f_d T_s)$ なので、$f_d T_s$ が小さいほどチャネルは「1 シンボルあたりゆっくり」です。

3 行目に注目してください。ミリ波では $f_d$ が桁違いに大きいのですが、サブキャリア間隔を 8 倍に広げて OFDM シンボルを短くすることで、$f_d T_s$ は LTE の低速ケースと同程度に抑えられています。5G NR が高い周波数帯で大きなサブキャリア間隔（ニューメロロジー）を用意している理由の 1 つが、まさにこのドップラー耐性です。数式が設計思想に直結している好例と言えます。

横軸が周波数選択性の強さ、縦軸が時間選択性の強さで、両者が独立だからこそ 4 つの象限すべてが存在しえます。LTE の 2 つの例（60 km/h と 300 km/h）はどちらも第 3 象限＝周波数選択性・低速フェージングにあり、OFDM で周波数選択性に対処すればチャネル推定自体は素直にできる領域だとわかります。一方、音声・低ビットレートで時速 300 km というケースは左上の第 2 象限に飛び出しており、帯域は狭いのにシンボル内でチャネルが動くという、等化器では救えないタイプの困難に直面します。

チャネル変動の「速さ」の測り方が整理できました。次は、変動の「深さと頻度」を表すもう 1 組の指標を見ておきましょう。

レベルクロス率と平均フェード時間

$T_c$ は相関関数から定義される平均的な尺度ですが、「1 秒間に何回、どのくらいの長さフェードに落ちるか」を直接知りたい場面もあります。誤り訂正符号のインターリーブ長やフレーム長を決めるときに必要な情報です。

包絡線 $r(t) = |g(t)|$ が、しきい値 $R$ を下から上に横切る単位時間あたりの回数をレベルクロス率 $N_R$ と呼びます。Rice の一般公式に Clarke スペクトルを代入すると、$\rho = R/R_{\text{rms}}$（実効値で正規化したしきい値）を用いて次の形になります。

$$ \begin{equation} N_R = \sqrt{2\pi}\,f_d\,\rho\,e^{-\rho^2} \end{equation} $$

この式の構造は読み解けます。$e^{-\rho^2}$ はレイリー分布の裾（しきい値が高いほど超える確率が下がる）、$\rho$ の因子は横切りの傾き、そして全体にかかる $f_d$ は「時間スケールの逆数」です。レベルクロス率は $f_d$ に正比例する — 速く動けば動くほどフェードの回数が増える、という当然の結果が明示されています。$N_R$ が最大になるのは $\rho = 1/\sqrt{2}$（$-1.5$ dB）のときで、$N_R^{\max} \approx 1.075 f_d$ です。フェードは 1 秒間におよそ $f_d$ 回起きると覚えておけば十分でしょう。

しきい値を下回っている時間の平均を平均フェード時間 $\bar{\tau}$ と呼びます。$r < R$ となる確率（レイリー分布では $1 - e^{-\rho^2}$）を $N_R$ で割れば求まります。

$$ \begin{equation} \bar{\tau} = \frac{1 – e^{-\rho^2}}{N_R} = \frac{e^{\rho^2} – 1}{\rho f_d\sqrt{2\pi}} \end{equation} $$

具体的に見てみましょう。$f_d = 100$ Hz、しきい値を平均電力より 10 dB 低い点（$\rho = 0.316$）に取ると、$N_R \approx 71.7$ 回/秒、$\bar{\tau} \approx 1.33$ ms です。1 秒あたり 70 回以上も 10 dB フェードに落ち、その 1 回あたりの長さが 1.3 ms。LTE のサブフレーム長が 1 ms であることを思い出すと、フェードは 1 サブフレームを丸ごと飲み込む長さであり、サブフレーム内だけで誤り訂正を完結させようとすると全滅する危険があるとわかります。だからこそインターリーブと HARQ が必要になるのです。

(a) の山の頂点が $\rho = 1/\sqrt{2}$（$-1.5$ dB）での $N_R^{\max} = 1.075 f_d$ で、しきい値を上げても下げても回数は減ります。緑の四角が $-10$ dB のしきい値で、$N_R = 0.717 f_d$（$f_d = 100$ Hz なら 71.7 回/秒）と読み取れます。(b) の平均フェード時間は単調増加で、$-10$ dB では $\bar\tau = 0.1327/f_d$（1.33 ms）です。深いフェードほど「回数は減るが 1 回は短い」という関係が、2 枚のグラフの形の違いとして現れています。

理論の道具立てはこれで揃いました。ここからは Python で実際にフェージング波形を作り、導いた式が本当に成り立つのかを一つずつ確かめていきます。

Python 実装（1）: 等方散乱からフェージング波形を作る

導出で使った仮定をそのままコードにします。$M$ 個の散乱体が一様ランダムな到来角 $\alpha_m \sim U[0, 2\pi)$ とランダム位相 $\phi_m \sim U[0,2\pi)$ をもつとして、複素チャネル利得を素朴に足し合わせます。

$$ g(t) = \frac{1}{\sqrt{M}}\sum_{m=1}^{M} e^{j\left(2\pi f_d t\cos\alpha_m + \phi_m\right)} $$

$M \to \infty$ で複素ガウス過程に収束し、期待値をとれば $\mathbb{E}[g^*(t)g(t+\tau)] = \mathbb{E}_\alpha[e^{j2\pi f_d\tau\cos\alpha}] = J_0(2\pi f_d\tau)$ が厳密に出ます。物理モデルとそのまま対応しているので、教育的にはこの素朴な形がいちばんわかりやすいでしょう。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

C_LIGHT = 3.0e8

def max_doppler(v_kmh, fc_hz):
    """最大ドップラー周波数 f_d = v * fc / c [Hz]"""
    return (v_kmh / 3.6) * fc_hz / C_LIGHT

def clarke_fading(fd, fs, n_samples, n_scatter=1024, seed=0, chunk=128):
    """等方散乱モデルで複素フェージング利得 g(t) を生成する"""
    rng = np.random.default_rng(seed)
    t = np.arange(n_samples) / fs
    alpha = rng.uniform(0, 2 * np.pi, n_scatter)   # 到来角（一様分布）
    phi = rng.uniform(0, 2 * np.pi, n_scatter)     # 初期位相（一様分布）
    g = np.zeros(n_samples, dtype=complex)
    # メモリ節約のため散乱体をチャンクに分けて加算する
    for i in range(0, n_scatter, chunk):
        a, p = alpha[i:i + chunk], phi[i:i + chunk]
        g += np.exp(1j * (2 * np.pi * fd * np.cos(a)[:, None] * t[None, :]
                          + p[:, None])).sum(axis=0)
    return t, g / np.sqrt(n_scatter)

fc = 2.0e9
for v in [10, 60, 120]:
    print(f"v = {v:3d} km/h -> f_d = {max_doppler(v, fc):7.2f} Hz, "
          f"T_c = {0.423 / max_doppler(v, fc) * 1e3:6.2f} ms")

実行すると、時速 10 km で $f_d = 18.52$ Hz（$T_c = 22.84$ ms）、時速 60 km で $f_d = 111.11$ Hz（$T_c = 3.81$ ms）、時速 120 km で $f_d = 222.22$ Hz（$T_c = 1.90$ ms）と表示されます。速度を 2 倍にすればコヒーレンス時間はきっちり半分になる、という比例関係が数値でも確認できます。搬送波が 2 GHz のとき、$f_d$ の値がおおよそ「時速の 1.85 倍」になっているのも覚えやすい目安です。

次に、生成した波形の包絡線を dB 表示で描いてみます。速度の違いが「揺れの細かさ」としてどう見えるかを確認しましょう。

fs = 4000.0          # サンプリング周波数 [Hz]
n = 4096             # 約 1.02 秒ぶん
fig, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True)

for ax, v in zip(axes, [10, 60, 120]):
    fd = max_doppler(v, fc)
    t, g = clarke_fading(fd, fs, n, n_scatter=1024, seed=3)
    env_db = 20 * np.log10(np.abs(g))
    ax.plot(t * 1e3, env_db, lw=0.9, color="tab:blue")
    ax.axhline(-10, color="tab:red", ls="--", lw=1.0, label="平均より10dB低い基準")
    ax.set_ylim(-40, 12)
    ax.set_ylabel("受信レベル [dB]")
    ax.set_title(f"時速 {v} km/h（$f_d$ = {fd:.1f} Hz, $T_c$ ≈ {0.423/fd*1e3:.2f} ms）")
    ax.legend(loc="lower right", fontsize=9)
    ax.grid(alpha=0.3)

axes[-1].set_xlabel("時間 [ms]")
fig.suptitle("移動速度によるレイリーフェージング波形の違い（搬送波 2 GHz）")
fig.tight_layout()
plt.show()

3 つの波形を見比べると、縦方向の落ち込みの深さはどれも同じなのに、横方向の細かさだけが速度に比例して変わることがはっきりわかります。時速 10 km では 1 秒間に十数回しか深いフェードが来ませんが、時速 120 km では 200 回近く訪れます。振幅の統計（レイリー分布）は速度に依存せず、変化する時間スケールだけが $f_d$ で決まる — これがドップラースプレッドの本質です。赤い破線（$-10$ dB）を下回る区間の長さが、前節で計算した平均フェード時間 $\bar\tau$ に対応します。

波形が用意できたので、次はその周波数中身を測って、導出したバスタブ形と照合します。

Python 実装（2）: ドップラースペクトルを測ってバスタブ形と比べる

Welch 法で複素信号の両側パワースペクトル密度を推定し、理論式 $S(f) = 1/(\pi f_d\sqrt{1-(f/f_d)^2})$ と重ねます。散乱体数が少ないとスペクトルが線スペクトルの集まりになってギザギザするので、散乱体を多め（4096 個）にし、複数試行で平均します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

fd = 100.0
fs = 1000.0
n = 2 ** 16
n_trial = 20

psd_acc = None
for s in range(n_trial):
    _, g = clarke_fading(fd, fs, n, n_scatter=4096, seed=100 + s)
    f_ax, psd = welch(g, fs=fs, nperseg=1024, return_onesided=False, detrend=False)
    psd_acc = psd if psd_acc is None else psd_acc + psd

f_ax = np.fft.fftshift(f_ax)
psd_est = np.fft.fftshift(psd_acc / n_trial)

# 理論値（端の特異点は評価しない）
f_th = np.linspace(-0.999 * fd, 0.999 * fd, 2000)
s_th = 1.0 / (np.pi * fd * np.sqrt(1 - (f_th / fd) ** 2))

plt.figure(figsize=(10, 5.5))
plt.plot(f_ax, psd_est, lw=1.0, color="tab:blue", label="シミュレーション（Welch推定）")
plt.plot(f_th, s_th, lw=2.0, color="tab:red", label=r"理論 $1/(\pi f_d\sqrt{1-(f/f_d)^2})$")
plt.axvline(-fd, color="gray", ls=":", lw=1.2)
plt.axvline(fd, color="gray", ls=":", lw=1.2, label="$\\pm f_d$（スペクトルの端）")
plt.xlim(-1.6 * fd, 1.6 * fd)
plt.ylim(0, 0.02)
plt.xlabel("周波数オフセット [Hz]")
plt.ylabel("パワースペクトル密度")
plt.title("ドップラースペクトル（バスタブ形）の理論と実測")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# 定量チェック
mask = np.abs(f_ax) < 0.95 * fd
ratio = psd_est[mask] / (1 / (np.pi * fd * np.sqrt(1 - (f_ax[mask] / fd) ** 2)))
print(f"帯域内の実測/理論 比の平均: {ratio.mean():.4f}")
trapz = getattr(np, "trapezoid", np.trapz)   # NumPy 2.x では trapezoid
rms = np.sqrt(trapz(f_ax ** 2 * psd_est, f_ax) / trapz(psd_est, f_ax))
print(f"RMSドップラースプレッド: 実測 {rms:.2f} Hz / 理論 {fd/np.sqrt(2):.2f} Hz")

推定スペクトルは、$f = 0$ で最小、$\pm f_d$ に近づくにつれて跳ね上がり、$|f| > f_d$ ではストンとゼロに落ちるという、まさにバスタブの断面図になります。帯域内の実測／理論比の平均は 1.003 とほぼ 1 で、RMS ドップラースプレッドは実測 70.9 Hz に対し理論値 $f_d/\sqrt{2} = 70.7$ Hz と 0.3 % 以内で一致しました。「一様な到来角＋ヤコビアン $f_d|\sin\theta|$」という導出が正しかったことが数値で裏づけられたわけです。

もう 1 つ注目したいのは、$|f| > f_d$ の外側にまったく電力が漏れていないことです。残っているのは数値計算の丸め誤差レベルの成分だけで、グラフでは完全にゼロに見えます。チャネル利得が厳密に帯域制限された過程であるという先ほどの主張が、ここで確認できます。この事実が、後のパイロット間隔の議論で標本化定理を持ち出す根拠になります。

スペクトルが確認できたので、そのフーリエ変換である自己相関がベッセル関数になっているかを直接測ってみましょう。

Python 実装（3）: 自己相関と $J_0$ の一致を確かめる

FFT を使って標本自己相関を計算し、$J_0(2\pi f_d\tau)$ と重ねます。$\tau$ 軸は $f_d\tau$（正規化遅れ）で表すと、速度によらず同じ曲線に乗るはずです。

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

def autocorr(x, n_lag):
    """FFT による標本自己相関（非バイアス推定）"""
    n = len(x)
    X = np.fft.fft(x - x.mean(), 2 * n)
    r = np.fft.ifft(np.abs(X) ** 2)[:n_lag] / np.arange(n, n - n_lag, -1)
    return r

fd = 100.0
fs = 2000.0
n = 2 ** 18
n_lag = int(3.0 / fd * fs)

acc = np.zeros(n_lag, dtype=complex)
for s in range(8):
    _, g = clarke_fading(fd, fs, n, n_scatter=2048, seed=200 + s)
    acc += autocorr(g, n_lag)
r_g = (acc / 8).real
r_g /= r_g[0]

lag = np.arange(n_lag) / fs
plt.figure(figsize=(10, 5.5))
plt.plot(fd * lag, r_g, lw=1.4, color="tab:blue", label="シミュレーションの自己相関")
plt.plot(fd * lag, j0(2 * np.pi * fd * lag), "--", lw=2.0, color="tab:red",
         label=r"理論 $J_0(2\pi f_d\tau)$")
plt.axhline(0, color="gray", lw=0.8)
plt.axvline(0.3827, color="tab:green", ls=":", lw=1.5, label="最初のゼロ点 $f_d\\tau=0.383$")
plt.xlabel(r"正規化遅れ $f_d\tau$")
plt.ylabel("正規化自己相関")
plt.title("時間自己相関はベッセル関数 $J_0$ になる")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

err = np.max(np.abs(r_g[:int(2 / fd * fs)] - j0(2 * np.pi * fd * lag[:int(2 / fd * fs)])))
print(f"f_d*tau < 2 の範囲での最大誤差: {err:.4f}")

2 本の曲線はほぼ完全に重なり、$f_d\tau < 2$ の範囲での最大誤差は 0.0064 と、相関値のスケール（最大 1）から見て 1 % 未満に収まります。特に、$f_d\tau = 0.383$ で相関がゼロを横切り、その後に負の値をとってから振動しながら減衰していく — という $J_0$ 特有の挙動が、シミュレーション波形からもきちんと再現されました。相関が単調減少ではなく振動する点は、指数相関を仮定した簡易モデル（AR(1) など）との決定的な違いです。

負の相関が現れる物理的な意味も押さえておきましょう。$\tau \approx 0.5/f_d$ 付近では、前方波と後方波の位相差がちょうど $\pi$ 程度ずれ、合成ベクトルが「反転しやすい」状態になります。2 波モデルの $\cos(2\pi f_d t)$ が負になるのと同じ話です。

自己相関が確認できたので、次はここからコヒーレンス時間を実測し、$9/(16\pi f_d)$ という式を検証します。

Python 実装（4）: 速度を振ってコヒーレンス時間を実測する

定義 2 に従い、電力 $|g(t)|^2$ の正規化自己相関が 0.5 を下回る時刻を $T_c$ と定義して測ります。理論では $J_0^2(2\pi f_d T_c) = 0.5$、すなわち $T_c f_d \approx 0.179$ になるはずです。

import numpy as np
import matplotlib.pyplot as plt

def measure_Tc(fd, n_sample=2 ** 17, n_trial=8, seed0=1000):
    """電力相関が 0.5 を下回る時刻としてコヒーレンス時間を実測する"""
    fs = 40 * fd                       # 十分なオーバーサンプリング
    n_lag = int(1.0 / fd * fs)
    acc = np.zeros(n_lag)
    for s in range(n_trial):
        _, g = clarke_fading(fd, fs, n_sample, n_scatter=2048, seed=seed0 + s)
        acc += autocorr(np.abs(g) ** 2, n_lag).real
    r = acc / acc[0]
    lag = np.arange(n_lag) / fs
    k = int(np.argmax(r < 0.5))        # 初めて 0.5 を下回るインデックス
    return np.interp(0.5, [r[k], r[k - 1]], [lag[k], lag[k - 1]])

fc = 2.0e9
speeds = np.array([10, 30, 60, 120, 240])
fds = np.array([max_doppler(v, fc) for v in speeds])
Tc_meas = np.array([measure_Tc(fd, seed0=1000 + int(v)) for fd, v in zip(fds, speeds)])

for v, fd, tc in zip(speeds, fds, Tc_meas):
    print(f"v={v:4d} km/h  f_d={fd:7.2f} Hz  T_c実測={tc*1e3:7.3f} ms  "
          f"9/(16pi f_d)={0.17905/fd*1e3:7.3f} ms  T_c*f_d={tc*fd:.4f}")

出力は次のようになります。時速 10 km で $T_c f_d = 0.181$、30 km/h で 0.179、60 km/h で 0.180、120 km/h で 0.181、240 km/h で 0.179 — 速度を 24 倍に振っても $T_c f_d$ は 0.179 前後で一定です。理論値 $9/(16\pi) = 0.1790$ との差は最大でも 1.3 %（時速 120 km の 0.1813）に収まっており、「コヒーレンス時間は $f_d$ に反比例する」という主張と、その比例定数の両方が確認できました。

この結果をグラフでも見ておきましょう。

plt.figure(figsize=(10, 5.5))
v_fine = np.linspace(5, 260, 200)
fd_fine = np.array([max_doppler(v, fc) for v in v_fine])
plt.plot(v_fine, 0.17905 / fd_fine * 1e3, "-", lw=2.0, color="tab:red",
         label=r"理論 $T_c = 9/(16\pi f_d)$")
plt.plot(v_fine, 0.423 / fd_fine * 1e3, "--", lw=2.0, color="tab:orange",
         label=r"経験則 $T_c = 0.423/f_d$（幾何平均）")
plt.plot(v_fine, 0.5 / fd_fine * 1e3, ":", lw=2.0, color="tab:green",
         label=r"ナイキスト間隔 $1/(2f_d)$")
plt.plot(speeds, Tc_meas * 1e3, "o", ms=9, color="tab:blue",
         label="シミュレーション実測（電力相関0.5）")
plt.yscale("log")
plt.xlabel("移動速度 [km/h]（搬送波 2 GHz）")
plt.ylabel("コヒーレンス時間 [ms]")
plt.title("コヒーレンス時間は移動速度に反比例する")
plt.legend()
plt.grid(alpha=0.3, which="both")
plt.tight_layout()
plt.show()

両対数に近い形（縦軸のみ対数）で描くと、実測点が理論直線にきれいに乗ります。3 本の線が上下に並んでいるのは、前節で述べた「定義の恣意性」の可視化です。もっとも厳しい $9/(16\pi f_d)$ と、もっとも緩いナイキスト間隔 $1/(2f_d)$ の間はわずか 2.8 倍しかなく、幾何平均の $0.423/f_d$ はその中間、ナイキスト間隔の 85 % に位置しています。定義の違いは桁を変えない、という感覚を持っておけば実務では困りません。

ここまでは統計量の検証でした。最後に、コヒーレンス時間が実際の受信機性能にどう効くのかを、パイロット間隔とチャネル推定誤差の関係で示します。

Python 実装（5）: パイロット間隔がコヒーレンス時間を超えると何が起こるか

受信機は $D$ シンボルおきに挿入されたパイロットでチャネルを測り、その間のデータシンボルのチャネルを補間で推定します。ここでは 2 つの方式を比べます。

• 前値保持（ZOH）: 直近のパイロットで測った値をそのまま使う
• 線形補間: 前後 2 つのパイロットを直線で結ぶ

評価指標は正規化平均二乗誤差 $\mathrm{NMSE} = \mathbb{E}|\hat{g}-g|^2 / \mathbb{E}|g|^2$ です。パイロット自体にも受信雑音が乗るので、SNR = 20 dB の白色雑音を加えます。

理論的な予測も立てておきましょう。雑音がなければ、前値保持で遅れ $\tau$ 経過したときの誤差は

$$ \mathbb{E}\left|g(t+\tau)-g(t)\right|^2 = 2\big(R_g(0) – \mathrm{Re}\,R_g(\tau)\big) = 2\big(1 – J_0(2\pi f_d\tau)\big) $$

です。$\tau$ が小さいときは $J_0 \approx 1 – (\pi f_d\tau)^2$ より $\mathrm{NMSE} \approx 2(\pi f_d\tau)^2$ と $\tau^2$ で増え、$\tau$ が大きくなると $J_0 \to 0$ で NMSE は 2 に飽和します。NMSE が 1 を超えるということは、「推定するくらいならゼロと答えたほうがマシ」という壊滅的な状態です。

import numpy as np
import matplotlib.pyplot as plt

fc, Ts, snr_db = 2.0e9, 71.4e-6, 20.0   # LTE 相当の OFDM シンボル長
sigma2 = 10 ** (-snr_db / 10)
rng = np.random.default_rng(7)

def pilot_nmse(v_kmh, D, n_sym=2 ** 16, n_trial=6):
    """パイロット間隔 D シンボルのときのチャネル推定 NMSE（ZOH / 線形補間）"""
    fd = max_doppler(v_kmh, fc)
    e_zoh, e_lin = [], []
    for s in range(n_trial):
        _, g = clarke_fading(fd, 1 / Ts, n_sym, n_scatter=1024, seed=int(v_kmh * 100 + D + s))
        idx = np.arange(0, n_sym, D)                       # パイロット位置
        noise = np.sqrt(sigma2 / 2) * (rng.standard_normal(len(idx))
                                       + 1j * rng.standard_normal(len(idx)))
        est = g[idx] + noise                               # パイロットでの観測値
        n = np.arange(n_sym)
        zoh = est[np.minimum(n // D, len(idx) - 1)]        # 前値保持
        lin = np.interp(n, idx, est.real) + 1j * np.interp(n, idx, est.imag)
        sl = slice(0, (len(idx) - 1) * D)                  # 外挿区間は除く
        p = np.mean(np.abs(g[sl]) ** 2)
        e_zoh.append(np.mean(np.abs(zoh[sl] - g[sl]) ** 2) / p)
        e_lin.append(np.mean(np.abs(lin[sl] - g[sl]) ** 2) / p)
    return fd, np.mean(e_zoh), np.mean(e_lin)

v = 60
Ds = [1, 2, 4, 7, 14, 28, 56, 112]
res = [pilot_nmse(v, D) for D in Ds]
fd = res[0][0]
Tc = 0.423 / fd
print(f"v={v} km/h, f_d={fd:.1f} Hz, T_c={Tc*1e3:.3f} ms = {Tc/Ts:.1f} シンボル")
for D, (_, eh, el) in zip(Ds, res):
    print(f" D={D:4d}  tau/T_c={D*Ts/Tc:6.3f}  ZOH={10*np.log10(eh):7.2f} dB  "
          f"線形補間={10*np.log10(el):7.2f} dB")

出力を見ると、時速 60 km・2 GHz では $T_c = 3.81$ ms、OFDM シンボル 53 個ぶんに相当します。パイロット間隔 $D$ を増やしていくと、NMSE は次のように推移します。

読み取れることは 3 つあります。第 1 に、線形補間は $\tau/T_c \lesssim 0.3$ の範囲では雑音制限（$-21$ dB 前後）のまま横ばいで、パイロットを増やしても誤差は減りません。むしろ $D = 1$（全シンボルがパイロット）のほうが $-19.9$ dB とわずかに悪いくらいです。これは、$D \ge 2$ では前後 2 つのパイロットを平均する形になり雑音が $3$ dB 近く抑えられる一方、$D = 1$ では平均する相手がいないからです。この領域ではパイロットのオーバーヘッドが純粋に無駄なので、間隔を広げてスループットを稼ぐべきです。第 2 に、$\tau$ が $T_c$ に近づくと誤差が急激に立ち上がります。$0.53 T_c$ で $-17.6$ dB、$1.05 T_c$ で $-8.3$ dB、$2.1 T_c$ でついに $0$ dB — $T_c$ を超えた途端に推定が使い物にならなくなるという直感が、数値で裏づけられました。第 3 に、前値保持は線形補間よりはるかに脆く、$0.26 T_c$ の時点ですでに $-10.9$ dB まで劣化しています。補間の賢さが実効的な許容間隔を 3〜4 倍広げているわけです。

$T_c$ を超えると壊滅するのは、単に「相関が下がるから」だけではありません。前に述べたとおり $T_c \approx 0.85 \times 1/(2f_d)$ なので、$\tau > T_c$ はドップラー過程をナイキスト率未満で標本化することにほぼ等しく、どんな補間器を使っても原理的に復元できない領域に入るのです。この点を押さえておくと、「もっと高性能な補間フィルタを作れば解決する」という誤った期待を持たずに済みます。

理論値との突き合わせもしておきましょう。

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

taus = np.array([D * Ts for D in Ds])
nmse_zoh = np.array([r[1] for r in res])
nmse_lin = np.array([r[2] for r in res])
# 前値保持の理論値：遅れ 0..D-1 シンボルの平均 + 雑音項
theory = np.array([np.mean([2 - 2 * j0(2 * np.pi * fd * k * Ts) for k in range(D)]) + sigma2
                   for D in Ds])

plt.figure(figsize=(10, 5.5))
plt.semilogx(taus / Tc, 10 * np.log10(nmse_zoh), "o-", label="前値保持（実測）")
plt.semilogx(taus / Tc, 10 * np.log10(theory), "k--", lw=1.2, label="前値保持（理論）")
plt.semilogx(taus / Tc, 10 * np.log10(nmse_lin), "s-", label="線形補間（実測）")
plt.axvline(1.0, color="tab:red", lw=1.8, ls="--", label="$\\tau = T_c$")
plt.axhline(0.0, color="gray", lw=1.0, ls=":")
plt.xlabel(r"パイロット間隔 / コヒーレンス時間 $\tau/T_c$")
plt.ylabel("チャネル推定 NMSE [dB]")
plt.title("パイロット間隔が $T_c$ を超えるとチャネル推定は破綻する（SNR 20 dB）")
plt.legend()
plt.grid(alpha=0.3, which="both")
plt.tight_layout()
plt.show()

前値保持の実測点は理論曲線 $2(1-J_0(2\pi f_d\tau))$ の平均に 0.2 dB 以内で乗っており、導いた自己相関の式がそのまま推定誤差の予測式として使えることがわかります。$\tau = T_c$ の縦線を境に、線形補間の曲線が水平から急な立ち上がりへ転じる「膝」が見えます。実際の規格でパイロット間隔を $T_c$ の 1/4 から 1/3 程度に設定するのは、この膝から十分な余裕をとるためです。たとえば LTE のセル固有参照信号は時間方向に約 0.25 ms 間隔で配置されており、$1/(2\times 0.25\ \text{ms}) = 2$ kHz までのドップラーに追随できます。時速 300 km・2.6 GHz の $f_d = 722$ Hz なら十分な余裕があるとわかります。

数値実験がひととおり終わりました。最後に、ここまでの結果が実際の設計判断にどう効くのかを整理しておきます。

設計への含意

OFDM のサブキャリア間干渉

OFDM では、1 シンボル長 $T$ の間チャネルが一定であることを前提にサブキャリアの直交性が保たれます。チャネルが動くと直交性が崩れ、サブキャリア間干渉（ICI）が発生します。その大きさを見積もってみましょう。

時変チャネル $g(t)$ を 1 シンボル内で DFT した係数を $H[\ell]$ とすると、$\ell = 0$ の成分が所望信号、$\ell \neq 0$ の成分が ICI です。全電力が保存されるので $\sum_\ell \mathbb{E}|H[\ell]|^2 = 1$ であり、ICI 電力は

$$ P_{\text{ICI}} = 1 – \mathbb{E}|H[0]|^2 = 1 – \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1} J_0\!\big(2\pi f_d (n-m) T_s’\big) $$

と書けます（$T_s’ = T/N$ は標本間隔）。$f_d T \ll 1$ のとき $J_0(z) \approx 1 – z^2/4$ を代入すると、和は $(\pi f_d T_s’)^2 \sum\sum (n-m)^2$ となります。ここで $\frac{1}{N^2}\sum\sum(n-m)^2 = 2\mathrm{Var}[n] \approx N^2/6$ を使えば、

$$ \begin{equation} P_{\text{ICI}} \approx \frac{(\pi f_d T)^2}{6}, \qquad \mathrm{SIR} \approx \frac{6}{(\pi f_d T)^2} \end{equation} $$

という簡潔な近似式が得られます。数値検証したところ、$f_d T = 0.01$ で厳密値 $1.645\times10^{-4}$ に対し近似値 $1.645\times10^{-4}$、$f_d T = 0.1$ でも厳密 $1.63\times10^{-2}$ に対し近似 $1.64\times10^{-2}$ と、実用範囲では小数点以下 2 桁まで一致します。

具体例で見ましょう。LTE（サブキャリア間隔 15 kHz、$T = 66.7$ μs）で時速 300 km・2.6 GHz なら $f_d T = 0.048$ となり、SIR は約 24 dB です。64QAM を使うには所要 SNR が 20 dB 前後なので、まだぎりぎり成立しますが余裕はありません。時速 500 km・3.5 GHz にすると $f_d T = 0.108$、SIR は 17 dB まで落ち、高次変調は現実的でなくなります。サブキャリア間隔を広げて $T$ を短くするのが唯一の根本対策であり、5G NR が高周波数帯で 60 kHz・120 kHz のニューメロロジーを用意している理由がここにあります。

(a) では厳密値（$J_0$ の二重和）と近似式が両対数上でほとんど区別できないほど重なっており、$f_d T \lesssim 0.2$ の実用範囲では近似で十分だとわかります。(b) の SIR 曲線には、LTE で時速 300 km・2.6 GHz のケース（緑、$f_d T = 0.048$、SIR 24.2 dB）と時速 500 km・3.5 GHz のケース（橙、$f_d T = 0.108$、SIR 17.2 dB）を置きました。後者は 64QAM の所要 SNR の目安である 20 dB（灰の点線）をすでに下回っており、変調次数を落とすかシンボル長 $T$ を短くするしかないことが図から直接読み取れます。

CSI フィードバックとスケジューリング

適応変調・符号化（AMC）やビームフォーミングは、端末が測ったチャネル情報を基地局に返し、基地局がそれに基づいて送信パラメータを決める仕組みです。測定してから実際に送信されるまでの遅延を $\Delta$ とすると、送信時のチャネルと測定時のチャネルの相関は $J_0(2\pi f_d\Delta)$ しかありません。$\Delta \gtrsim T_c$ なら、フィードバック情報はランダムな値と大差なくなります。

LTE/NR の CSI フィードバック遅延は数ミリ秒オーダーです。時速 60 km・2 GHz の $T_c = 3.8$ ms と比べると、すでに際どい水準です。高速移動時に AMC ゲインが消え、周波数選択スケジューリングが機能しなくなるのはこのためで、対策としては (1) 平均的なチャネル統計に基づく保守的な設定に切り替える、(2) チャネル予測でリードタイムを稼ぐ、(3) 開ループ方式（送信ダイバーシチ）に切り替える、といった選択肢がとられます。

時間ダイバーシチとインターリーブ

一方で、チャネルが速く変わることは悪いことばかりではありません。符号化とインターリーブを組み合わせれば、深いフェードに落ちたシンボルの誤りを、別の時刻の良好なシンボルで救えます。これが時間ダイバーシチです。

効果を得るには、インターリーブ長が相関の切れる時間、つまり $\tau_0 \approx 0.38/f_d$ を十分に超える必要があります。時速 5 km の歩行者（$f_d = 9.3$ Hz）だと $\tau_0 \approx 41$ ms なので、意味のある時間ダイバーシチを得るには 100 ms 以上のインターリーブが必要になり、遅延要求の厳しい音声通信では成立しません。低速移動こそが時間ダイバーシチにとっては最悪ケースという、やや逆説的な事実です。だからこそ低速時には周波数ダイバーシチや空間ダイバーシチに頼ることになります。

測定側から見た $f_d$ の推定

ここまでは $f_d$ が既知という前提でしたが、実際の受信機は $f_d$ そのものを推定してパイロット間隔や補間フィルタ帯域を適応的に切り替えます。代表的な手法は 2 つです。

1. レベルクロス率法: 包絡線が実効値を横切る回数を数え、$N_R = \sqrt{2\pi}f_d\rho e^{-\rho^2}$ から逆算する。実装が軽い反面、雑音や見通し成分があるとバイアスが乗る
2. 相関法: パイロットから得た $\hat{R}_g(\tau)$ を短い遅れで測り、$J_0(2\pi f_d\tau) \approx 1-(\pi f_d\tau)^2$ の 2 次係数から $f_d$ を求める。雑音に強いが、遅れの選び方に注意が必要

いずれも、本記事で導いた $J_0$ の形が推定器の設計根拠になっています。理論式が単なる説明用の道具ではなく、そのままアルゴリズムになるという良い例でしょう。

まとめ

本記事では、移動によって生じるドップラースプレッドと、その時間領域の対応物であるコヒーレンス時間を、導出から実測まで通して解説しました。

• 到来角 $\theta$ の波は $f = f_d\cos\theta$ のドップラーシフトを受ける。$f_d = v f_c/c$ は速度と搬送波周波数の積に比例する
• 到来角が一様分布のとき、変数変換のヤコビアン $|df/d\theta| = f_d\sqrt{1-(f/f_d)^2}$ と 2 つの逆像から、バスタブ形のドップラースペクトル $S(f) = 1/(\pi f_d\sqrt{1-(f/f_d)^2})$ が得られる。両端が発散するのは $\cos\theta$ の停留点で角度→周波数の写像が潰れるため
• RMS ドップラースプレッドは $\sigma_D = f_d/\sqrt{2}$
• $S(f)$ の逆フーリエ変換はベッセルの積分表示により $R_g(\tau) = J_0(2\pi f_d\tau)$ になる。最初のゼロ点は $\tau = 0.383/f_d$
• 電力の相関は厳密に $J_0^2(2\pi f_d\tau)$ であり、それが 0.5 になる時刻から $T_c \approx 9/(16\pi f_d) \approx 0.179/f_d$、$1/f_d$ との幾何平均をとると広く使われる目安 $T_c \approx 0.423/f_d$ が出る
• $T_c$ はドップラー過程のナイキスト間隔 $1/(2f_d)$ の約 85 % にあたる。パイロット間隔が $T_c$ を超えるとチャネル推定が原理的に破綻する
• シミュレーションでは、スペクトル・自己相関・$T_c f_d = 0.179$・パイロット間隔と NMSE の関係のすべてが理論と一致した
• 設計面では、OFDM の ICI（$P_{\text{ICI}} \approx (\pi f_d T)^2/6$）、CSI フィードバック遅延、インターリーブ長のいずれもがコヒーレンス時間で律速される

ドップラースプレッドとコヒーレンス時間は、遅延スプレッドとコヒーレンス帯域幅と対をなす、時変チャネルを記述する 2 本柱の一方です。この 2 組 4 つの量を押さえておけば、どんな無線システムの仕様書を読んでも「なぜこの数字なのか」が推測できるようになります。

次のステップとして、以下の記事も参考にしてください。

• Jakes モデルによるフェージングシミュレーション
• マルチパス伝搬 — 反射・回折・散乱がつくる受信波の重ね合わせ
• ライスフェージングチャネル — 見通し成分があるときの統計
• フェージングチャネルのアウテージ確率
• ダイバーシチ技術によるフェージング対策
• OFDM の原理と実装