MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

高分解能到来方向推定の仕組み
これでわかる MUSIC 法と ESPRIT 法
唐沢好男

一列等間隔に並んだ M 素子アレーアンテナがある。このアレーアンテナに、信号源も到来方
向も異なる L 個（L<M）の電波が入射している。このアレーアンテナの受信信号を解析して、電波の到来方向を高精度に求めるにはどうしたらよいかと言う問題を考える。先人たちはこれに取り組み、信号処理の傑作と言われる二つの方法を編み出している。MUSIC 法[1]と ESPRIT 法[2]である。共にアレー出力の相関行列の固有値解析を行うのであるが、面白いのは目の付け所が正反対であるという点である。一方（MUSIC）は、信号を除いた雑音だけのベクトル空間（雑音部分空間）を、一方（ESPRIT）は、雑音を除いた信号だけのベクトル空間（信号部分空間）を探索することによって、それぞれ目的とするものを手に入れている。この二つの手法を対比しながら、それぞれの仕組みをできるだけ分かりやすく（＝目から鱗が落ちたと言う気持ちが湧いてくるように）解説してみたい。
なお、両手法やその後の到来方向推定法の発展については、丁寧にまとめられた解説本も出ており[3],[4]、本稿作成においても参考にしている。また、本技術レポートシリーズの“TR-YK-006：
空間信号処理アンテナ”[5]で MUSIC 法を取り上げており、その部分とは重複がある。

１．命題
図１に示す素子数 M の等間隔直線状アレーアンテナ（uniform linear array: ULA）を考える。信号が独立な送信源（完全な無相関信号でなくてもよい）が L 個あり、それぞれが角度l（l=1, 2,…, L）で入射している。条件は L<M-1。受信機雑音を含んだアレー出力信号 xm（m=1, 2, …, M）
を解析して、電波の到来方向l を精度良く求めたい。これが命題である。（注：本稿ではこれを命題とするが、MUSIC 法や ESPRIT 法では、より広い条件に適用できる。
）

２．基本方針
最も簡単な例として、１波を２素子のアレーで受信することを考える。この場合、到来方向推定には、大きく分けて以下の三つの方法がある。
方法１：一方のアンテナに移相器を入れて位相を制御し、受信強度が最大になる位相から到来方向を求める（フェーズドアレー：ピークサーチ）

1

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

図１ 本稿での命題：送信源と到来方向が異なる L 個の電波を M 素子等間隔線形アレー（L<M-1）で受信し、その受信信号から到来方向を高精度に求めるには？

方法２：方法１と同じアンテナの構成で移相器の位相を制御し、受信強度が 0 になる位相から到来方向を求める（フェーズドアレー：ヌルサーチ）
方法３：二つの素子の受信信号の位相差を求め、アンテナ間隔と周波数により到来方向を計算により求める（干渉計方式）
方法１と方法２は同じフェーズドアレーで実現できるが、アンテナパターンを書いてみると分かるように、主ビームのピーク付近は位相停留点（移相器の位相を微少量ずらしても出力が変化しない点）になっていて、方向の推定精度が甘くなる。一方、方法２のヌルサーチでは、ヌル
点の落ち込みは急峻であるので、方向推定の精度が非常によい。図２はこのイメージである。この比較から、到来方向推定には方法１より方法２が勝ると言える。この考えの下で、方法２を多数の到来波が存在する環境に一般化した方法が MUSIC（Multiple Signal Classification）法であり、考案者は Ralph O. Schmidt [1] (1986; 発明当時はスタンフォード大学)である。
方法２をベースとする MUSIC 法は画期的な方法ではあったが、それにもいくつかの弱点がある。最大の問題は到来方向を求めるのに、最後の段階において、角度空間をサーチする必要があり、求めたい到来方向が数値計算によって直接に出てこないことにある。
（もう一つ、キャリブ
レーション問題も有るが、これは後述する）。そこで、到来方向を計算によって求めることができる方法３に着目し、その一般化を行って得た方法が ESPRIT（Estimation of Signal Parameter via
Rotational Invariance Techniques）法であり、考案者は同じスタンフォード大学の Richard Roy とThomas Kailath [2] (1989)である。
現時点において、MUSIC 法が ESPRIT 法に取って代わられているわけではなく、それぞれの特徴が活かされる分野で使い分けられている。

2

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

図２

方法１（ピークサーチ）と方法２（ヌルサーチ）

