IUGONET研究会2019
September 13, 2019京都大学宇治キャンパス

短時間データに含まれる複数信号の周波数推定手法ESPRIT (Estimation ofSignal Parameters via RotationalInvariance Techniques)の紹介
山本衛 (京都大学生存圏研究所)

ESPRIT法とは• ある信号を等間隔に配置された複数の素子で観測するとき、素子間の位相差を推定する手法。
• アレイアンテナの電波到来方向の推定手法として開発されたらしい。しかしながら、• 電波到来方向の推定＝複数アンテナ間の信号位相差• 周波数推定＝時間サンプル間の信号位相差（の微分）• 距離推定＝複数周波数で得られたレーダーエコーの位相差から求める（例えばFM-CWレーダー）• 全てに適用できる。特長は、• 少数サンプル（短時間のデータ）から推定できる。
• 複数の信号を分離して推定できる。
• フーリエ変換に基づく方法のように、周波数が信号長の整数分割に縛られることがない。

• 問題点は、• 等間隔サンプリングが必要• （有効な）信号の数を決めてから解析する。
• 信号強度は求まらない（求めにくい？）

2

ESPRIT法（問題設定）• アレイアンテナの電波到来方向の推定として説明• M本のアンテナが一列に並んでいる。（1次元アレイ）• 各アンテナで得られるデータ: , ,… ( )• 雑音: , ,… ( )• Ｌ個の方向（角度 , ,… ）から電波が到来する。
求めたいものはこれ。

3

ESPRIT法（問題設定）• 信号 ( )をベクトル表記すると

=, ,…
=+
ここで、
=( ), ( ), …
( ) は、L個の方向から来る電波の信号。さらに、求めたいもの• =,, … ( ) （方向行列）•• ( )=•=( ),,( ), …

( ), …

