レーダーで2機の航空機を検出したとき、2機が十分離れていれば個別に識別できますが、接近していると1つの目標として誤認される可能性があります。従来のビームフォーミング（遅延和法）では、アンテナアレイのビーム幅以下の角度差を持つ信号源を分離できません。ビーム幅は $\lambda/D$（$D$ はアレイ開口長）に比例するため、分解能の改善にはアレイを大きくするしかないように思えます。
しかし、信号処理のアプローチによってこの限界を超えることができます。その代表的な手法がCapon法（MVDR法）とMUSIC法です。これらは信号の統計的性質を利用して、従来法の数分の1のビーム幅で信号源の到来方向（DOA: Direction of Arrival）を推定できます。
適応ビームフォーミングの理論を理解することは、以下の場面で大きな力を発揮します。
- レーダーシステム: 密集した複数の目標を分離・追尾するために、MUSIC法やCapon法が広く使われています。特に航空管制レーダーやミサイル誘導レーダーでは、角度分解能が生命線です
- 無線通信（5G/6G）: Massive MIMOにおいて、多数のユーザーの到来方向を高精度に推定し、ユーザーごとにビームを向ける技術の基盤となります
- 音響工学: マイクロホンアレイによる音源方向推定（音源定位）にも同じ理論が適用されます。スマートスピーカーやロボットの音声認識に不可欠な技術です
本記事の内容
- 従来ビームフォーミング（遅延和法）の分解能限界
- アレイ信号モデル（ステアリングベクトルと共分散行列）
- Capon法（MVDR）の導出と性質
- MUSIC法の原理（信号部分空間と雑音部分空間）
- 固有値分解に基づくDOA推定
- Pythonシミュレーション（2信号源の分離、3手法の比較）
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
従来ビームフォーミングの分解能限界
遅延和法（Delay-and-Sum）
$M$ 素子の等間隔リニアアレイ（ULA: Uniform Linear Array）を考えます。素子間隔を $d$、波長を $\lambda$ とします。角度 $\theta$（アレイの法線からの角度）からの平面波が到来するとき、各素子の受信信号には位相差が生じます。
$m$ 番目の素子（$m = 0, 1, \ldots, M-1$）の受信信号は
$$ x_m(t) = s(t) \, e^{-j m k d \sin\theta} + n_m(t) $$
ここで $s(t)$ は信号、$n_m(t)$ は雑音、$k = 2\pi/\lambda$ は波数です。
遅延和ビームフォーミングは、走査方向 $\theta_0$ に対して各素子の信号の位相を補正し、コヒーレントに加算する処理です。
$$ y(\theta_0) = \sum_{m=0}^{M-1} x_m \, e^{j m k d \sin\theta_0} = \bm{w}^H \bm{x} $$
ここで重みベクトル $\bm{w}$ と受信ベクトル $\bm{x}$ は
$$ \bm{w} = \frac{1}{M}\begin{pmatrix} 1 \\ e^{-jkd\sin\theta_0} \\ \vdots \\ e^{-j(M-1)kd\sin\theta_0} \end{pmatrix}, \quad \bm{x} = \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_{M-1} \end{pmatrix} $$
出力パワーは
$$ P_{\text{CBF}}(\theta_0) = E[|y(\theta_0)|^2] = \bm{w}^H \bm{R} \bm{w} $$
ここで $\bm{R} = E[\bm{x}\bm{x}^H]$ は受信信号の共分散行列（covariance matrix）です。
分解能限界
遅延和法のアレイパターンは、$\theta = 0$ にメインビームを向けた場合
$$ \text{AF}(\theta) = \frac{1}{M}\frac{\sin(M\pi d\sin\theta/\lambda)}{\sin(\pi d\sin\theta/\lambda)} $$
この関数のメインビームの半値幅（HPBW）は、$d = \lambda/2$ のとき近似的に
$$ \text{HPBW} \approx \frac{0.886\lambda}{Md} = \frac{1.77}{M} \quad [\text{rad}] $$
レイリー限界（Rayleigh criterion）によれば、2つの信号源を分離するためには、角度差がHPBW以上でなければなりません。$M = 8$, $d = \lambda/2$ のとき HPBW $\approx 12.7°$ であり、これ以下の角度差の信号源は分離できません。
この限界は、アレイの物理的な開口長 $D = (M-1)d$ で決まる回折限界に相当します。しかし、信号の統計的性質を利用すれば、この回折限界を超えた分解能が達成可能です。その最初のアプローチがCapon法です。
アレイ信号モデル
ステアリングベクトル
DOA推定の理論を展開するために、まず信号モデルを厳密に定式化しましょう。
$K$ 個の信号源がそれぞれ方向 $\theta_1, \theta_2, \ldots, \theta_K$ から到来しているとします。方向 $\theta_k$ からの平面波に対するステアリングベクトル（steering vector）は
$$ \bm{a}(\theta_k) = \begin{pmatrix} 1 \\ e^{-jkd\sin\theta_k} \\ e^{-j2kd\sin\theta_k} \\ \vdots \\ e^{-j(M-1)kd\sin\theta_k} \end{pmatrix} $$
です。ステアリングベクトルは、方向 $\theta_k$ からの到来波が各素子でどのような位相関係を持つかを表しています。
受信信号ベクトル
$K$ 個の信号と雑音を含む受信信号ベクトルは
$$ \bm{x}(t) = \sum_{k=1}^{K} \bm{a}(\theta_k) s_k(t) + \bm{n}(t) = \bm{A}\bm{s}(t) + \bm{n}(t) $$
と表されます。ここで
$$ \bm{A} = [\bm{a}(\theta_1), \bm{a}(\theta_2), \ldots, \bm{a}(\theta_K)] $$
は $M \times K$ のステアリング行列、$\bm{s}(t) = [s_1(t), \ldots, s_K(t)]^T$ は信号ベクトル、$\bm{n}(t)$ は雑音ベクトルです。
共分散行列
受信信号の共分散行列は
$$ \bm{R} = E[\bm{x}\bm{x}^H] = \bm{A}\bm{R}_s\bm{A}^H + \sigma_n^2\bm{I} $$
ここで $\bm{R}_s = E[\bm{s}\bm{s}^H]$ は信号の共分散行列（$K \times K$）、$\sigma_n^2$ は雑音電力、$\bm{I}$ は $M \times M$ の単位行列です。雑音は各素子で無相関かつ等電力（白色雑音）と仮定しています。
実際の測定では、共分散行列は $N$ 個のスナップショット（時間サンプル）から標本共分散行列として推定します。
$$ \hat{\bm{R}} = \frac{1}{N}\sum_{n=1}^{N} \bm{x}(t_n)\bm{x}^H(t_n) $$
スナップショット数 $N$ が大きいほど推定精度が向上します。
信号モデルを定式化したところで、Capon法の導出に進みましょう。
Capon法（MVDR法）
発想
Capon法（1969年）は、最小分散無歪み応答（MVDR: Minimum Variance Distortionless Response）とも呼ばれます。
従来ビームフォーミングの問題は、メインビーム以外の方向からの信号（干渉波）もサイドローブを通じて出力に混入し、分解能を低下させることです。Capon法の発想は、注目方向 $\theta_0$ からの信号は歪みなく通過させながら、他の方向からの干渉とノイズの出力パワーを最小化する重みベクトルを求めることです。
これは、料理で特定の味（信号）を保ちつつ、他の雑味（干渉・雑音）を最小化するフィルターに例えられます。
数学的定式化
Capon法の最適化問題は次のように定式化されます。
$$ \min_{\bm{w}} \bm{w}^H \bm{R} \bm{w} \quad \text{subject to} \quad \bm{w}^H \bm{a}(\theta_0) = 1 $$
ここで、目的関数 $\bm{w}^H \bm{R} \bm{w}$ は出力パワー、制約条件 $\bm{w}^H \bm{a}(\theta_0) = 1$ は方向 $\theta_0$ からの信号を歪みなく通過させることを保証します。
ラグランジュの未定乗数法による解
ラグランジアンを構成します。
$$ \mathcal{L}(\bm{w}, \mu) = \bm{w}^H \bm{R} \bm{w} – \mu(\bm{w}^H \bm{a}(\theta_0) – 1) $$
$\bm{w}$ で微分してゼロとおくと
$$ \frac{\partial \mathcal{L}}{\partial \bm{w}^*} = \bm{R}\bm{w} – \mu\bm{a}(\theta_0) = 0 $$
$\bm{R}\bm{w} = \mu\bm{a}(\theta_0)$ を $\bm{w}$ について解くと
$$ \bm{w} = \mu \bm{R}^{-1}\bm{a}(\theta_0) $$
制約条件 $\bm{w}^H \bm{a}(\theta_0) = 1$ に代入して $\mu$ を求めると
$$ \mu^* \bm{a}^H(\theta_0) \bm{R}^{-1} \bm{a}(\theta_0) = 1 $$
$$ \mu = \frac{1}{\bm{a}^H(\theta_0) \bm{R}^{-1} \bm{a}(\theta_0)} $$
したがって、Capon法の最適重みベクトルは
$$ \bm{w}_{\text{Capon}} = \frac{\bm{R}^{-1}\bm{a}(\theta_0)}{\bm{a}^H(\theta_0)\bm{R}^{-1}\bm{a}(\theta_0)} $$
Caponスペクトル
Capon法の出力パワー（空間スペクトル）は
$$ P_{\text{Capon}}(\theta_0) = \bm{w}_{\text{Capon}}^H \bm{R} \bm{w}_{\text{Capon}} = \frac{1}{\bm{a}^H(\theta_0)\bm{R}^{-1}\bm{a}(\theta_0)} $$
この $P_{\text{Capon}}(\theta_0)$ を $\theta_0$ の関数としてプロットすると、信号源の方向にピークが現れます。
Capon法の特性
Capon法は従来法に比べて以下の利点があります。
高い分解能: Capon法のビーム幅は従来法より狭く、接近した信号源を分離できる可能性が高くなります。ただし、分解能はSNR（信号対雑音比）とスナップショット数に依存します。
サイドローブの低減: Capon法は干渉方向にヌル（零点）を自動的に形成するため、サイドローブが効果的に抑制されます。
一方で、以下の制約もあります。
共分散行列の逆行列が必要: $M \times M$ の行列の逆行列を計算するため、計算量が $O(M^3)$ です。また、スナップショット数 $N$ が素子数 $M$ 以上でないと、標本共分散行列がランク落ちして逆行列が存在しません。
コヒーレント信号に弱い: マルチパス環境などで信号同士が相関（コヒーレント）の場合、共分散行列のランクが低下し、分解能が劣化します。
Capon法は従来法より高い分解能を持ちますが、さらに高い分解能を持つ手法がMUSIC法です。次にその原理を見ていきましょう。
MUSIC法 — 部分空間に基づくDOA推定
部分空間の概念
MUSIC法（MUltiple SIgnal Classification, 1986年 R.O. Schmidt）は、共分散行列の固有値分解に基づくDOA推定法です。
MUSIC法の鍵となるアイデアは、共分散行列の固有ベクトルが2つの直交する部分空間に分かれるという事実です。
共分散行列 $\bm{R}$ の固有値分解を考えます。
$$ \bm{R} = \bm{A}\bm{R}_s\bm{A}^H + \sigma_n^2\bm{I} $$
$\bm{R}$ は $M \times M$ のエルミート行列であるため、$M$ 個の実固有値 $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_M$ と対応する直交固有ベクトル $\bm{e}_1, \bm{e}_2, \ldots, \bm{e}_M$ を持ちます。
信号部分空間と雑音部分空間
$K$ 個の信号源がある場合（$K < M$）、固有値には明確な構造があります。
- $K$ 個の大きい固有値: $\lambda_1 \geq \cdots \geq \lambda_K > \sigma_n^2$。これらの固有値に対応する固有ベクトルは信号部分空間（signal subspace）$\mathcal{S}$ を張ります。
- $M – K$ 個の小さい固有値: $\lambda_{K+1} = \cdots = \lambda_M = \sigma_n^2$。これらの固有ベクトルは雑音部分空間（noise subspace）$\mathcal{N}$ を張ります。
信号部分空間はステアリングベクトル $\bm{a}(\theta_k)$（$k = 1, \ldots, K$）が張る空間と一致します。つまり
$$ \text{span}(\bm{e}_1, \ldots, \bm{e}_K) = \text{span}(\bm{a}(\theta_1), \ldots, \bm{a}(\theta_K)) $$
この事実から、重要な直交関係が導かれます。雑音部分空間の固有ベクトルはすべてのステアリングベクトルと直交するということです。
$$ \bm{e}_i^H \bm{a}(\theta_k) = 0 \quad (i = K+1, \ldots, M; \; k = 1, \ldots, K) $$
MUSICスペクトル
MUSIC法はこの直交性を利用してDOA推定を行います。雑音部分空間の固有ベクトルを列に持つ行列を $\bm{E}_n = [\bm{e}_{K+1}, \ldots, \bm{e}_M]$ とすると、MUSICスペクトル（MUSIC pseudo-spectrum）は
$$ P_{\text{MUSIC}}(\theta) = \frac{1}{\bm{a}^H(\theta)\bm{E}_n\bm{E}_n^H\bm{a}(\theta)} $$
で定義されます。
この式の分母 $\bm{a}^H(\theta)\bm{E}_n\bm{E}_n^H\bm{a}(\theta)$ は、ステアリングベクトル $\bm{a}(\theta)$ の雑音部分空間への射影のノルムの2乗です。
$\theta$ が信号源の方向 $\theta_k$ に一致すると、$\bm{a}(\theta_k)$ は雑音部分空間と直交するため分母がゼロ（理想的には）になり、$P_{\text{MUSIC}}$ は無限大のピークを示します。$\theta$ が信号源の方向から離れると、$\bm{a}(\theta)$ は雑音部分空間への射影を持つため、分母は正の値をとり、$P_{\text{MUSIC}}$ は有限の値になります。
MUSIC法の手順
MUSIC法の実行手順をまとめます。
- データ収集: $N$ スナップショットの受信データ $\bm{x}(t_n)$ を収集します
- 標本共分散行列の推定: $\hat{\bm{R}} = \frac{1}{N}\sum_{n=1}^{N}\bm{x}(t_n)\bm{x}^H(t_n)$
- 固有値分解: $\hat{\bm{R}}$ を固有値分解し、固有値を降順に並べます
- 信号数の推定: 固有値の段差から信号数 $K$ を推定します（MDL基準やAIC基準を使用）
- 雑音部分空間の構成: $K+1$ 番目以降の固有ベクトルから $\bm{E}_n$ を構成します
- MUSICスペクトルの計算: $\theta$ を走査して $P_{\text{MUSIC}}(\theta)$ を計算します
- ピーク検出: $P_{\text{MUSIC}}$ の $K$ 個の最大ピークの角度がDOA推定値です
MUSIC法の特性
MUSIC法は超解像（super-resolution）能力を持ちます。理論上、SNRが十分高くスナップショット数が十分多ければ、任意に接近した信号源を分離できます（レイリー限界を超えた分解能）。
ただし、以下の条件が必要です。
- 信号数 $K$ が素子数 $M$ より少ない: $K < M$ でないと雑音部分空間が存在しません
- 信号が無相関: コヒーレントな信号（マルチパス等）では雑音部分空間の構造が崩れます。空間スムージングなどの前処理が必要です
- 正確な信号数の推定: $K$ の推定を誤ると、雑音部分空間が正しく構成されず、偽のピークや見落としが生じます
MUSIC法の理論を理解したところで、Pythonシミュレーションで3つの手法を比較しましょう。
Pythonシミュレーション — 3手法の比較
2つの接近した信号源に対して、従来ビームフォーミング、Capon法、MUSIC法のDOA推定性能を比較します。
import numpy as np
import matplotlib.pyplot as plt
def steering_vector(theta, M, d, wavelength):
"""
ULAのステアリングベクトルを計算する。
Parameters
----------
theta : float 到来角 [rad]
M : int 素子数
d : float 素子間隔 [m]
wavelength : float 波長 [m]
"""
k = 2 * np.pi / wavelength
m = np.arange(M)
return np.exp(-1j * m * k * d * np.sin(theta))
def generate_array_data(thetas, M, d, wavelength, N, SNR_dB):
"""
ULAの受信データを生成する。
Parameters
----------
thetas : list 信号源の到来角 [rad]
M : int 素子数
d : float 素子間隔 [m]
wavelength : float 波長 [m]
N : int スナップショット数
SNR_dB : float SNR [dB]
"""
K = len(thetas)
A = np.column_stack([steering_vector(t, M, d, wavelength) for t in thetas])
# 信号（無相関ランダム信号）
S = (np.random.randn(K, N) + 1j * np.random.randn(K, N)) / np.sqrt(2)
# 雑音
sigma_n = 10**(-SNR_dB / 20)
noise = sigma_n * (np.random.randn(M, N) + 1j * np.random.randn(M, N)) / np.sqrt(2)
X = A @ S + noise
return X
# --- パラメータ設定 ---
M = 8 # 素子数
d = 0.5 # 素子間隔 (λ/2に正規化)
wavelength = 1.0
N = 200 # スナップショット数
SNR_dB = 20 # SNR [dB]
# 2つの信号源（5°差 — 従来法では分離困難）
theta1_deg, theta2_deg = 10, 15
theta1 = np.radians(theta1_deg)
theta2 = np.radians(theta2_deg)
# データ生成
np.random.seed(42)
X = generate_array_data([theta1, theta2], M, d, wavelength, N, SNR_dB)
# 標本共分散行列
R_hat = (X @ X.conj().T) / N
# 走査角度
theta_scan_deg = np.linspace(-90, 90, 1000)
theta_scan = np.radians(theta_scan_deg)
# --- (1) 従来ビームフォーミング ---
P_cbf = np.zeros(len(theta_scan))
for i, t in enumerate(theta_scan):
a = steering_vector(t, M, d, wavelength)
w = a / M
P_cbf[i] = np.real(w.conj() @ R_hat @ w)
P_cbf_dB = 10 * np.log10(P_cbf / np.max(P_cbf) + 1e-10)
# --- (2) Capon法 ---
R_inv = np.linalg.inv(R_hat)
P_capon = np.zeros(len(theta_scan))
for i, t in enumerate(theta_scan):
a = steering_vector(t, M, d, wavelength)
P_capon[i] = np.real(1.0 / (a.conj() @ R_inv @ a))
P_capon_dB = 10 * np.log10(P_capon / np.max(P_capon) + 1e-10)
# --- (3) MUSIC法 ---
eigenvalues, eigenvectors = np.linalg.eigh(R_hat)
# 固有値を降順にソート
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
K_est = 2 # 信号数（既知と仮定）
En = eigenvectors[:, K_est:] # 雑音部分空間
P_music = np.zeros(len(theta_scan))
for i, t in enumerate(theta_scan):
a = steering_vector(t, M, d, wavelength)
proj = a.conj() @ En @ En.conj().T @ a
P_music[i] = np.real(1.0 / (proj + 1e-10))
P_music_dB = 10 * np.log10(P_music / np.max(P_music) + 1e-10)
# --- 可視化 ---
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) 3手法の比較
ax = axes[0, 0]
ax.plot(theta_scan_deg, P_cbf_dB, '#45B7D1', linewidth=2, label='CBF (Conventional)')
ax.plot(theta_scan_deg, P_capon_dB, '#4ECDC4', linewidth=2, label='Capon (MVDR)')
ax.plot(theta_scan_deg, P_music_dB, '#FF6B6B', linewidth=2, label='MUSIC')
ax.axvline(x=theta1_deg, color='gray', linestyle='--', alpha=0.7, label=f'True DOA ({theta1_deg}°, {theta2_deg}°)')
ax.axvline(x=theta2_deg, color='gray', linestyle='--', alpha=0.7)
ax.set_xlabel('Angle [degrees]', fontsize=12)
ax.set_ylabel('Spatial Spectrum [dB]', fontsize=12)
ax.set_title(f'DOA Estimation: M={M}, SNR={SNR_dB}dB, N={N}', fontsize=13)
ax.set_xlim(-30, 50)
ax.set_ylim(-40, 3)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (b) 固有値プロット
ax = axes[0, 1]
ax.stem(range(1, M+1), 10*np.log10(eigenvalues), linefmt='c-',
markerfmt='co', basefmt='gray')
ax.axhline(y=10*np.log10(eigenvalues[K_est]),
color='red', linestyle='--', alpha=0.7,
label=f'Signal/Noise boundary (K={K_est})')
ax.set_xlabel('Eigenvalue Index', fontsize=12)
ax.set_ylabel('Eigenvalue [dB]', fontsize=12)
ax.set_title('Eigenvalue Spectrum', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (c) 信号源間隔を変化させた場合の分解能比較
ax = axes[1, 0]
separations = [3, 5, 8, 15]
for sep in separations:
t1, t2 = np.radians(10), np.radians(10 + sep)
X_temp = generate_array_data([t1, t2], M, d, wavelength, N, SNR_dB)
R_temp = (X_temp @ X_temp.conj().T) / N
ev, evec = np.linalg.eigh(R_temp)
idx_temp = np.argsort(ev)[::-1]
evec = evec[:, idx_temp]
En_temp = evec[:, 2:]
P_temp = np.zeros(len(theta_scan))
for i, t in enumerate(theta_scan):
a = steering_vector(t, M, d, wavelength)
proj = a.conj() @ En_temp @ En_temp.conj().T @ a
P_temp[i] = np.real(1.0 / (proj + 1e-10))
P_temp_dB = 10 * np.log10(P_temp / np.max(P_temp) + 1e-10)
ax.plot(theta_scan_deg, P_temp_dB, linewidth=1.5,
label=f'$\\Delta\\theta$ = {sep}°')
ax.set_xlabel('Angle [degrees]', fontsize=12)
ax.set_ylabel('MUSIC Spectrum [dB]', fontsize=12)
ax.set_title('MUSIC: Resolution vs Signal Separation', fontsize=13)
ax.set_xlim(-10, 40)
ax.set_ylim(-40, 3)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (d) SNR依存性
ax = axes[1, 1]
snr_values = [0, 5, 10, 20]
for snr in snr_values:
X_temp = generate_array_data([theta1, theta2], M, d, wavelength, N, snr)
R_temp = (X_temp @ X_temp.conj().T) / N
ev, evec = np.linalg.eigh(R_temp)
idx_temp = np.argsort(ev)[::-1]
evec = evec[:, idx_temp]
En_temp = evec[:, 2:]
P_temp = np.zeros(len(theta_scan))
for i, t in enumerate(theta_scan):
a = steering_vector(t, M, d, wavelength)
proj = a.conj() @ En_temp @ En_temp.conj().T @ a
P_temp[i] = np.real(1.0 / (proj + 1e-10))
P_temp_dB = 10 * np.log10(P_temp / np.max(P_temp) + 1e-10)
ax.plot(theta_scan_deg, P_temp_dB, linewidth=1.5,
label=f'SNR = {snr} dB')
ax.set_xlabel('Angle [degrees]', fontsize=12)
ax.set_ylabel('MUSIC Spectrum [dB]', fontsize=12)
ax.set_title('MUSIC: Effect of SNR', fontsize=13)
ax.set_xlim(-10, 40)
ax.set_ylim(-40, 3)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
シミュレーション結果から、3つの手法の分解能の違いが明確に確認できます。
- 3手法の比較（左上）: 5°の角度差を持つ2信号源に対して、従来ビームフォーミング（CBF）は1つの太いピークしか示さず、2信号源を分離できていません。Capon法は2つのピークが見え始めますがまだ不明瞭です。MUSIC法は2つの鋭いピークを明確に示し、高い分解能で2信号源を分離しています。
- 固有値スペクトル（右上）: 2つの大きな固有値と6つの小さな固有値の間に明確な段差があり、信号数 $K = 2$ が正しく判定できます。この段差が信号部分空間と雑音部分空間の境界です。
- 角度分離依存性（左下）: MUSIC法は3°の角度差でも2信号源をかろうじて分離できていますが、分解能はSNRとスナップショット数に依存します。15°以上離れていれば従来法でも十分分離可能です。
- SNR依存性（右下）: SNRが0 dBではMUSIC法でも分解能が大幅に低下し、ピークが不明瞭になります。SNR 10 dB以上で良好な分解能が得られ、20 dBでは非常に鋭いピークが観測されます。
Pythonによるスナップショット数の影響の分析
スナップショット数 $N$ がDOA推定精度に与える影響を調べます。
import numpy as np
import matplotlib.pyplot as plt
def steering_vector(theta, M, d, wavelength):
k = 2 * np.pi / wavelength
m = np.arange(M)
return np.exp(-1j * m * k * d * np.sin(theta))
def music_spectrum(R_hat, M, K, d, wavelength, theta_scan):
ev, evec = np.linalg.eigh(R_hat)
idx = np.argsort(ev)[::-1]
evec = evec[:, idx]
En = evec[:, K:]
P = np.zeros(len(theta_scan))
for i, t in enumerate(theta_scan):
a = steering_vector(t, M, d, wavelength)
proj = a.conj() @ En @ En.conj().T @ a
P[i] = np.real(1.0 / (proj + 1e-10))
return P
M = 8
d = 0.5
wavelength = 1.0
SNR_dB = 15
theta1, theta2 = np.radians(10), np.radians(15)
theta_scan_deg = np.linspace(-30, 50, 500)
theta_scan = np.radians(theta_scan_deg)
# 異なるスナップショット数
N_values = [10, 50, 200, 1000]
colors = ['#FF6B6B', '#FFD93D', '#4ECDC4', '#45B7D1']
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# (a) スナップショット数の影響
ax = axes[0]
np.random.seed(42)
for N_snap, color in zip(N_values, colors):
A = np.column_stack([steering_vector(theta1, M, d, wavelength),
steering_vector(theta2, M, d, wavelength)])
S = (np.random.randn(2, N_snap) + 1j*np.random.randn(2, N_snap)) / np.sqrt(2)
sigma_n = 10**(-SNR_dB/20)
noise = sigma_n * (np.random.randn(M, N_snap) + 1j*np.random.randn(M, N_snap)) / np.sqrt(2)
X = A @ S + noise
R_hat = (X @ X.conj().T) / N_snap
P = music_spectrum(R_hat, M, 2, d, wavelength, theta_scan)
P_dB = 10 * np.log10(P / np.max(P) + 1e-10)
ax.plot(theta_scan_deg, P_dB, color=color, linewidth=1.5,
label=f'N = {N_snap}')
ax.axvline(x=10, color='gray', linestyle='--', alpha=0.5)
ax.axvline(x=15, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Angle [degrees]', fontsize=12)
ax.set_ylabel('MUSIC Spectrum [dB]', fontsize=12)
ax.set_title('Effect of Snapshot Number on MUSIC', fontsize=13)
ax.set_xlim(-10, 40)
ax.set_ylim(-40, 3)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
# (b) DOA推定のRMSE vs スナップショット数（モンテカルロ）
ax = axes[1]
N_range = np.logspace(1, 3, 20).astype(int)
n_trials = 50
rmse_music = np.zeros(len(N_range))
rmse_capon = np.zeros(len(N_range))
for j, N_snap in enumerate(N_range):
errors_music = []
errors_capon = []
for trial in range(n_trials):
A = np.column_stack([steering_vector(theta1, M, d, wavelength),
steering_vector(theta2, M, d, wavelength)])
S = (np.random.randn(2, N_snap) + 1j*np.random.randn(2, N_snap)) / np.sqrt(2)
sigma_n = 10**(-SNR_dB/20)
noise = sigma_n * (np.random.randn(M, N_snap) + 1j*np.random.randn(M, N_snap)) / np.sqrt(2)
X = A @ S + noise
R_hat_mc = (X @ X.conj().T) / N_snap
# MUSIC
P_m = music_spectrum(R_hat_mc, M, 2, d, wavelength, theta_scan)
# 2つのピークを検出
peaks_idx = []
for ii in range(1, len(P_m)-1):
if P_m[ii] > P_m[ii-1] and P_m[ii] > P_m[ii+1]:
peaks_idx.append(ii)
if len(peaks_idx) >= 2:
peaks_sorted = sorted(peaks_idx, key=lambda x: P_m[x], reverse=True)[:2]
est_angles = sorted([theta_scan_deg[p] for p in peaks_sorted])
error = np.sqrt(((est_angles[0]-10)**2 + (est_angles[1]-15)**2) / 2)
errors_music.append(error)
rmse_music[j] = np.mean(errors_music) if errors_music else np.nan
ax.semilogx(N_range, rmse_music, 'c-o', linewidth=2, markersize=5, label='MUSIC')
ax.set_xlabel('Number of Snapshots $N$', fontsize=12)
ax.set_ylabel('RMSE [degrees]', fontsize=12)
ax.set_title('DOA Estimation Accuracy vs Snapshot Number', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.show()
スナップショット数の影響の分析から、以下の重要な知見が得られます。
- スナップショット数が少ない場合の劣化（左図）: $N = 10$（素子数 $M = 8$ よりわずかに多い）では、MUSICスペクトルのピークが広がり、分解能が低下しています。$N = 50$ 以上でピークが鋭くなり、$N = 200$ で十分安定します。
- RMSEの収束（右図）: スナップショット数の増加に伴い、DOA推定の二乗平均平方根誤差（RMSE）が減少します。$N > 100$ で RMSE が 1° 以下に収束しています。
- 実用上の指針: MUSIC法を実用的に使うためには、最低でも $N \geq 2M$（素子数の2倍）のスナップショットが必要です。精度を保証するには $N \geq 4M$ 程度が推奨されます。
まとめ
本記事では、適応ビームフォーミングの理論として、Capon法とMUSIC法によるDOA推定について解説しました。
- 従来ビームフォーミングの分解能限界: ビーム幅 $\approx \lambda/(Md)$ 以下の角度差を持つ信号源は分離できません。これは物理的な回折限界に相当します
- 信号モデル: ULAの受信信号はステアリングベクトルと信号の積として表され、共分散行列が信号の統計的性質を集約します
- Capon法（MVDR）: 注目方向の信号を歪みなく通過させつつ、出力パワーを最小化する適応的な重みを $\bm{w} = \bm{R}^{-1}\bm{a}/(\bm{a}^H\bm{R}^{-1}\bm{a})$ で計算します。従来法より高い分解能を持ちます
- MUSIC法: 共分散行列の固有値分解により信号部分空間と雑音部分空間を分離し、ステアリングベクトルと雑音部分空間の直交性を利用して超解像DOA推定を行います
- 分解能の比較: MUSIC法 > Capon法 > 従来法の順に高い分解能を持ちますが、MUSIC法はSNRとスナップショット数に敏感です
これらの手法は、レーダー、通信、音響など幅広い分野の基盤技術であり、現代の信号処理において最も重要な成果の一つです。
次のステップとして、以下の記事も参考にしてください。