３．アレーアンテナの相関行列と固有値
次節以降において MUSIC 法と ESPRIT 法の原理を述べる。そのためには、数学的準備、すな
わちアレーアンテナの相関行列と固有値の理解が必要になる。ここでは、両方式に共通する部分をまとめる。
アレーアンテナでの受信信号はベクトルと行列を用いて表すと簡易な形になる。本稿では、ベクトルを英文字やギリシャ文字の小文字の太字、行列を大文字の太字で、ベクトルや行列の転置を上付き文字 T、複素共役転置を H で表す。アレーアンテナの受信信号（図１）は次式である。

受信信号x(t ) = As(t ) + n(t )(1)

受信信号ベクトル： x = ( x1
信号源ベクトル： s = ( s1
雑音ベクトル： n = (n1方向行列： A = ( a1x2xM ) T

s2

sL )T

n2

nM )T
aL )

a2

方向ベクトル（スティアリングベクトル）：al = ( e j (l ,d1 )e j (l ,d2 )e j (l ,d M ) )T

素子間隔d の ULA の場合

 ( , d ) = −kd cos →  (l , d m ) = −k (m − 1)d cosl3

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

雑音電力： ni2

信号源電力： sl= PN , ni ni* = 0 (i, i = 1,2,

, M ; i  i)= PSl , sl sl* = 0 (l , l  = 1,2,

, L; l  l )2

（ここでは、無相関としているが、一部に独立な信号が含まれていればよい、すなわち、完全相関以外の有相関であってよい）。
相関行列
相関行列は以下の式で与えられる。

Rxx = x(t ) x H (t ) = ( As(t ) + n(t ) )( As(t ) + n(t ) )H= ASAH + PN I M
S = diag ( PS1

(2)

PS 2

PSL ) （信号源間無相関の場合）

相関行列は空間情報（＝通信路情報）の統計的性質が（ほぼ）100%詰まっている宝の山である。ここから、どう財宝を掘り出すかが腕の見せ所になる。

固有値と固有ベクトル
相関行列 Rxx に対する以下の固有値方程式から、固有値とそれに属する固有ベクトル e を求める。

Rxx e =  e(3)

相関行列は、ij 成分と ji 成分が複素共役の関係にあるエルミート行列であり、かつ、対角要素が全て非負の実数であるので、M 個存在する固有値も全て非負の実数になる。この固有値を大
きい順に1～M とし、それぞれの固有値に属する固有ベクトルを e1～eM とする。雑音電力が 0（PN =0）のときは、L 番目までの固有値が正の値を持ち、L+1～M 番目の固有値は 0 になる。すなわち Rxx の行列のランクは L である。このことから、入射波の電力は、L 番目までの固有値に全て取り込まれていることになる。
固有ベクトルのうち、上位 L 個の固有値に属する固有ベクトル（e1～eL）は雑音の有無によらない。残りの M-L 個の固有ベクトルは、固有値の値が同じなので（いわゆる縮退しているので）、固有ベクトルが一意に定まらないが、実際には、測定値から得られる相関行列なので、完全縮退にはなっておらず、M-L 個の直交ベクトルを定めることができ、問題は無い。
信号電力と雑音電力の差が大きい場合には、信号源数 L の判定は容易であるが、実際の環境では、受信強度の弱い信号源もあり、判定が難しい場合もある。そのような場合には、統計的な評価手法（例えば AIC のような）を用いて推測することになる。ここでは、波数 L は正しく求4

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

められているとして議論を進める。
以下、固有値の関係式を整理してまとめておく。

Rxx = E  E H(4)

 = diag ( 1 2

1  2 M ) L  L+1 =

E = ( e1 e2

= M

eM )

相関行列を信号成分と雑音成分に分解すると
M

Rxx = Rss + Rnn =  i ei eiH

(5)

i =1

L

Rss = ASAH =  Si ei eiH , Rnn = PN I M
i =1

雑音がない場合（PN=0）の固有値は

 Si
 0

(1  i  L)

i = 

(L + 1  i  M )雑音を含めた固有値は

 Si + PN
 PN

(1  i  L)
(L + 1  i  M )i = 

(6)

式(6)から分かるように固有値の次元は電力であり、物理量としては平均受信電力である。固有ベクトルはこの受信電力を得るためのアレー合成ウェイトに対応する。