( )( )(方向ベクトル）(雑音ベクトル）

4

ESPRIT法（考え方）• アンテナ列を2分割し、その間の位相差を求める。
• サブアレイ１：アンテナ1番～M-1番。。
• サブアレイ２：アンテナ2番～Ｍ番。。
• 関係性は、 Ψ =
ここで、
• =|,
•

=

• Ψ = diag, |,|,…,

5

ESPRIT法（考え方）
M
M-1J1

Ⅼ

MMⅬⅬ

Ψ

＝M-1

J2A

ⅬM

A

• アンテナ列を2分割し、その間の位相差を求める。
• サブアレイ１：アンテナ1番～M-1番。。
• サブアレイ２：アンテナ2番～Ｍ番。。
• 関係性は、• =|•=• Ψ = diag, |Ψ=,,|

ここで、,…,

この式をΨについて解きたい。
しかし、A（方向行列）もまた未知なので、このままでは解けない。
6

信号の相関行列と固有値分解• アンテナ（M本）から得られる信号の相関行列は、≜ E ( ) ( ) 、 ここでEは期待値を求める操作。実際には時間平均または空間平均を行う。 は共役転置（以後同じ）。

• 相関行列の固有値を求め、固有値を大きい方から並べて固有値分解する。（ 注：はエルミート行列）= ΛM

M RxxM

M＝ MEMΛMMEH• 固有ベクトルを並べた行列と固有値は、信号部分空間と雑音部分空間に区別できる（次ページ）

7

信号部分空間と雑音部分空間
M
LM M 信号成分 Es

M
M-L

雑音成分
En

λ1

λ2

M…

0λL
λn0…
λn固有ベクトル行列Esの構造

固有値の対角行列Λの構造
λ1 > λ2 > … > λL > λn > ….最後まで
信号部分

雑音部分

8

ESPRIT法（解くべき方程式）• （天下りですが）信号部分空間の固有値ベクトル行列Esと、方向行列Ａの関係=を満たす正則行列Tがある。
• これを用いると解くべき方程式は以下の通り。
Φ=
Φ=
ΨL

J1 Es ＝ M-1
M Es

L

J2 Es ＝
各列の1行目からM
M-1行目までを取る

Es
M-1

各列の2行目からM行目までを取る• 上の方程式を最小二乗法（TLS (Total-Least-Square)法で解いてΦを求め、さらにその固有値として、素子間の位相情報であるΨが求められる（少しややこしい行列演算を行うため、ここでは省略）。
9相関行列の計算
時系列信号から瞬時周波数の推定を目指す。
信号X(t) = X1, X2, X3, … XN (N個の等間隔の複素時系列信号とします）
最大ラグをM (=<N)として、自己相関関数Rxxを求める。
M
E[X1X1*] E[X1X2*] …
E[X2X1*] E[X2X2*] …
Rxx = M …
……………………………… E[XMXM*]

RxxはM x Mサイズの複素数のエルミート行列になる。上記のE[ ]はアンサン
ブル平均、この場合は、M < N なので、自己相関関数の自然な計算で平均化できる。
Rxxを固有展開する。固有値λを求めて大 小とならべて固有値の対角行列Λをつくる。対応する固有ベクトルを並べて固有ベクトル行列Eを作る。
固有値の大きいL個が信号(s)成分、それ以外はノイズ(n)成分である。
10

解析例
• 条件

• 複素信号、サンプリング間隔＝1ms• 雑音：ガウシアン乱数（rms = 1）
• 3つの正弦波信号を重畳する。
• 周波数 1Hz、振幅 = 5• 周波数 -13Hz 、振幅 = 1• 周波数 399Hz 、振幅 = 1• 調べ方• 3つの正弦波が到来していると仮定する（L=3）• 最大ラグ数＝30点（固定）• データ長を1024点（約1秒）、128点（約1.3秒）、64点（0.64秒）と変化させて、時系列データ、FFTパワースペクトル、ESPRIT法の結果を示す。
11

信号点数1024（約1秒間）

Esprit法の結果＝20回試行の平均値（標準偏差）
1.002 Hz （0.017 Hz）-13.15 Hz （0.374 Hz）
399.01 Hz （0.117 Hz）

12信号点数128（約1.3秒間）

Esprit法の結果＝20回試行の平均値（標準偏差）1.214 Hz （0.123 Hz）-17.75 Hz （2.667 Hz）
399.11 Hz （0.349 Hz）

13信号点数64（約0.64秒間）

Esprit法の結果＝20回試行の平均値（標準偏差）
1.991 Hz （0.330 Hz）-20.29 Hz （6.214 Hz）
399.16 Hz （0.822 Hz）

14

1st COSMIC-2 Beacon Signal (1)September 6, 2019 11:23UT-, Bangkok, Thailand
965.6 MHz

400.9 MHz

Monitoring 400.9 MHz and 965.6 MHz signal from FORMOSAT-7/COSMIC-2• Started on September 1, 2019
• Bangkok (Thailand), BKK, 13.7278 deg N, 100.7724 deg E• Ho Chi Minh (Vietnam), HCM, 10.8486 deg N, 106.5598 deg E• Kototabang (EAR site, Indonesia), EAR, 0.2044 deg S, 100.3203 deg E1st COSMIC-2 Beacon Signal (2)September 6, 2019 11:23UT-, Bangkok, Thailand

Signal processing scheme• Frequency finding and down samplingwith preserving phase (difference)
information. Three-step processing.
• Sampling speed decreases fromoriginal 100 kHz (last page)(top left) 1 kHz(top right) 50 Hz … for S4 index(bottom right) 10 Hz … for TEC

ESPRIT法とは• ある信号を等間隔に配置された複数の素子で観測するとき、素子間の位相差を推定する手法。
• アレイアンテナの電波到来方向の推定手法として開発されたらしい。しかしながら、• 電波到来方向の推定＝複数アンテナ間の信号位相差• 周波数推定＝時間サンプル間の信号位相差（の微分）• 距離推定＝複数周波数で得られたレーダーエコーの位相差から求める（例えばFM-CWレーダー）• 全てに適用できる。特長は、• 少数サンプル（短時間のデータ）から推定できる。
• 複数の信号を分離して推定できる。
• フーリエ変換に基づく方法のように、周波数が信号長の整数分割に縛られることがない。

• 問題点は、• 等間隔サンプリングが必要• （有効な）信号の数を決めてから解析する。
• 信号強度は求まらない（求めにくい？）

17