複数のアンテナを並べたアレイアンテナに、電波がどの方向から飛来しているかを知りたい場面を考えます。レーダーで複数の目標を分離したい、5G基地局でユーザーの方向を特定してビームを向けたい、電波天文学で天体の位置を高精度に測りたい――こうした「到来方向推定（DOA: Direction Of Arrival estimation）」は、現代の無線システムの中核をなす技術です。
素朴に考えると、アンテナアレイにビームを向けて、受信電力が最大になる角度を探せばよさそうです。これが古典的なビームフォーマ法です。しかし、ビームフォーマには「角度分解能がアンテナ開口で決まるビーム幅に縛られる」という決定的な限界があります。2つの目標がビーム幅より近い角度にいると、1つの山にしか見えず分離できないのです。
ここで登場するのが本記事の主役、MUSIC（MUltiple SIgnal Classification）法です。MUSICは受信信号の統計的な構造――より正確には受信共分散行列の固有値分解――を巧みに利用し、ビーム幅という物理的な制約を突破して、近接した複数波を鋭く分離します。応用先は広く、フェーズドアレイレーダーの多目標分離、MIMO通信の到来角推定、ソナーや地震波解析、さらには電波天文学の干渉計まで、波の到来方向を測るあらゆる場面で使われています。
本記事の内容
- アレイ受信信号のモデル化とステアリングベクトル
- 受信共分散行列 $\bm{R} = E[\bm{x}\bm{x}^H]$ の構造
- 共分散行列の固有値分解と信号部分空間・雑音部分空間への分解
- ステアリングベクトルが雑音部分空間と直交することの導出
- MUSICスペクトル $P(\theta) = 1/(\bm{a}^H \bm{E}_n \bm{E}_n^H \bm{a})$ の導出
- なぜビームフォーマより高分解能になるのかの直感的理解
- Pythonによる近接2波・相関波に対するMUSICとビームフォーマの分解能比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
固有値分解、エルミート行列、共分散行列といった線形代数・統計の基礎が出てきますが、本文中で必要な性質はその都度説明します。
到来方向推定とは
夜、海岸に立って遠くの船の汽笛が「右の方」から聞こえるとき、私たちは無意識に左右の耳に音が届く時間差を使って方向を判断しています。アレイアンテナによる到来方向推定も、原理はこれと全く同じです。複数のアンテナ素子に電波が届くわずかな時間差（位相差）から、波がどの方向から来たかを逆算するのです。
イメージとしては、海面を一定方向に進む平面波を考えてください。波の山が斜めにアンテナの列を横切っていくとき、列の端のアンテナと反対側の端のアンテナとでは、波の山が到達するタイミングがずれます。このずれ（位相差）は到来角に一意に対応するので、各素子で観測した位相差のパターンを読めば、到来角が分かります。
問題は、複数の波が同時に飛んでくるときです。各素子で観測される信号は、すべての波の位相差パターンが重なり合ったものになります。この「混ざった」信号から、何個の波が、それぞれどの方向から来ているかを分離して取り出すことが到来方向推定の核心です。MUSICは、この分離を信号の統計的な構造を使って鮮やかに行います。
まずは、アレイが受け取る信号を数式でモデル化するところから始めましょう。
アレイ受信信号のモデル化
等間隔リニアアレイとステアリングベクトル
最も基本的な配置として、$M$ 個のアンテナ素子を間隔 $d$ で一直線に並べた等間隔リニアアレイ（ULA: Uniform Linear Array）を考えます。アレイの法線方向（正面）から測った到来角を $\theta$ とします。
遠方からの平面波がアレイに角度 $\theta$ で入射すると、隣り合う素子の間には伝搬距離の差 $d\sin\theta$ が生じます。波長を $\lambda$ とすると、この距離差は位相差
$$ \phi = \frac{2\pi}{\lambda} d \sin\theta $$
に対応します。波長 $\lambda$ あたり位相が $2\pi$ 回るので、距離 $d\sin\theta$ では $2\pi d\sin\theta / \lambda$ だけ位相が進む、という関係です。
第0素子を基準（位相0）にとると、第 $m$ 素子（$m = 0, 1, \dots, M-1$）では基準より $m\phi$ だけ位相が進みます。1つの到来波に対する各素子の応答をまとめたベクトルをステアリングベクトルと呼び、
$$ \bm{a}(\theta) = \begin{bmatrix} 1 \\ e^{j\phi} \\ e^{j2\phi} \\ \vdots \\ e^{j(M-1)\phi} \end{bmatrix}, \quad \phi = \frac{2\pi d \sin\theta}{\lambda} $$
と表します。$\bm{a}(\theta)$ は「角度 $\theta$ から来る波が、各素子にどんな位相パターンを刻むか」を表す $M$ 次元の複素ベクトルです。到来角 $\theta$ を変えると、このベクトルは $M$ 次元複素空間の中で連続的に動いていきます。この「ステアリングベクトルが角度ごとに描く軌跡」が、後でMUSICの鍵になります。
複数波の重ね合わせとスナップショット
いま、互いに異なる方向 $\theta_1, \theta_2, \dots, \theta_K$ から $K$ 個の波（信号源）が同時に到来しているとします。各波の複素振幅（時刻 $t$ における信号波形）を $s_1(t), \dots, s_K(t)$ とすると、線形性により、素子で観測される信号はすべての波の応答の重ね合わせになります。
第 $m$ 素子の受信信号 $x_m(t)$ は
$$ x_m(t) = \sum_{k=1}^{K} a_m(\theta_k) \, s_k(t) + n_m(t) $$
ここで $a_m(\theta_k)$ はステアリングベクトル $\bm{a}(\theta_k)$ の第 $m$ 成分、$n_m(t)$ は第 $m$ 素子の熱雑音です。これを全素子についてベクトルでまとめると、すっきりした行列形式になります。
$$ \bm{x}(t) = \bm{A}\,\bm{s}(t) + \bm{n}(t) $$
各記号の意味を整理します。
$$ \bm{x}(t) = \begin{bmatrix} x_0(t) \\ \vdots \\ x_{M-1}(t) \end{bmatrix}, \quad \bm{A} = [\bm{a}(\theta_1), \dots, \bm{a}(\theta_K)], \quad \bm{s}(t) = \begin{bmatrix} s_1(t) \\ \vdots \\ s_K(t) \end{bmatrix}, \quad \bm{n}(t) = \begin{bmatrix} n_0(t) \\ \vdots \\ n_{M-1}(t) \end{bmatrix} $$
$\bm{x}(t)$ は $M \times 1$ の観測ベクトル、$\bm{A}$ は $K$ 本のステアリングベクトルを横に並べた $M \times K$ のアレイ多様体行列（array manifold matrix）、$\bm{s}(t)$ は $K \times 1$ の信号ベクトル、$\bm{n}(t)$ は $M \times 1$ の雑音ベクトルです。ある時刻 $t$ で全素子を同時にサンプリングして得た1本の観測ベクトル $\bm{x}(t)$ をスナップショットと呼びます。
このモデルが教えてくれる重要な事実は、観測ベクトル $\bm{x}(t)$ の信号成分 $\bm{A}\bm{s}(t)$ は、$K$ 本のステアリングベクトル $\bm{a}(\theta_1), \dots, \bm{a}(\theta_K)$ の線形結合に過ぎないということです。つまり信号成分は、$M$ 次元空間の中で、これら $K$ 本のベクトルが張る高々 $K$ 次元の部分空間の中にしか存在できません。この事実が、後で「信号部分空間」という概念に結実します。
次は、このモデルから受信信号の統計量――共分散行列――を作り、その内部構造を調べていきましょう。
受信共分散行列の構造
共分散行列の定義
DOA推定で本質的なのは、1本のスナップショットそのものではなく、多数のスナップショットから得られる統計的な相関構造です。観測ベクトル $\bm{x}(t)$ の自己相関を表す $M \times M$ の受信共分散行列を
$$ \bm{R} = E[\bm{x}(t)\bm{x}^H(t)] $$
で定義します。$\bm{x}^H$ は $\bm{x}$ の共役転置（エルミート転置）、$E[\cdot]$ は時間に関する期待値（アンサンブル平均）です。$\bm{R}$ の $(i,j)$ 成分は第 $i$ 素子と第 $j$ 素子の信号の相関 $E[x_i \, x_j^*]$ であり、素子間でどれだけ似た信号が観測されているかを表します。
共分散行列を信号項と雑音項に分解する
モデル $\bm{x} = \bm{A}\bm{s} + \bm{n}$ を代入して $\bm{R}$ を展開しましょう。ゴールは、$\bm{R}$ を「信号が作る項」と「雑音が作る項」のきれいな和に分けることです。
$$ \bm{R} = E\!\left[(\bm{A}\bm{s} + \bm{n})(\bm{A}\bm{s} + \bm{n})^H\right] $$
共役転置を展開すると $(\bm{A}\bm{s}+\bm{n})^H = \bm{s}^H\bm{A}^H + \bm{n}^H$ なので、積を展開して期待値の線形性を使うと、4つの項に分かれます。
$$ \bm{R} = \bm{A}\,E[\bm{s}\bm{s}^H]\,\bm{A}^H + \bm{A}\,E[\bm{s}\bm{n}^H] + E[\bm{n}\bm{s}^H]\,\bm{A}^H + E[\bm{n}\bm{n}^H] $$
ここで2つの自然な仮定を置きます。第一に、信号と雑音は無相関である（$E[\bm{s}\bm{n}^H] = \bm{0}$）。雑音は受信機内部で発生する熱雑音であり、外から来る信号とは独立だからです。第二に、各素子の雑音は互いに無相関な白色雑音で、同じ電力 $\sigma^2$ を持つ（$E[\bm{n}\bm{n}^H] = \sigma^2 \bm{I}$）。すると中央の2項が消え、
$$ \bm{R} = \bm{A}\,\bm{R}_s\,\bm{A}^H + \sigma^2 \bm{I} $$
という基本式が得られます。ここで $\bm{R}_s = E[\bm{s}\bm{s}^H]$ は $K \times K$ の信号共分散行列です。$\bm{R}_s$ の対角成分は各信号源の電力、非対角成分は信号源どうしの相関を表します。信号源が互いに独立なら $\bm{R}_s$ は対角行列になります。
信号項の階数（ランク）
この $\bm{R} = \bm{A}\bm{R}_s\bm{A}^H + \sigma^2\bm{I}$ という分解こそ、MUSICの出発点です。第1項 $\bm{A}\bm{R}_s\bm{A}^H$ に注目しましょう。$\bm{A}$ は $M \times K$ 行列で、$M > K$（素子数が信号源数より多い）のとき、その階数は高々 $K$ です。$\bm{R}_s$ がフルランク（$K$）なら、$\bm{A}\bm{R}_s\bm{A}^H$ の階数はちょうど $K$ になります。
つまり信号項 $\bm{A}\bm{R}_s\bm{A}^H$ は $M \times M$ の行列でありながら、その「有効な次元」は $K$ しかありません。残りの $M – K$ 次元方向には信号のエネルギーが全く乗っていないのです。一方で雑音項 $\sigma^2\bm{I}$ は全方向に均等にエネルギーを乗せます。この「信号は一部の次元にしかいない／雑音は全方向に均一」という非対称性を、固有値分解で炙り出すのが次のステップです。
固有値分解と部分空間への分解
共分散行列の固有値分解
$\bm{R}$ はエルミート行列（$\bm{R}^H = \bm{R}$）です。実際、$(\bm{A}\bm{R}_s\bm{A}^H)^H = \bm{A}\bm{R}_s^H\bm{A}^H = \bm{A}\bm{R}_s\bm{A}^H$（$\bm{R}_s$ もエルミート）であり、$\sigma^2\bm{I}$ も明らかにエルミートだからです。エルミート行列は次の重要な性質を持ちます。
- すべての固有値は実数である
- 固有ベクトルは互いに直交するように選べる（正規直交基底をなす）
- 半正定値なので固有値は非負である
したがって $\bm{R}$ は実固有値 $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_M \geq 0$ と、対応する正規直交固有ベクトル $\bm{e}_1, \dots, \bm{e}_M$ を持ち、
$$ \bm{R} = \sum_{i=1}^{M} \lambda_i \, \bm{e}_i \bm{e}_i^H $$
とスペクトル分解できます。各 $\bm{e}_i$ は $\bm{R}\bm{e}_i = \lambda_i \bm{e}_i$ を満たす固有ベクトルです。
固有値が信号と雑音に分かれることの導出
ここからが核心です。固有値 $\lambda_i$ と固有ベクトル $\bm{e}_i$ が、信号と雑音にどう関係するかを導きます。
まず、信号部分空間の外、すなわち全ステアリングベクトルに直交する方向を考えます。$\bm{e}$ がすべてのステアリングベクトルに直交する（$\bm{A}^H\bm{e} = \bm{0}$）固有ベクトルだとすると、
$$ \bm{R}\bm{e} = \bm{A}\bm{R}_s\bm{A}^H\bm{e} + \sigma^2\bm{e} = \bm{A}\bm{R}_s(\bm{A}^H\bm{e}) + \sigma^2\bm{e} = \bm{A}\bm{R}_s\bm{0} + \sigma^2\bm{e} = \sigma^2\bm{e} $$
となります。途中で $\bm{A}^H\bm{e} = \bm{0}$ を使うと第1項が消えるのがポイントです。この式は「$\bm{e}$ が固有値 $\sigma^2$ の固有ベクトルである」ことを意味します。$\bm{A}$ の列空間に直交する方向は $M – K$ 次元あるので、固有値 $\sigma^2$ が $M – K$ 個現れることが分かります。
一方、信号部分空間の中にある固有ベクトルに対しては、第1項 $\bm{A}\bm{R}_s\bm{A}^H\bm{e}$ が正の寄与をします。$\bm{A}\bm{R}_s\bm{A}^H$ は階数 $K$ の半正定値行列なので、その値域に対応する $K$ 個の固有値は正です。これらの固有ベクトルに対する $\bm{R}$ の固有値は
$$ \lambda_i = (\text{信号項の固有値}) + \sigma^2 > \sigma^2 \quad (i = 1, \dots, K) $$
となり、必ず $\sigma^2$ より大きくなります。整理すると、固有値は次のようにきれいに二分されます。
$$ \underbrace{\lambda_1 \geq \dots \geq \lambda_K}_{> \sigma^2 \text{（信号＋雑音）}} > \underbrace{\lambda_{K+1} = \dots = \lambda_M = \sigma^2}_{\text{（雑音のみ）}} $$
大きい方の $K$ 個の固有値には信号のエネルギーが乗っており、小さい方の $M – K$ 個の固有値はすべて雑音電力 $\sigma^2$ に等しくなるのです。これは驚くべき事実です。固有値の大きさを見るだけで、何個の波が来ているか（信号源数 $K$）を推定できることを意味します。
信号部分空間と雑音部分空間
固有ベクトルを固有値の大小で2グループに分けます。
- 信号部分空間（signal subspace）: 大きい $K$ 個の固有値に対応する固有ベクトル $\bm{e}_1, \dots, \bm{e}_K$ が張る空間。これらをまとめた $M \times K$ 行列を $\bm{E}_s = [\bm{e}_1, \dots, \bm{e}_K]$ とします。
- 雑音部分空間（noise subspace）: 小さい $M – K$ 個の固有値に対応する固有ベクトル $\bm{e}_{K+1}, \dots, \bm{e}_M$ が張る空間。これらをまとめた $M \times (M-K)$ 行列を $\bm{E}_n = [\bm{e}_{K+1}, \dots, \bm{e}_M]$ とします。
固有ベクトルは正規直交基底をなすので、信号部分空間と雑音部分空間は互いに直交し、両者を合わせると $M$ 次元空間全体になります。「雑音部分空間」という名前ですが、雑音はあらゆる方向に存在するので、ここに信号がいないという意味で雑音的、と理解すると正確です。
ここまでで、共分散行列を固有値分解すると空間が2つの直交する部分空間に分かれることが分かりました。次は、この直交構造が到来角の情報をどう保持しているかを見ます。
ステアリングベクトルと雑音部分空間の直交性
直交性の導出
MUSICの心臓部となる命題を導きます。主張はこうです。真の到来角 $\theta_k$（$k=1,\dots,K$）におけるステアリングベクトル $\bm{a}(\theta_k)$ は、雑音部分空間に完全に直交する。 式で書くと
$$ \bm{E}_n^H \, \bm{a}(\theta_k) = \bm{0}, \quad k = 1, \dots, K $$
です。これを示しましょう。
まず、信号部分空間 $\bm{E}_s$ の張る空間と、アレイ多様体行列 $\bm{A}$ の列空間（ステアリングベクトルが張る空間）が一致することを確認します。前節で見たように、信号部分空間は信号項 $\bm{A}\bm{R}_s\bm{A}^H$ のエネルギーが乗る $K$ 次元空間でした。$\bm{A}\bm{R}_s\bm{A}^H$ の値域は $\bm{A}$ の列空間そのものなので（$\bm{R}_s$ がフルランクであれば）、
$$ \text{span}\{\bm{e}_1, \dots, \bm{e}_K\} = \text{span}\{\bm{a}(\theta_1), \dots, \bm{a}(\theta_K)\} $$
が成り立ちます。つまり信号部分空間は、真の到来角のステアリングベクトルがちょうど張る空間なのです。
次に、雑音部分空間は信号部分空間の直交補空間でした。したがって、信号部分空間に含まれる任意のベクトル――特に各ステアリングベクトル $\bm{a}(\theta_k)$――は、雑音部分空間のすべての基底ベクトル $\bm{e}_{K+1}, \dots, \bm{e}_M$ と直交します。これを行列でまとめれば
$$ \bm{E}_n^H \, \bm{a}(\theta_k) = \bm{0}, \quad k = 1, \dots, K $$
が得られます。導出完了です。雑音部分空間への射影がゼロになる、という極めて鋭い条件が、真の到来角でのみ成立するのです。
直交性を「角度のスキャン」に変える
この直交性を実際の角度推定に使うには、視点を逆転させます。私たちは固有値分解から $\bm{E}_n$ を計算できますが、真の到来角 $\theta_k$ は知りません。そこで、角度 $\theta$ をテスト用に連続的に動かしながら、各 $\theta$ のステアリングベクトル $\bm{a}(\theta)$ が雑音部分空間とどれだけ直交しているかを測ります。
雑音部分空間への射影のノルムの二乗
$$ \| \bm{E}_n^H \bm{a}(\theta) \|^2 = \bm{a}^H(\theta)\,\bm{E}_n\bm{E}_n^H\,\bm{a}(\theta) $$
を考えると、これは $\theta$ が真の到来角に一致したとき理論上ゼロになり、それ以外ではゼロにならない正の値をとります。$\bm{E}_n\bm{E}_n^H$ は雑音部分空間への射影行列です。
この量の逆数をとれば、真の到来角で「無限大に発散する」鋭いピークを持つ関数が作れます。次節でいよいよMUSICスペクトルを定義します。
MUSICスペクトルの導出
MUSICスペクトルの定義
前節で導いた直交性 $\bm{E}_n^H\bm{a}(\theta_k) = \bm{0}$ をピーク探索に変換するため、雑音部分空間への射影のノルム二乗の逆数をMUSICスペクトルとして定義します。
$$ 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{E}_n^H\bm{a}(\theta)\|^2$ は、テスト角 $\theta$ のステアリングベクトルを雑音部分空間に射影したベクトルの長さの二乗です。
- $\theta$ が真の到来角 $\theta_k$ のとき、直交性より分母が（理論上）ゼロになり、$P_{\text{MUSIC}}(\theta) \to \infty$ と発散する
- $\theta$ が真の到来角からずれていると、分母は正の値をとり、$P_{\text{MUSIC}}(\theta)$ は有限の小さな値になる
したがって、$P_{\text{MUSIC}}(\theta)$ を $\theta$ について計算してプロットすると、真の到来角 $\theta_1, \dots, \theta_K$ にちょうど $K$ 個の鋭いピークが立ちます。ピークの位置を読み取れば到来角が推定できる、というのがMUSICの全体像です。
なお、$P_{\text{MUSIC}}(\theta)$ は受信電力そのものを表すわけではないことに注意してください。これは「角度ごとの雑音部分空間との直交の度合い」を可視化した指標であり、その意味で「擬スペクトル（pseudo-spectrum）」と呼ばれます。ピークの高さは信号電力ではなく、推定の鋭さ（直交性の良さ）を反映します。
実データでの計算手順
実際には期待値 $\bm{R} = E[\bm{x}\bm{x}^H]$ は得られないので、$N$ 個のスナップショット $\bm{x}(t_1), \dots, \bm{x}(t_N)$ から標本共分散行列で近似します。
$$ \hat{\bm{R}} = \frac{1}{N}\sum_{n=1}^{N} \bm{x}(t_n)\,\bm{x}^H(t_n) $$
スナップショット数 $N$ が大きいほど $\hat{\bm{R}}$ は真の $\bm{R}$ に近づきます。MUSICの計算手順をまとめると次のようになります。
- $N$ 個のスナップショットから標本共分散行列 $\hat{\bm{R}}$ を作る
- $\hat{\bm{R}}$ を固有値分解し、固有値を降順に並べる
- 固有値の落差から信号源数 $K$ を推定する（または既知とする）
- 小さい $M-K$ 個の固有値に対応する固有ベクトルから雑音部分空間 $\bm{E}_n$ を構成する
- テスト角 $\theta$ をスキャンしながら $P_{\text{MUSIC}}(\theta)$ を計算する
- ピーク位置を到来角の推定値とする
有限サンプルでは分母は厳密にはゼロにならないので、ピークは無限大ではなく「非常に高い鋭い山」になります。それでも分離能は古典法を大きく上回ります。
ここまでで、MUSICがなぜ機能するかの理論はそろいました。では、肝心の「なぜビームフォーマより高分解能なのか」を、直交性の観点から整理しましょう。
なぜビームフォーマより高分解能なのか
ビームフォーマ法の復習
比較対象として、古典的なビームフォーマ法（遅延和ビームフォーマ、Bartlett法）を振り返ります。ビームフォーマは、テスト角 $\theta$ にビームを向けたときの出力電力
$$ P_{\text{BF}}(\theta) = \bm{a}^H(\theta)\,\bm{R}\,\bm{a}(\theta) $$
を角度に対してスキャンし、電力が最大になる角度を到来角とみなします。これは「アレイの感度を $\theta$ 方向に集中させ、そこから来る電力を測る」という直感的な方法です。
分解能の本質的な違い
ビームフォーマの分解能は、メインビーム幅で決まります。ビーム幅はアレイの開口（素子数 $M$ と間隔 $d$ の積）に反比例し、おおよそ $\lambda / (Md)$ ラジアン程度です。2つの到来波がこのビーム幅より近いと、それぞれのビームパターンが重なり合い、出力電力スペクトルは1つの大きな山に融合してしまいます。ビームフォーマは信号エネルギーを「角度方向に広げて」測るため、近接した波を分離できないのです。これは開口で決まる物理的な制約であり、SN比をいくら上げても本質的には改善しません。
一方MUSICは、電力ではなく部分空間の直交性を使います。鍵は、雑音部分空間が真のステアリングベクトルと「厳密に直交する」という性質です。テスト角 $\theta$ が真の到来角からほんのわずかでもずれると、$\bm{a}(\theta)$ は雑音部分空間と直交しなくなり、$\|\bm{E}_n^H\bm{a}(\theta)\|^2$ は急速に増加します。逆数をとった $P_{\text{MUSIC}}$ は真の角度で発散し、隣の真の角度との間では深く落ち込みます。この「直交か非直交か」という条件は二値的に鋭く、ビーム幅のような連続的な広がりを持ちません。
結果として、MUSICは2つの波がビーム幅よりはるかに近くても、SN比とスナップショット数が十分なら別々の鋭いピークとして分離できます。このため、MUSICのような部分空間法は超分解能（super-resolution）法と呼ばれます。違いを一言でまとめると、ビームフォーマは「信号エネルギーの山を探す」のに対し、MUSICは「信号部分空間に直交する穴（雑音部分空間との直交点）を探す」のです。穴の方が山よりはるかに鋭く位置を特定できる、というわけです。
MUSICの前提と限界
ただし、MUSICの超分解能は無条件ではありません。重要な前提として、信号源が互いに無相関（または弱相関）であり、信号共分散行列 $\bm{R}_s$ がフルランクであることが要求されます。マルチパスのように波が互いに強く相関すると $\bm{R}_s$ がランク落ちし、信号部分空間の次元が真の信号源数より小さくなって、直交性の前提が崩れます。この場合、空間平滑化（spatial smoothing）という前処理で相関を緩和する必要があります。後のPython実装で、相関波がMUSICを劣化させる様子と空間平滑化の効果を実際に確認します。
理論が出そろったので、ここからは具体的な数値例とPython実装で、MUSICの威力を体感していきましょう。
具体例: 2素子・1波の場合
抽象的な議論を具体化するため、最小構成の $M=2$ 素子、$K=1$ 波の例で固有値分解を手で追ってみます。半波長間隔（$d = \lambda/2$）のアレイに、角度 $\theta_1$ から電力 $p$ の波が来ているとします。$d = \lambda/2$ のとき位相差は $\phi = \pi\sin\theta_1$ なので、ステアリングベクトルは
$$ \bm{a}(\theta_1) = \begin{bmatrix} 1 \\ e^{j\phi} \end{bmatrix} $$
です。信号共分散は $\bm{R}_s = [p]$（スカラー）、雑音電力を $\sigma^2$ とすると、共分散行列は
$$ \bm{R} = p\,\bm{a}(\theta_1)\bm{a}^H(\theta_1) + \sigma^2\bm{I} = p\begin{bmatrix} 1 & e^{-j\phi} \\ e^{j\phi} & 1 \end{bmatrix} + \sigma^2\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} $$
となります。$\bm{a}(\theta_1)\bm{a}^H(\theta_1)$ は階数1の行列で、その固有値は $\|\bm{a}(\theta_1)\|^2 = 2$ と $0$ です。これに $\sigma^2\bm{I}$ を足すと全固有値が $\sigma^2$ だけ持ち上がるので、$\bm{R}$ の固有値は
$$ \lambda_1 = 2p + \sigma^2, \quad \lambda_2 = \sigma^2 $$
になります。大きい方 $\lambda_1$ が信号＋雑音、小さい方 $\lambda_2 = \sigma^2$ が雑音のみ、という前章の一般論がそのまま現れています。固有値の比 $\lambda_1/\lambda_2 = 1 + 2p/\sigma^2$ はSN比に直結し、信号が強いほど2つの固有値の差が開きます。
$\lambda_2 = \sigma^2$ に対応する固有ベクトル $\bm{e}_2$（雑音部分空間）は $\bm{a}(\theta_1)$ に直交し、$\bm{e}_2 = \frac{1}{\sqrt{2}}[1, -e^{j\phi}]^T$ ととれます。実際に内積を計算すると
$$ \bm{e}_2^H\bm{a}(\theta_1) = \frac{1}{\sqrt{2}}\left(1\cdot 1 + (-e^{-j\phi})\cdot e^{j\phi}\right) = \frac{1}{\sqrt{2}}(1 – 1) = 0 $$
となり、確かに真の到来角で直交します。MUSICスペクトルの分母 $|\bm{e}_2^H\bm{a}(\theta)|^2$ は $\theta=\theta_1$ でゼロ、それ以外で正となり、$P_{\text{MUSIC}}$ に鋭いピークが立つことがこの最小例からも見て取れます。
この手計算の感触を、より現実的な多素子・複数波の設定でコードによって確かめましょう。
Pythonでの実装
アレイ受信信号の合成
まず、ULAの受信スナップショットをSN比つきで合成する関数を作ります。複素ガウス信号と複素ガウス雑音を用いて、指定した角度から指定した電力の波を到来させます。
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(0)
def steering_vector(theta_deg, M, d_over_lambda=0.5):
"""到来角theta_deg[度]に対するULAのステアリングベクトルを返す"""
theta = np.deg2rad(theta_deg)
m = np.arange(M) # 素子インデックス 0..M-1
phi = 2 * np.pi * d_over_lambda * np.sin(theta)
return np.exp(1j * m * phi) # 形状(M,)
def generate_snapshots(angles_deg, M, N, snr_db, d_over_lambda=0.5,
correlated=False):
"""複数波の受信スナップショット行列X (M x N) を生成する"""
K = len(angles_deg)
# アレイ多様体行列 A (M x K)
A = np.column_stack([steering_vector(a, M, d_over_lambda)
for a in angles_deg])
# 信号波形 S (K x N): 複素ガウス, 単位電力
if correlated:
# 1つの信号源を共通成分として相関を持たせる
base = (rng.standard_normal(N) + 1j*rng.standard_normal(N)) / np.sqrt(2)
S = np.array([base * np.exp(1j*0.3*k) for k in range(K)])
else:
S = (rng.standard_normal((K, N)) + 1j*rng.standard_normal((K, N))) / np.sqrt(2)
# 雑音電力をSN比から決定 (信号電力=1基準)
noise_power = 10 ** (-snr_db / 10)
Ncomp = np.sqrt(noise_power/2) * (rng.standard_normal((M, N))
+ 1j*rng.standard_normal((M, N)))
X = A @ S + Ncomp # 受信信号 (M x N)
return X
この関数は、引数 angles_deg で指定した各角度から到来する波を重ね合わせ、snr_db で指定したSN比の白色雑音を加えた受信信号行列を返します。correlated=True のときは全信号源を共通成分から作って強い相関を持たせ、後でMUSICが劣化する様子を再現します。これで「真の到来角が分かっている合成データ」を自由に作れるようになりました。
共分散行列の固有値構造を見る
理論で導いた「固有値が信号と雑音の2グループに分かれる」ことを、実際の合成データで確かめます。
import numpy as np
import matplotlib.pyplot as plt
# 設定: 8素子ULA, 真の到来角 -10度と20度の2波, SNR=10dB
M = 8
true_angles = [-10.0, 20.0]
K = len(true_angles)
N = 200 # スナップショット数
X = generate_snapshots(true_angles, M, N, snr_db=10)
# 標本共分散行列 R_hat = (1/N) X X^H
R = X @ X.conj().T / N
# 固有値分解 (エルミート行列なのでeighを使用, 固有値は昇順)
eigvals, eigvecs = np.linalg.eigh(R)
eigvals = eigvals[::-1] # 降順に並べ替え
plt.figure(figsize=(8, 5))
plt.stem(np.arange(1, M+1), eigvals)
plt.axhline(eigvals[K:].mean(), color='red', linestyle='--',
label='noise level (mean of small eigenvalues)')
plt.xlabel('Eigenvalue index (descending)')
plt.ylabel('Eigenvalue')
plt.title('Eigenvalues of the sample covariance matrix (M=8, K=2)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print("eigenvalues:", np.round(eigvals, 3))
このグラフから、固有値が明確に2つのグループに分かれることが読み取れます。大きい上位2個（$=K$）の固有値が信号部分空間に対応し、残りの6個（$=M-K$）はほぼ同じ小さな値（雑音電力 $\sigma^2$ に相当する赤い破線）に揃っています。理論が予言したとおり「大きい $K$ 個＋小さい $M-K$ 個」の構造が現れており、固有値の落差を見るだけで信号源数が2だと判定できることが分かります。
MUSICスペクトルとビームフォーマスペクトルの計算
雑音部分空間からMUSICスペクトルを、共分散行列から直接ビームフォーマスペクトルを計算し、角度スキャンを行います。
import numpy as np
def music_spectrum(R, K, M, scan_deg, d_over_lambda=0.5):
"""共分散行列RからMUSIC擬スペクトルを計算する"""
eigvals, eigvecs = np.linalg.eigh(R)
# 昇順なので小さい方 M-K 本が雑音部分空間
En = eigvecs[:, :M-K] # (M x (M-K))
P = np.zeros(len(scan_deg))
for i, ang in enumerate(scan_deg):
a = steering_vector(ang, M, d_over_lambda)
# 分母 = ||En^H a||^2
proj = En.conj().T @ a
P[i] = 1.0 / np.real(proj.conj() @ proj)
return P
def beamformer_spectrum(R, M, scan_deg, d_over_lambda=0.5):
"""Bartlettビームフォーマの出力電力スペクトルを計算する"""
P = np.zeros(len(scan_deg))
for i, ang in enumerate(scan_deg):
a = steering_vector(ang, M, d_over_lambda)
P[i] = np.real(a.conj() @ R @ a)
return P
music_spectrum は固有値分解で雑音部分空間 En を取り出し、各スキャン角について $\|\bm{E}_n^H\bm{a}(\theta)\|^2$ の逆数を計算します。beamformer_spectrum は $\bm{a}^H\bm{R}\bm{a}$ をそのまま評価します。この2つを同じデータに適用して比較できる準備が整いました。
近接2波に対する分解能の比較
MUSICの真価が問われる「ビーム幅以下に近接した2波」を分離できるか試します。8素子ULAでは半値ビーム幅はおよそ12〜13度なので、真の到来角を5度と11度（差6度）に設定して、両手法を比べます。
import numpy as np
import matplotlib.pyplot as plt
M = 8
true_angles = [5.0, 11.0] # ビーム幅より近い2波
K = len(true_angles)
N = 300
scan = np.linspace(-40, 40, 1601)
X = generate_snapshots(true_angles, M, N, snr_db=10)
R = X @ X.conj().T / N
P_music = music_spectrum(R, K, M, scan)
P_bf = beamformer_spectrum(R, M, scan)
# 見やすさのためdB正規化
P_music_db = 10*np.log10(P_music / P_music.max())
P_bf_db = 10*np.log10(P_bf / P_bf.max())
plt.figure(figsize=(10, 5))
plt.plot(scan, P_music_db, label='MUSIC', linewidth=2)
plt.plot(scan, P_bf_db, label='Beamformer (Bartlett)', linewidth=2)
for a in true_angles:
plt.axvline(a, color='gray', linestyle='--', alpha=0.7)
plt.xlabel('Angle [deg]')
plt.ylabel('Normalized pseudo-spectrum [dB]')
plt.title('MUSIC vs Beamformer: two closely spaced sources (5deg, 11deg)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-40, 2)
plt.tight_layout()
plt.show()
このグラフが本記事のハイライトです。灰色の破線が真の到来角（5度と11度）です。MUSICは2つの角度にきれいに分離した2本の鋭いピークを立てており、真値とほぼ一致しています。一方ビームフォーマは2つの波を分離できず、両者の中間あたりに1つのなだらかな山があるだけで、2波があることすら判別できません。これがビーム幅の壁です。MUSICが部分空間の直交性を使うことで、開口で決まる物理的な分解能限界を突破していることが視覚的に確認できます。
SN比とスナップショット数の影響
MUSICの分解能は無条件ではなく、SN比とスナップショット数に依存します。同じ近接2波を、SN比を変えながら推定してみます。
import numpy as np
import matplotlib.pyplot as plt
M = 8
true_angles = [5.0, 11.0]
K = 2
N = 300
scan = np.linspace(-40, 40, 1601)
snr_list = [-5, 0, 5, 15]
plt.figure(figsize=(10, 5))
for snr in snr_list:
X = generate_snapshots(true_angles, M, N, snr_db=snr)
R = X @ X.conj().T / N
P = music_spectrum(R, K, M, scan)
P_db = 10*np.log10(P / P.max())
plt.plot(scan, P_db, label=f'SNR = {snr} dB', linewidth=1.8)
for a in true_angles:
plt.axvline(a, color='gray', linestyle='--', alpha=0.6)
plt.xlabel('Angle [deg]')
plt.ylabel('Normalized MUSIC spectrum [dB]')
plt.title('Effect of SNR on MUSIC resolution (5deg, 11deg)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-40, 2)
plt.tight_layout()
plt.show()
このグラフから、SN比が下がるほどMUSICのピークが鈍り、分離が難しくなることが読み取れます。SN比15dBや5dBでは2本のピークが明瞭に分かれていますが、SN比が-5dBまで落ちると2つのピークが融合しかけ、谷が浅くなります。雑音が強いと標本共分散行列の推定誤差が増え、雑音部分空間の推定がぶれて直交性が崩れるためです。MUSICの超分解能は「十分なSN比とスナップショット数」という前提のもとで発揮される、という重要な但し書きをこの結果は教えてくれます。
相関波の問題と空間平滑化
最後に、理論パートで触れた「信号が相関するとMUSICが劣化する」現象と、その対策である空間平滑化を確認します。空間平滑化は、$M$ 素子アレイを重なり合う複数のサブアレイに分割し、各サブアレイの共分散行列を平均することで、相関による $\bm{R}_s$ のランク落ちを回復する前処理です。
import numpy as np
import matplotlib.pyplot as plt
def spatial_smoothing(X, M, L):
"""前方空間平滑化: サイズLのサブアレイで共分散を平均する"""
Msub = M - L + 1 # サブアレイ数
Rss = np.zeros((L, L), dtype=complex)
N = X.shape[1]
for i in range(Msub):
Xi = X[i:i+L, :] # i番目のサブアレイ (L x N)
Rss += Xi @ Xi.conj().T / N
return Rss / Msub
M = 10
true_angles = [0.0, 12.0]
K = 2
N = 300
scan = np.linspace(-40, 40, 1601)
# 相関の強い2波を生成
X = generate_snapshots(true_angles, M, N, snr_db=15, correlated=True)
R_plain = X @ X.conj().T / N
P_plain = music_spectrum(R_plain, K, M, scan)
# 空間平滑化後 (サブアレイサイズ L=7)
L = 7
R_ss = spatial_smoothing(X, M, L)
P_ss = music_spectrum(R_ss, K, L, scan)
plt.figure(figsize=(10, 5))
plt.plot(scan, 10*np.log10(P_plain/P_plain.max()),
label='MUSIC (correlated, no smoothing)', linewidth=2)
plt.plot(scan, 10*np.log10(P_ss/P_ss.max()),
label='MUSIC + spatial smoothing', linewidth=2)
for a in true_angles:
plt.axvline(a, color='gray', linestyle='--', alpha=0.6)
plt.xlabel('Angle [deg]')
plt.ylabel('Normalized MUSIC spectrum [dB]')
plt.title('Correlated sources: plain MUSIC vs spatial smoothing')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-40, 2)
plt.tight_layout()
plt.show()
このグラフから、相関波に対する素のMUSICの限界と空間平滑化の効果が読み取れます。平滑化なしのMUSIC（青）は、2波が強く相関しているため信号共分散 $\bm{R}_s$ がランク落ちし、2つのピークが正しく立たず到来角を見誤っています。一方、空間平滑化を施したMUSIC（オレンジ）は相関が緩和されて $\bm{R}_s$ のフルランク性が回復し、真の到来角（0度と12度）に2本のピークが復活しています。マルチパス環境のように波が相関する実用場面では、空間平滑化のような前処理が不可欠であることが分かります。
これらの実験を通じて、MUSICの威力と前提条件の両方を実感できたはずです。最後に要点を整理しましょう。
まとめ
本記事では、MUSIC法による到来方向推定の理論と導出、Python実装を解説しました。
- 受信信号モデル: アレイの観測ベクトルは $\bm{x}(t) = \bm{A}\bm{s}(t) + \bm{n}(t)$ と書け、信号成分はステアリングベクトルが張る $K$ 次元部分空間にのみ存在する
- 共分散行列の分解: $\bm{R} = \bm{A}\bm{R}_s\bm{A}^H + \sigma^2\bm{I}$ と信号項と雑音項に分かれ、信号項は階数 $K$ である
- 固有値の二分: $\bm{R}$ の固有値分解で、大きい $K$ 個（信号＋雑音）と小さい $M-K$ 個（雑音のみ $=\sigma^2$）に分かれ、固有値の落差から信号源数が分かる
- 直交性: 真の到来角のステアリングベクトルは雑音部分空間に直交する（$\bm{E}_n^H\bm{a}(\theta_k)=\bm{0}$）
- MUSICスペクトル: $P_{\text{MUSIC}}(\theta) = 1/(\bm{a}^H\bm{E}_n\bm{E}_n^H\bm{a})$ は真の到来角で発散する鋭いピークを持つ擬スペクトル
- 超分解能: ビームフォーマが信号エネルギーの山（ビーム幅に縛られる）を探すのに対し、MUSICは直交性という鋭い条件を使うため、ビーム幅以下の近接2波を分離できる
- 前提と限界: 信号源が無相関で $\bm{R}_s$ がフルランクであることが必要。相関波では空間平滑化などの前処理が有効
MUSICは部分空間法の代表格であり、ここから発展してESPRIT（回転不変性を使い、スキャン不要で直接到来角を求める）、Root-MUSIC（多項式の根として角度を解く）、最尤推定法などへとつながります。アレイ信号処理の世界を体系的に理解する出発点として、本記事の固有値分解と部分空間の考え方は今後も繰り返し登場します。
次のステップとして、以下の記事も参考にしてください。