信号部分空間と雑音部分空間
信号部分空間 ES：固有ベクトル e1～eL の L 個のベクトルで張る空間
雑音部分空間 EN：固有ベクトル eL+1～eM の M-L 個のベクトルで張る空間

ES = ( e1 e2

EN = ( eL+1 eL+2

eL )

(7a)

eM )

(7b)

方向行列 A は信号部分空間にあり、構成要素である方向ベクトル al（l=1, 2, …. L）は固有ベ
クトル e1～eL の線形結合で表される。

5

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

４．MUSIC 法
発想
たとえ話から入りたい。ある村を突風が襲った。どの家も甚大な被害を受けたが、一軒だけび
くともしない家があった。その不思議を探ろうとこの家を調査してみたところ、一つの方向に極めて頑丈な壁があることが分かった。そうか、風はこの方向から吹き付けていたのだと知ることができた。これから学ぶ MUSIC 法も同じである。
図１で示した M 素子のアレーアンテナ（ULA）で受信する場合を考える。前章で学んだ相関
行列 Rxx は、式(4)より M 個の固有値とそれぞれに属する M 個の固有ベクトルに分解できる。
普通、電波の到来方向を調べようと思うときは、信号部分空間にある L 個の固有値を調べに行きたいと思うであろう。MUSIC 法は逆に、電波の電力を含まない M-L 個の固有値を生み出す
雑音部分空間に着目するのである。固有値が電波の影響を受けないと言うことは、電波が来ているのに、何も受信しないということ。電波の影響を受けない M-L 個の固有ベクトルが作り出す
アンテナパターンを調べようと言うことである。この M-L 個の固有ベクトルから M-L 種のアン
テナパターンが描ける。それを全部重ねてみる。実際に電波が到来している方向はその全てでヌル点になっているはずである。アンテナパターンのヌル点は鋭い落ち込みになっている。このアンテナパターンを上下反転させれば、その落ち込みがピークになり、電波の到来方向を指していることになる。このようにして、電波の到来方向が同定できる。受信強度の最大方向を求めるのとは比べ物にならない精度になる。頭の丸いピーク方向に答えを求めるのではなく、楔形のヌル点方向を見ると言う発想である。図３は上述の MUSIC 法の原理をまとめている。

方法
L+1 番目から M 番目の固有値L+1～M は PN（測定値では多少のばらつきがあり全く同じにはならない）となり、それらの固有ベクトル eL+1～eM で与えられるアンテナ合成ウェイト wL+1～wM
は、入射波を受け付けないアンテナパターンを作る。これらの固有ベクトルが張る M-L 次の空
間は雑音部分空間である。前述のとおり、その M-L 個のアンテナパターンを一つ一つ描いてみ
れば、電波到来方向については、全てでヌルが形成されていることになる。ただし、その一つ一
つには、到来方向以外にも偶然にヌルができている可能性がある。このため、共通にヌルをもつ
全部の電力パターンを重ね合わせることにより、真の到来方向のみにヌルを作ることができる。
このアンテナパターンは以下の式で与えられる。
M

G( )   emH a( )2(8)m = L +1

到来方向をピークで指し示すには、上記アンテナパターンのヌル点方向がピークになるようにすればよく、これはすなわち、アンテナパターンを上下反転させればよい。最終的に求めたい
到来方向推定は、aH()a()を掛けて正規化し、以下の計算式による。

6

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

図３

P( ) =a H ( )a( )M e a( )m = L +1=

MUSIC 法早分かり

Hm

2

a H ( )a( )a H ( ) EN ENH a( )(9)

この方式の欠点は、到来方向が数値で求まらず、式(9)の計算によってピークサーチをしなければならない点である。また、ピーク角度はあくまで方向を指し示すものであって、入射波の強度の情報が得られないことである。ただし後者については、以下の方法によって信号強度を同定
7

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

することができる[1], [3]。
入射波の数 L と到来角度1～L が、上記の方法によって定まったとする。これより、式(1)の方
向行列 A が定まる。各到来波の強度は測定値から求められる相関行列と雑音電力を用いてQ = ( AH A)−1 AH ( Rxx − PN I M ) A( AH A)−1(10)を求める。この行列 Q の対角成分に各到来波の受信電力が現れる。さらに、アンテナ出力ベク
トル x(t)と方向行列 A を用いて、入射信号 s(t)を最尤推定法により分離できる。

s (t ) = ( AH A)−1 AH x(t )(11)

式(11)で、x(t)の係数である(AHA)-1AH は信号を分離するためのウェイト係数行列 WH であり、これを(w1, w2, …, wL)H とすれば、それぞれのウェイトが L 個の到来波を選択受信するアレーウ
ェイトになる。
（この WH は行列 A の一般逆行列（Moore-Penrose の逆行列）である。
）
推定結果の一例（[5]の再掲）
半波長間隔で直線状に配置された６素子アレー（M=6）を用いる。これに、到来角度①50°, ②120°,③ 130°から３波が入射する。信号は±１の乱数系列で、それぞれの振幅を 1 (0dB), 2 (6dB),
4 (12dB)とする。素子ごとの平均雑音電力 PN は-10 dB とする。SN 比に換算すると、上記の３波
は、10dB, 16dB, 22dB である。
MUSIC 法での信号処理は、受信アンテナの出力信号 x(t)を得てから後のことであるが、その信号つくりのために、式(1)で定める方向行列 A を用いて x(t)の作成からはじめる。実際には、スナップショット的にデータを作るので x(i) (i=1, 2, …, N)とする。
スナップショット数 N をパラメータにして<x(i)xH(i)>で相関行列 Rxx を作る。次に、Rxx の固有値と固有ベクトルを求める。このケースでは、固有値は N=100 程度でほぼ安定した値になり、
110, 11.2, 4.9, 0.13, 0.10, 0.09 が得られた。第３固有値までが、雑音電力（0.1）に比べて十分大きな値になっており、到来波数 L は 3 であると判定できる。そこで、雑音部分空間を張る固有ベク
トル e4, e5, e6 を用いて、式(9)により MUSIC スペクトルを求める。図４がその結果である。(a)はスナップショット数が十分でない代表として N=10 を、(b)は十分な数の例として N=1000 の場合を示している。両図とも、与えたとおりの到来方向（①50°, ②120°,③ 130°）を指し示しているが、スナップショット数が足りないと、②の推定に見られるように分解能が甘くなっている
様子がわかる。要は、スナップショット数が少ないと、入力データの部分的な相関や雑音の影響
で、平均値として求められる相関行列 Rxx が安定しないためである。このため、実際の信号処理
では、スナップショットの数を大きめにとるのが良いであろう。N=100 と 1000 の場合に、正しく推定された到来角度を用いて、(10)式で入力信号の電力を求めると、①16.07 (N=100), 16.02(N=1000)、②3.89 (100), 3.97 (1000)、③0.99 (100), 1.00 (1000)となり、設定電力が①16, ②4, ③1 であったので、N=100 程度で十分な推定ができることがわかる。図５は、(11)式により③の到来波
8

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

を選択受信するウェイトでアンテナパターンを書いたものである。図より、①と②を受信せず、
かつ、最尤法によりアレーアンテナの残りの自由度を使って最大利得で③波を受信するパターンが実現できていることがわかる。
ここでは、入力データをシミュレーションによって作っているので、誤差要因がない理想的なデータを使ったことになる。実際には、アレーアンテナの素子の特性の不ぞろいやアンテナ間のカップリングなどにより MUSIC スペクトルには多かれ少なかれシミュレーションには含まれない誤差が生じる。そのため、測定に際しては、既知の方向の電波を受信する等によるキャリブレーションが重要になる[4]。

(a)
図４

図５

N=10

(b)

N=1000MUSIC 法に到来方向推定（６素子半波長間隔アレーの場合）

信号選択のためのアンテナパターン（①、②を除去し③を選択）

9

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

５．ESPRIT 法
ESPRIT は干渉計の原理を、到来波数とアレー素子数が多い場合に拡張したものである。
ESPRIT の原著論文[2]やそれを解説する専門書[3], [4]では、特性が同じ２本のアンテナ（ダブレット(doublet)と呼ばれる）を一対にし、これを複数配置した二つのサブアレー（図６）で原理説明が行われている。ダブレットの対でアンテナパターンが同じであれば、ダブレット毎のパターンは自由でよい。この場合は、(ダブレットの数-1)>L の到来方向推定が可能である。MUSIC 法
では L<M-1 であるので、必要なアンテナ素子数が２倍になる。しかしながら、図１に示す等間
隔線形アレー（ULA）の構成にすると、MUSIC 法と同様、M-1 以内の到来方向を推定することができ、この方が実際的な構成でもある。そこで、本稿では ULA の構成について原理説明を行う。

図６ ESPRIT が想定するアレーアンテナ系（複数のダブレットとそれらで構成する二つのサブアレー）（ただし、本稿での原理説明では図７の ULA を用いる）

干渉計の原理干渉計の原理をおさらいするために、以下の二つのケースについて、干渉計原理で到来方向を求めてみよう。
(i) M=2, L=1
受信信号

x1 (t ) = s1 (t ) + n1 (t )x2 (t ) = e− jk d cos1 s1 (t ) + n2 (t )10

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

相関演算

x1 (t ) x2* (t ) = e jk d cos1 PS1
測定値<x1x2*>から位相差 kdcos1 が得られ、到来方向1 が容易に求められる。ただし、1 の角度範囲によって、位相差が±を超えないように素子間距離d を設定しないと、角度に不確定性が生じる。1 が 0～ の範囲にあるならば、d <0.5 波長が条件になる。
(ii) M=3, L=2
受信信号

x1 (t ) = s1 (t ) + s2 (t ) + n1 (t )x2 (t ) = e− jk d cos1 s1 (t ) + e− jk d cos2 s2 (t ) + n2 (t )x3 (t ) = e− j 2k d cos1 s1 (t ) + e− j 2k d cos2 s2 (t ) + n3 (t )
相関演算x1 (t ) x1* (t ) = PS1 + PS 2 + PN = x2 (t ) x2* (t ) = x3 (t ) x3* (t )x1 (t ) x2* (t ) = e jk d cos1 PS1 + e jk d cos2 PS 2 = x2 (t ) x3* (t )x1 (t ) x3* (t ) = e j 2k d cos1 PS1 + e j 2k d cos2 PS 2
相関行列 Rxx を構成する９個の要素のうち、測定値で独立なものは、相関演算の左辺の３つの量。一方、右辺では PS1, PS2, PN, 1, 2 の５つの量が未知のものになる。到来方向推定では、1, 2が求められれば良く、他のものを求める必要はないが、３つの連立方程式から1, 2 が求められるのであろうか？これを直接に求めようとすると、結構大変そうに見える（*）
。
【注*：一般化したケースについて ESPRIT は答えを見つけるのであるから、当然、このケースでも解はある。直接的には、ESPRIT で得られた一般解で M=3, L=2 とすればよいはず。それを書き下してまとめれば、この場合の答えになる。その確認は計算が得意な人にお任せしたい。
】。

ESPRIT
これを一般化してゆくのは、簡単ではなさそうと言うことが理解できたと思うので、以下、
ESPRIT による具体的な解法を見てゆきたい。
図１の構成において、図７のように二つの M-1 素子アレーがd だけずれたサブアレー１，２
を考える。
11

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

図７間隔d の M 素子 ULA と２つのサブアレー

サブアレー１，２の方向行列 A1, A2 は、(1)式の方向行列 A を、次式のように共有する。

(11)
したがって、A1 と A2 の関係は

A2 = A1 e − jk d cos1

0
 =


0


(12)

0
e

− jk d cos 2

0






e − jk d cos L 
00

ここで、対角行列  は回転因子を要素にしており ESPRIT の名前の由来である Rotational
Invariance （回転不変性）は、(12)式の関係を指している。
A1 と A2 が測定で求められれば、対角行列 が定まり、それより、目的とする到来角度l が求
まるが、A1, A2 は未知であるので、行列 を求めることができない。
そこで、相関行列の信号部分空間に目を向けてみる。ここでも、信号部分空間 ES から以下のように二つの部分空間 ES1, ES2 を作る。

12

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa(13)
先に述べたように、方向行列 A の空間と信号部分空間 ES は同じ空間内にあるので、L×L の正則行列 T を用いて

A = ES

(14)と A と ES を結びつけることができる。同様に A1 と ES1、A2 と ES2 も、式(11), (12), (14)より、以下で結び付けられる。

A1 = ES1A2 = ES 2 = A1
これを整理し、ES 2 = A2 −1 = A1 −1 = ES1 −1(15)となる。ここで L×L 次元の行列 TT-1 を =  −1(16)とおくと、ES1 と ES2 は次式で結び付けられる。

ES 2 = ES1

(17)

ES1 と ES2 は相関行列から得られる量であり、工夫をすれば、 は求められる量になる。かりに 求まったとしても、まだ T が未知であり、目的の にはたどり着いていないが、網は絞り込まれてきている。
ところで、行列の対角化において「正則行列 に対して、対角行列 と正則行列 T を用いて、T-1T= と対角化ができれば、言い換えると、 =T T-1 で表されるとき、 の対角要素は、の固有値である」と言う性質がある。故に、行列 の固有値：1’, 2’, …, L’が求めたかった対角行列 の要素であったと言うことになる。これより、

13

MUSIC と ESPRITl  = e− jk d cosTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

(18a)

l

→

 arg   l k d l = − cos−1 (18b)として、到来角度を定めることができる。
ESPRIT はこのようにして、相関行列から計算のみにより到来方向l を推定することができる。
これで、原理の説明は終わるのであるが、途中であいまいにしてきたことが残っている。式(17)から行列 を求める演算である。行列 を求めるのに、式(17)を解けばよいが、ES1 は正方行列ではないので、単純に逆行列を求めて について解くことはできない。これを求める方法として、二つの方法が提案されている。LS (Least-Squares)-ESPRIT と TLS (Total- Least-Squares)-ESPRITである。二つを比べると LS の方が簡易であり、精度は TLS の方が高い。
LS-ESPRIT
ES1, ES2 は共に(M-1)×L の行列であり、未知数よりも方程式の数の方が多い。このような場合は|| ES2- ES1 ||（||・||はフロベニウスノルム）の最小化を行う最小二乗法が用いられ、結果のみ示すと次式になる。右辺 ES2 にかかる係数部分 ( ESH1ES1 )−1 ESH1 は ES1 の一般逆行列（Moore-Penrose の逆行列）である（同じ形の式は(11)式にも）。

 = ( ESH1ES1 ) ESH1ES 2
−1

(19)

TLS-ESPRITLS-ESPRIT で規範とした|| ES2- ES1 ||の に関する最小化では、最小二乗法の原理により ES2に含まれる誤差の最小化である。誤差は ES1 にも ES2 にも同等に有るので、両者の全体としての誤差の最小化から を求める方法が TLS-ESPRIT である。導出の詳細と結果式については[2]あるいは[3], [4]をみてほしい。
（実用的には LS-ESPRIT の適用で十分の場合が多い(p. 15 の表１)）

ESPRIT のキーポイント
ESPRIT 法は、干渉計の原理から導かれる A2 = A1 の中の  に、求めたい情報が含まれるが、方向行列 A が未知なので、直接これを解くことはできない。そこで、A と同じ部分空間を構成する信号部分空間に着目し、 A2 = A1 の関係をここに求めると、 ES 2 = ES1 と書き表され、このは測定値から得られる ES1 と ES2 を用いて解ける。そして、求めたかった  は を対角化したもの（＝ の固有値を要素とする対角行列）であった、と言うストーリーである。このようなス
トーリーにまとめ上げた Roy & Kailath の慧眼にたいして脱帽するしかない。

推定結果の一例
ESPRIT 法による推定結果の一例を示す。環境設定は、４節の MUSIC 法で示した例（図４）と同じである。すなわち、半波長間隔で直線状に配置された６素子アレー（M=6）を用いる。これに、到来角度①50°, ②120°,③ 130°から３波が入射する。信号は±１の乱数系列で、それ14

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

ぞれの振幅を 1 (0dB), 2 (6dB), 4 (12dB)とする。素子ごとの平均雑音電力 PN は-10 dB とする。SN比に換算すると、上記の３波は、10dB, 16dB, 22dB である。(17)式で与えられる行列 を求める
に際しては(19)式の LS-ESPRIT を用いた。
シミュレーションでは、スナップショット数 N を 10, 30, 100, 300, 1000 の５種類とした。表１
は設定値（）に対する推定結果の一試行例（雑音の乱数を替えて多数回試行した後の平均値ではない）を示している。表より、この程度の SN 比であれば、データサンプル数 N がずいぶん少ないと思われる 10 程度であっても、誤差 1°以内の精度で、また、N が 100 以上では、ほぼ 0.1°以内の精度で推定できることが分かる。また、行列 を求めるに際して、簡易な計算式である
LS-ESPRIT（(19)式）でまったく問題ないことも分かる。なお、行列 の固有値l’を求めると、その絶対値が 0.99～1.01 程度にばらつき、求めたい exp (-jkd cos l)の 1.0 とは完全に一致しない（すなわち、(18a)式が厳密に満たされることは無い）。その場合でも、arg(l’)として、(18b)式
中の位相を求めることには、この不一致はほとんど影響していないように見える。
求めたい到来角度が一連の手順でずばり計算できる ESPRIT は、有効性が非常に高いと実感できるであろう。手法を理解するには、それなりの時間がかかるであろうが、step-by-step で計算式を書き下してみると極めてシンプルで、行列計算を得意とする計算ソフト（Matlab など）に非常に相性がよい。
表１ ESPRIT でのシミュレーション結果の一例（M=6, L=3, d=0.5 波長；表内の数値の単位は[ °]）６．両手法のまとめ
二つの手法：MUSIC 法と ESPRIT 法について、信号処理として共通する前処理と、それぞれの個性が発揮される本処理に分けて整理したものを表２に示す。
以上見てきたとおり、高精度到来各推定として、MUSIC 法と ESPRIT はディジタル信号処理の傑作である。MUSIC 法が最後の段階において、角度空間をサーチして到来方向を特定するのに対して、ESPRIT は計算だけで求められることが特徴になる。ESPRIT において、図６のようなダブレットを用いる測定では、ダブレット間の校正は不要と言うのも特徴であるが、図１のような ULA 構成による測定では、両手法共に、事前のアレー校正は必須である[4]。

15

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

表２

MUSIC 法と ESPRI 法の特徴と手順の比較

ここでは、送信源に相関がないことを前提に述べたが、完全相関でなければ本稿でまとめた手順がそのまま通用する。マルチパスのように、同一の送信源からの到来波の場合（いわゆるコヒーレント波の入射）に対しては、空間平均法と言った前処理が必要になる[3]。
MUSIC 法では、到来波を平面波として角度空間でサーチしたが、近距離等で素子アンテナ毎に入射角が違う場合でも、波源位置を含む面でサーチすれば、点推定も可能である。
本稿では、高分解能到来方向推定法のレジェンドである MUSIC と ESPRIT について、発想の源と基本原理を解説した。その後、両手法とも、種々の改良・高機能化も進み、さらには、新たな手法も提案されて、今日の無線システムの中の様々な場面で活躍している[3], [4]。

余談
このレポートを書いている今日この時点（R02.07.31）において、コロナウィルス感染の第２波が日本全体に押し寄せ、世の中、大変な不安に掻き立てられている。以下、たとえ話。感染が深刻なある都市に、二人の医師、MUSIC 博士と ESPRIT 博士が乗り込み、対策に当たることにした。ESPRIT 博士は、感染した人から症状や経緯を聞いて対策の構想を練った。MUSIC 博士は、感染しなかった人だけを集めて、なぜかからなかったかを綿密に調べた。それぞれのデータに基づき、二人は、まったく異なる対策を提案して実行したが、その両方共に抜群の効果が認められて、その後の感染を防ぐことができた。感染者・非感染者両方の人の話を聞いて総合的に考えれば、もっとよい対策があるのではと思うところもあるが、対策自体は、目的としてはもうそれで十分なようである。
16

MUSIC と ESPRITTechnical Report YK-049-rev
August 18, 2020
Y. Karasawa

参考文献[1] R.O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas and Propagat.,vol. AP-34, no. 3, pp. 276-280, 1986.
[2] R. Roy and T. Kailath, “ESPRIT- Estimation of signal parameters via rotational invariance techniques,” IEEE Trans.
Accous., Speech, Signal Proc., vol. 37, no.7, pp. 984-995, 1989.
[3] 菊間信良, アレーアンテナによる適応信号処理, 科学技術出版, 1998.
[4] 山田寛樹, 高分解到来波推定の基礎と実際, 信学会, AP研, アンテナ・伝搬における設計・解析ワークシ
ョップ, No. 33, pp. 1-133, 2006.10.
[5] 唐沢好男, “スペースダイバーシチ・アダプティブアレー・到来方向推定法:アレーアンテナによる空間信号処理の三兄弟,” 私報：TR-YK-006, pp. 1-21, 2018. 01,
http://www.radio3.ee.uec.ac.jp/ronbun/TR_YK_006_Adapytive_Array.pdf

読者の皆さんへ
このレポートには続編（TR-YK-050: F/B型空間平均法）があります。
さらに、TR-YK-049とTR-YK-050を一つにまとめて整理した資料：TR-YK-070:
Open Access Book 2: 固有値が開く無線技術の扉（第５章 電波到来方向の高分解能推定）を公開しています。

17