MUSIC 法による複数正弦波のパラメタ推定と Python 実装

はじめに

本記事は最近のマイブームとなっている「正弦波のパラメタ推定シリーズ」のうちの一つである。

今回は、高分解能な周波数推定手法の一つである MUSIC (MUltiple SIgnal Classification) 法 [1-5] について、2つのバリエーションである spectral-MUSIC と root-MUSIC を Python で実装した。両者は、信号の共分散行列から推定した雑音部分空間を利用するという共通の原理に基づいている。さらに周波数だけでなく、最小二乗法に基づいて振幅と位相も推定する方法についても実装した。

本記事では、それぞれの理論的背景を概説した後、検証実験を通じてその性能と特性を比較した結果を報告する。

はじめに
MUSIC 法の共通原理：部分空間の直交性
周波数推定の2つのアプローチ
- spectral-MUSIC
- root-MUSIC
最小二乗法による振幅・位相の推定
実装
実験的評価
考察
おわりに
参考文献

MUSIC 法の共通原理：部分空間の直交性

信号モデル

信号モデルとして、 $M$ 個の複素正弦波と白色雑音 $w[n]$ の和で構成される離散信号 $x[n]$ を考える。 $F_s$ は標本化周波数である。

$\begin{align*} x[n] = \sum_{k=1}^{M} A_k \exp\left\{j\left( \frac{2\pi f_k n}{F_s} + \phi_k\right)\right\} + w[n] \end{align*}$

$j$ は虚数単位、 $f_k$ は物理的な周波数 [Hz] である。この信号モデルは、正規化角周波数 $\omega_k = 2\pi f_k / F_s$ を用いて、より簡潔に表現できる。

$\begin{align*} x[n] = \sum_{k=1}^{M} c_k \exp (j \omega_k n) + w[n] \end{align*}$

$\omega_k$ は「1サンプルあたりの位相変化量」[rad / sample] であり、信号の周期性を特徴づける重要な量である。 $c_k =A_k \exp(\phi_k)$ は複素振幅である。

共分散行列の構築

信号を以下のような $L$ 次元のベクトル群に変換する。

$\mathbf{v}_0 = [x[0], x[1], ..., x[L-1]]^{\top}$
$\mathbf{v}_1 = [x[1], x[2], ..., x[L]]^{\top}$
$\mathbf{v}_2 = [x[2], x[3], ..., x[L+1]]^{\top}$
...

これら $L$ 次元のベクトル $\mathbf{v}_{i}$ は、窓長 $L$ の矩形窓で信号を切り出した $L$ サンプル分の断片たちである。アレイ信号処理分野の用語を借りて、これらを「スナップショット」（あるいはスナップショットベクトル）と呼ぶことにする。窓長 $L$ は後述する部分空間の次元 $L$ に対応しており、重要なパラメタである。

共分散行列 $\mathbf{R}_x$ は、これらのスナップショットベクトル $\mathbf{v}_{i}$ たちが、平均的にどのような相関を持っているかを表す行列である（サイズは $L \times L$ ）である。

$\begin{align*} \mathbf{R}_x = \mathbb{E}[\mathbf{v}_{i} \mathbf{v}_{i}^{*} ] \end{align*}$

ここで $\mathbf{v}_{i}^{*}$ は $\mathbf{v}_{i}$ の随伴（共役転置）である。定義は無限個のスナップショットに渡る期待値であり、したがって無限長の信号を必要とするが、実際には無限長の信号は得られない。あくまで観測された信号は有限長であり、そこから得られるのは有限個のスナップショット（ $N$ 個とする）のみである。そこで、期待値を標本平均で近似（推定）する。

$\begin{align*} \mathbf{R}_x \approx \widehat{\mathbf{R}}_x \triangleq \frac{1}{N} \sum_{i=0}^{N-1} \mathbf{v}_{i} \mathbf{v}_{i}^{*} \end{align*}$

ここで $\widehat{\mathbf{R}}_x$ は標本共分散行列である。さらに、標本平均の計算を効率的に行うために、Hankel 行列が用いられる。全てのスナップショットベクトル $\mathbf{v}_{i}$ を列として並べたデータ行列 $\mathbf{X}$ を考えると、これが Hankel 行列になっている：

$\begin{align*} \mathbf{X} &= [\mathbf{v}_{0}, \mathbf{v}_{1}, \mathbf{v}_{2}, \ldots, \mathbf{v}_{N-1} ] \\ & = \begin{bmatrix} x[0] & x[1] & x[2] & \cdots & x[N-1] \\ x[1] & x[2] & x[3] & \cdots & x[N] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x[L-1] & x[L] & x[L+1] & \cdots & x[N + L -2] \\ \end{bmatrix} \end{align*}$

ただし、正方行列ではない Hankel 行列の定義を採用した。こうすると scipy の hankel 関数を使うときに都合がよい（Link）。この Hankel 行列 $\mathbf{X}$ を使うと、標本共分散行列は次式に示す簡潔な形で計算できる。

$\begin{align*} \widehat{\mathbf{R}}_x = \frac{1}{N} \mathbf{X} \mathbf{X}^{*} \end{align*}$

実装上、for文を使わずに行列積のみで計算するので、高速化が期待できるというわけである。

補足：前方後方平均化による共分散行列の推定精度向上

MUSIC 法の推定精度は、共分散行列の推定精度に大きく依存する。信号の観測サンプル数が限られている場合、この共分散行列の推定誤差が性能低下の主な原因となる。この問題に対処するための効果的な手法が、前方後方平均化 (Forward-Backward Averaging) [5] である。

正弦波信号は、時間反転させても、その統計的性質（具体的には自己相関）が変わらないという性質を持つ。

前方 (Forward) データ: 元の信号 $x[n]$
後方 (Backward) データ: 元の信号を時間反転し、複素共役を取った信号 $x_b[n] = \overline{x[N-1-n]} = x[N-1-n]^{*}$

前方後方平均化は、これら2つの等価な信号からそれぞれ共分散行列を計算し、それらを平均する。これにより、少ないサンプル数から、よりロバストな共分散行列の推定値を得ようという技術である。

これまで計算してきた通常の共分散行列を前方共分散行列 $\mathbf{R}_f$ と呼ぶ。

$\begin{align*} \mathbf{R}_f = \mathbb{E} [x x^{*}] \approx \widehat{\mathbf{R}}_f \end{align*}$

続いて、時間反転・複素共役した信号 $x_b$ から後方共分散行列 $\mathbf{R}_{b}$ を定義する。

$\begin{align*} \mathbf{R}_b \triangleq \mathbb{E} [x_b x_b^{*}] \approx \widehat{\mathbf{R}}_b \end{align*}$

$\mathbf{R}_b$ は元の前方共分散行列 $\mathbf{R}_f$ を用いて、

$\begin{align*} \mathbf{R}_b = \mathbf{J} \overline{\mathbf{R}}_{f} \mathbf{J} \end{align*}$

とも書ける（標本共分散行列も同様）。 $\mathbf{J}$ は対角線がすべて 1 で、他の要素が 0 の行列（単位行列を左右反転させたもの）である。 $\overline{\mathbf{R}}_{f}$ は $\mathbf{R}_{f}$ の各成分について複素共役を取った行列である。

最後に、これら2つの行列を単純に平均して、 前方後方平均化共分散行列 $\mathbf{R}_{fb}$ の推定値 $\widehat{\mathbf{R}}_{fb}$ を計算する。

$\begin{align*} \mathbf{R}_{fb} & \triangleq \frac{\mathbf{R}_f + \mathbf{R}_b}{2} = \frac{\mathbf{R}_f + \mathbf{J} \overline{\mathbf{R}_{f}} \mathbf{J}}{2} \\ & \approx \frac{\widehat{\mathbf{R}}_f + \mathbf{J} \widehat{\overline{\mathbf{R}}}_{f}\mathbf{J}}{2} \\ & \triangleq \widehat{\mathbf{R}}_{fb} \end{align*}$

この $\widehat{\mathbf{R}}_{fb}$ は、 $\widehat{\mathbf{R}}_f$ 単独よりも統計的に安定した推定値となり、結果として MUSIC 法の推定精度を向上させる。特に短時間の信号や低 SNR の信号に対する推定性能が改善されることが期待できる。 $\widehat{\mathbf{R}}_{fb}$ がより「統計的に安定した推定値となる」ことの詳しい説明は、記事を改めて行いたい。

固有値分解と部分空間分離

共分散行列 $\mathbf{R}_x$ を固有値分解する。

$\begin{align*} \mathbf{R}_x = \mathbf{E} \boldsymbol{\Lambda} \mathbf{E}^{*} \end{align*}$

ここで、 $\mathbf{E}$ は固有ベクトルを列に持つ行列、 $\boldsymbol{\Lambda}$ は固有値を対角成分に持つ対角行列であり、 $\mathbf{E}^{*}$ は $\mathbf{E}$ の随伴（共役転置）である。固有ベクトルは、信号成分がどのような波形の組み合わせで表現されるかという、いわば基本パタンの集まりを表していると見なせる。また各固有値（の絶対値）は、対応する固有ベクトルの、信号全体のパワーに対する寄与率を表していると見なせる。

信号成分は特定の周波数を持つ波形であり、そのパワーは特定の固有ベクトル方向に集中するはずである。雑音成分はあらゆる周波数成分を均等に含む、完全にランダムな波形である。そのパワーは、特定の方向に集中せず、全ての固有ベクトル方向に均等に分散するはずである。

それゆえ、 $\mathbf{R}_x$ の固有値（ $L$ 個存在）はまた、2つのグループに分かれる。 $M$ 個の正弦波に対応する $2M$ 個の固有値※ は「信号のパワー＋雑音のパワー」を反映しており、大きな値を持つ。残りの $L - 2M$ 個の固有値は「雑音のパワーのみ」を反映しており、ほぼ全て小さな値（雑音の分散）のオーダーになる。

この固有値の値に従って、固有ベクトル全体もまた2つのグループに分かれる。

信号部分空間 $\mathbf{E}_s$ : $2M$ 個の大きな固有値に対応する固有ベクトルが張る部分空間
雑音部分空間 $\mathbf{E}_n$ : $L - 2M$ 個の小さな固有値に対応する固有ベクトルが張る部分空間

これらの部分空間が（近似的に）互いに直交している点が重要である（後述）。

ちなみに、 $\mathbf{E}_s$ は $2M$ 本の固有ベクトル（ $L$ 次元）が並んだ行列として保持する。ゆえに $\mathbf{E}_s$ のサイズは $L \times 2M$ である。実装上、 $L$ は信号長の $1/2 \sim 1/3$ 程度に取れば十分である。例えば、標本化周波数 44100 Hz で 0.1 秒の信号を考えると、信号長の 1/3 に取るならば $L = 4410 / 3 = 1470$ であり、1/2 ならば $L = 4410 / 2 = 2205$ である。

※：固有値が $2M$ 個なのは、後述するように実信号という制約から正の周波数と負の周波数が $M$ 個ずつ同数現れるからである。

直交性の原理

MUSIC 法の拠り所となる、信号部分空間 $\mathbf{E}_s$ と雑音部分空間 $\mathbf{E}_n$ が近似的に直交すると見なせるという事実、いわば「直交性の原理」を以下では説明する。

雑音のない共分散行列を考える。

$\begin{align*} \mathbf{R}_{s} = \mathbf{A} \mathbf{S} \mathbf{A}^{*} \end{align*}$

ここで、 $\mathbf{A}$ はステアリングベクトルを並べた行列（サイズは $L \times M$ ）である。

$\begin{align*} \mathbf{A} &= [\mathbf{a}(\omega_1), \mathbf{a}(\omega_2), \ldots, \mathbf{a}(\omega_{M})] \\ \mathbf{a}(\omega_k) &= [1, \exp(j\omega_k), \exp(j2\omega_k), ..., \exp(j(L-1)\omega_k)]^{\top} \end{align*}$

また $\mathbf{S}$ は信号源（つまり各正弦波）の共分散行列（サイズは $M \times M$ ）である。

$\begin{align*} \mathbf{S} = \mathbb{E}[\mathbf{s}(n) \mathbf{s}(n)^{*} ] \end{align*}$

$\mathbf{s}(n) = [s_1(n), s_2(n), ..., s_M(n)]^{\top}$ は、サンプル時刻 $n$ における各正弦波を並べたベクトルである。この行列 $\mathbf{A} \mathbf{S} \mathbf{A}^{*}$ の値域（列空間） $\mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*})$ は、 $\mathbf{A}$ の列ベクトル、すなわちステアリングベクトル $\mathbf{a}(\omega_k)$ たちが張る空間であり、真の信号部分空間である（詳しくは下記の補足1を参照）。

定義から、 $\mathbf{R}_{s}$ の左零空間 $\mathcal{N}(\mathbf{R}_{s})$ は列空間 $\mathcal{C} (\mathbf{R}_{s})$ と互いに直交する。 $\mathbf{v}^{*} \mathbf{R}_{s} = \mathbf{0}$ となる $\mathbf{v} \in \mathcal{N}(\mathbf{R}_{s})$ は $\mathbf{v}^{*} \mathbf{a}(\omega_k) = 0$ を満たす（ $\mathbf{A}^{*}$ の階数が $M$ （補足1を参照）、かつ $\mathbf{S}$ が正則だから（補足2を参照））。つまり $\mathcal{N}(\mathbf{R}_{s})$ は、すべてのステアリングベクトル $\mathbf{a}(\omega_k)$ と直交する空間となる。

雑音のある共分散行列 $\mathbf{R}_{x} = \mathbf{R}_{s} + \sigma^{2} \mathbf{I}$ について、固有ベクトル $\mathbf{e}_{i}$ を考える。信号成分に対応する固有ベクトル（ $i \leq 2M$ ）は $\mathcal{C} (\mathbf{R}_{s})$ （信号部分空間 $\mathbf{E}_s$ ）を近似的に張り、雑音成分にのみ対応する固有ベクトル（ $i \gt 2M$ ）は $\mathcal{N}(\mathbf{R}_{s})$ （雑音部分空間 $\mathbf{E}_n$ ）を近似的に張る。

ゆえに、 $\mathbf{E}_n$ を張る固有ベクトルはステアリングベクトル $\mathbf{a}(\omega_k)$ と近似的に直交する。

$\begin{align*} \mathbf{a}(\omega_k)^{*} \mathbf{e}_{i} \approx 0 \end{align*}$

つまり $\mathbf{a}(\omega_k)^{*} \mathbf{E}_n \approx 0$ と見なすことができ、 $\mathbf{a}(\omega_k) \in \mathbf{E}_s$ であるから、 $\mathbf{E}_s$ と $\mathbf{E}_n$ が近似的に直交しているとも見なせる。このような直交性を利用する点が MUSIC 法のキモである。

補足1：共分散行列の列空間とステアリングベクトルたちが張る空間が一致することの証明

まず、一般の行列 $\mathbf{B}$ （サイズ $m \times n$ ）の列空間の定義を確認しておく。

$\begin{align*} \mathcal{C}(\mathbf{B}) = \{ \mathbf{y} \in \mathbb{C}^{m} \mid \mathbf{y} = \mathbf{B}\mathbf{x}, \; \mathbf{x} \in \mathbb{C}^{n} \} \end{align*}$

証明したいのは $\mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*}) = \mathcal{C}(\mathbf{A})$ である。

1． $\mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*}) \subset \mathcal{C}(\mathbf{A})$

任意の $\mathbf{y} \in \mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*})$ を取る。あるベクトル $\mathbf{x}$ が存在して

$\begin{align*} \mathbf{y} = (\mathbf{A} \mathbf{S} \mathbf{A}^{*}) \mathbf{x} \end{align*}$

と書ける。ここで、 $\mathbf{S} \mathbf{A}^{*} \mathbf{x}$ の部分を $\mathbf{z} = \mathbf{S} \mathbf{A}^{*} \mathbf{x}$ と置くと、

$\begin{align*} \mathbf{y} = \mathbf{A} \mathbf{z} \end{align*}$

とも書ける。すなわち $\mathbf{y} \in \mathcal{C}(\mathbf{A})$ であり、 $\mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*}) \subset \mathcal{C}(\mathbf{A})$ が示された。

2． $\mathcal{C}(\mathbf{A}) \subset \mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*})$

任意の $\mathbf{y} \in \mathcal{C}(\mathbf{A})$ を取る。あるベクトル $\mathbf{x}$ が存在して

$\begin{align*} \mathbf{y} = \mathbf{A} \mathbf{x} \end{align*}$

と書ける。信号源の共分散行列 $\mathbf{S}$ はサイズ $M \times M$ の対角行列であり、対角成分は各信号のパワー（つまり正値）である（補足2を参照）。ゆえに正則であり、逆行列が存在する。また $\mathbf{A}$ （サイズ $L \times M$ ）の列ベクトルは $M$ 本のステアリングベクトルたちであり、それぞれ異なるものとしているので線形独立である。ゆえに $\mathbf{A}$ の階数は $M$ である。 $\mathbf{A}$ の随伴 $\mathbf{A}^{*}$ （サイズ $M \times L$ ）も階数は $M$ である。

以上を頭に入れて、変形を進める。

$\begin{align*} \mathbf{y} &= \mathbf{A} \mathbf{x} \\ &= \mathbf{A} \mathbf{S} \mathbf{S}^{-1} \mathbf{x} \\ &= \mathbf{A} \mathbf{S} ( \mathbf{S}^{-1} \mathbf{x}) \end{align*}$

そこで $\mathbf{w} = \mathbf{S}^{-1} \mathbf{x}$ と置いて、 $\mathbf{y} = \mathbf{A} \mathbf{S} \mathbf{w}$ となる（ $\mathbf{w}$ は $M$ 次元ベクトル）。さらに、 $\mathbf{A}^{*}$ は階数 $M$ なので、 $\mathbf{A}^{*} \mathbf{v} = \mathbf{w}$ となるベクトル $\mathbf{v}$ が存在する（ $\mathbf{A}^{*}$ の像空間に関して全射）。この関係式を代入して、

$\begin{align*} \mathbf{y} & = \mathbf{A} \mathbf{S} \mathbf{A}^{*} \mathbf{v} \end{align*}$

を得る。ゆえに $\mathbf{y} \in \mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*})$ であり、 $\mathcal{C}(\mathbf{A}) \subset \mathcal{C}( \mathbf{A} \mathbf{S} \mathbf{A}^{*})$ が示された。

補足2：正弦波（信号源）の共分散行列が対角行列になることについて

いま正弦波 $s_k(n) = A_k \exp(j(\omega_k n + \phi_k))$ において、振幅と初期位相とが共に未知であり、 $A_k$ は $a_k^{(l)}, a_k^{(u)}$ の範囲で、 $\phi_k$ は $[-\pi, \pi]$ の範囲で一様に分布する確率変数であると仮定する。振幅の限界値 $a_k^{(l)} \: a_k^{(u)}$ はあらかじめ決めておくものとする。この仮定により、正弦波を一種の定常確率過程として扱うことができる。さらに追加の仮定として、

各正弦波 $k$ において、その振幅 $A_k$ と位相 $\phi_k$ は独立である
正弦波 $k$ のパラメータ $(A_k, \phi_k)$ は、正弦波 $l$ のパラメータ $(A_l, \phi_l)$ と互いに独立である（ただし $k \neq l$ ）

を設ける。これらの仮定の下で、まず $k \neq l$ の場合の相関の期待値を計算する。

$\begin{align*} \mathbb{E} [s_k(n) s_l(n)^{*}] &= \mathbb{E}[ A_k \exp(j(\omega_k n + \phi_k)) A_l \exp(-j(\omega_l n + \phi_l)) ] \\ &= \mathbb{E}[ A_k] \mathbb{E}[ A_l] \mathbb{E}[\exp(j(\omega_k n + \phi_k))] \mathbb{E}[\exp(-j(\omega_l n+ \phi_l))] \\ &= \mathbb{E}[ A_k] \mathbb{E}[ A_l] \exp(j\omega_k n) \mathbb{E}[\exp(j\phi_k)] \exp(-j\omega_l n) \mathbb{E}[\exp(-j\phi_l)] \end{align*}$

初期位相 $\phi_k$ , $\phi_l$ がそれぞれ $[-\pi, \pi]$ で一様分布していると仮定すれば、 $\mathbb{E}[\exp(j\phi_k)] =0$ , $\mathbb{E}[\exp(j\phi_l)]=0$ であるから、積全体もゼロである。

$\begin{align*} \mathbb{E} [s_k(n) s_l(n)^{*}] = 0 \end{align*}$

ある一つの正弦波における振幅と位相が確率変数であり、かつそれらが他の確率変数と独立ならば、異なる正弦波どうしは互いに無相関である、ということである。

$k = l$ の場合の相関の期待値はすなわち正弦波のパワー $P_k$ である：

$\begin{align*} P_k &= \mathbb{E} [ | s_k(n) |^{2} ] \\ &= \mathbb{E} [ |A_k \exp(j(\omega_k n + \phi_k))|^{2} ] \\ &= \mathbb{E} [ |A_k|^{2} ] \end{align*}$

以上から、共分散行列 $\mathbf{S} = \mathbb{E}[\mathbf{s}(n) \mathbf{s}(n)^{*} ]$ は対角行列であり、対角成分には各正弦波のパワーが並ぶことが分かった。

周波数推定の2つのアプローチ

雑音部分空間との直交性を利用して周波数を推定するにあたり、2つのアプローチが存在する。

spectral-MUSIC

直交性を評価するため、以下の MUSIC スペクトル $\widehat{P}_{MU}(\omega)$ の値を、あらかじめ用意した離散的な周波数グリッド上で計算する。

$\begin{align*} \widehat{P}_{MU}(\omega) = \frac{1}{\mathbf{a}(\omega)^{*} \mathbf{E}_n \mathbf{E}_n^{*} \mathbf{a}(\omega)} \end{align*}$

真の正規化角周波数 $\omega_k$ の位置 $\omega=\omega_k$ 付近では、 $\mathbf{a}(\omega)$ と $\mathbf{E}_n$ との直交性により分母がゼロに近づき、 $\widehat{P}_{MU}(\omega)$ に鋭いピークが立つ。 MUSIC スペクトルのピーク位置を探索することで、正規化角周波数 $ω_k$ の推定値 $\widehat{\omega}_k$ が得られる。

最後に、これを物理的な周波数 $\widehat{f}_k$ [Hz] の推定値に変換する。

$\begin{align*} \widehat{f}_k = \frac{F_s}{2\pi} \widehat{\omega}_k \end{align*}$

以上が spectral-MUSIC の概要である。

spectral-MUSIC の利点と欠点をまとめておく。

利点
- 直感的で理解しやすい、スペクトルを可視化できる。
- スペクトルの分母の計算に FFT（高速フーリエ変換）を活用することで、root-MUSIC よりも計算コストは小さく済む。
欠点
- 推定精度と計算コストが、周波数グリッドの密度に依存するトレードオフの関係にある。つまりグリッドを密に取れば高精度だが、その分だけ時間がかかるというわけである。
- 推定結果は、あらかじめ用意した周波数候補（グリッド）の中から最も近いものが選ばれるため、真の周波数がグリッドとグリッドの中間にあった場合、原理的に誤差（量子化誤差）が生じる。

root-MUSIC

root-MUSIC は、MUSIC スペクトルの分母 $Q(\omega)$ がゼロになる周波数を探索する問題を、ある多項式の根を求める問題に置き換えるアプローチである。

ステアリングベクトルを $z = \exp(j\omega)$ の関数 $\mathbf{a}(z) = [1, z^{-1}, ..., z^{-(L-1)}]^{\top}$ として表す。要するに $\omega$ から $z$ への変数変換である。 $L$ は部分空間の次元である。これを用いると、MUSIC スペクトルの分母 $Q(\omega)$ は以下のように書ける。

$\begin{align*} Q(\omega) &= D(z)\\ &= \mathbf{a}(z)^{*} \mathbf{C} \mathbf{a}(z) \end{align*}$

ここで $\mathbf{C}$ は雑音部分空間から計算される行列 $\mathbf{E}_{n} \mathbf{E}_{n}^{*}$ である。 $\mathbf{a}(z)^{*}$ は $[1, z, ..., z^{(L-1)}]$ となるため、 $D(z)$ を展開すると $z$ の負のべき乗と正のべき乗を含むローラン多項式となる。標準的な多項式ソルバは、負のべき乗を含まない通常の多項式を扱うため、 $D(z)$ に $z^{(L-1)}$ を掛けて、次数が $0$ から $2L-2$ までの多項式 $F(z)$ に変換する。

$\begin{align*} F(z) = z^{(L-1)} D(z) \end{align*}$

多項式 $F(z)$ の $k$ 次の項の係数 $d_k$ は、行列 $\mathbf{C}$ の要素 $C_{i,j}$ のうち、 $i-j=k -(L-1)$ となる全ての要素の和に等しい（下記の補足参照）。

ソルバによって求めた $F(z) = 0$ の根は $z_k = \exp(j\omega_k)$ の形をしているので、偏角を計算して正規化角周波数 $ω_k$ の推定値 $\widehat{\omega}_k$ を得る。あとは spectral-MUSIC と同様に、物理的な周波数 $\widehat{f}_k$ [Hz] に変換し、推定が完了する。

$\begin{align*} \widehat{f}_k = \frac{F_s}{2\pi} \widehat{\omega}_k \end{align*}$

root-MUSIC の利点と欠点をまとめておく。

利点
- スペクトル探索が不要なため、計算コストが周波数グリッドの密度に依存しない。
- 多項式の根を計算することで、周波数を連続値として直接算出するため、原理的に高い精度が期待できる。
欠点
- 信号に雑音が含まれる場合に、計算された多数の根の中から真の信号に対応する根を安定して選別するための、実装上の工夫（ヒューリスティクス）が必要となる。
- 標準的な多項式ソルバ（Numpyの polyroots など）の計算コストが高い。実際、計算量は多項式次数の3乗に比例するオーダーである。多項式次数 $L$ はスナップショットの窓長、すなわち部分空間の次元数であり、具体的には $L = 4410 / 3 = 1470$ 程度の値を取る。ゆえに、窓長 $L$ を短くすればするほど計算時間削減の効果は大きくなるが、同時にまた周波数推定に使えるデータの数が減少するので、推定精度も低下する。したがって、あまり $L$ を小さくしすぎない注意が必要である。

根の選別とヒューリスティクス

信号に雑音が含まれず、真の共分散行列が与えられるという理想的な条件では、所望の根は全て複素平面における単位円周上に乗っている。しかしながら、実際は有限サンプルから計算される共分散行列の推定精度や雑音の影響を受けるため、根の配置は円周上から外れる。

根の選別のための一つのヒューリスティクスを挙げると、円周から外れた互いに共役な根（正の周波数と負の周波数に対応）は単位円を挟んで配置されるので、まず単位円の内側から単位円周に最も近いものを順に $M$ 個以上選出し、その後にフィルタリングする方法がある。つまり雑音の影響で負の周波数になってしまったものを除く作業である。

その他のヒューリスティクス（今回の実装で採用）には、単純に単位円に最も近いものから順に複数個、つまり少なくとも $2M$ 個を選び出した上で、正の周波数のみを取り出す方法がある。雑音の入らないクリーンな信号に対しては、理想的には $2M$ 個ちょうどの根を選び出せば十分である（ $M$ 個の正周波数と $M$ 個の負周波数のペア）。しかしながら、現実には雑音混入の影響や推定精度の限界から、本来の周波数成分以外の成分、いわば偽の (spurious, false) 周波数成分に対応する余分な根が、単位円のかなり近くにしばしば発生する。それらの根が混入すると、抽出すべきであった所望の根が選出から漏れてしまう（根の見逃し）。つまり $2M$ 個の根を選出するだけでは不十分であり、所望の根が全て含まれるように過剰な個数の根、例えば $4M$ 個を選出しておく。このように根の候補を十分に確保したうえで、後処理によって根をフィルタリングする戦略である。

補足：root-MUSIC 多項式の係数導出

まず、 $D(z) = \mathbf{a}(z)^{*} \mathbf{C} \mathbf{a}(z)$ を展開する。

$\begin{align*} D(z) = \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} C_{i,j} (\mathbf{a}(z)^{*})_i (\mathbf{a}(z))_j \end{align*}$

ここで、 $C_{i,j}$ は行列 $\mathbf{C}$ の $(i, j)$ 成分である。 $(\mathbf{a}(z)^{*})_i$ と $\mathbf{a}(z)_j$ はそれぞれ、 $\mathbf{a}(z)^{*}$ の $i$ 番目の要素と $\mathbf{a}(z)$ の $j$ 番目の要素であり、それぞれ $z^{i}$ と $z^{-j}$ である。ゆえに $z$ のべき乗のみで改めて書くと、

$\begin{align*} D(z) = \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} C_{i,j} z^{(i-j)} \end{align*}$

となる。上式は $D(z)$ が $z$ の $-(L-1)$ 次から $(L-1)$ 次までの項を持つローラン多項式であることを示している。

この $D(z)$ に $z^{(L-1)}$ を掛けて、 $F(z)$ を作る。

$\begin{align*} F(z) &= z^{(L-1)} \left(\sum_{i=0}^{L-1} \sum_{j=0}^{L-1} C_{i,j} z^{(i-j)} \right) \\ &= \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} C_{i,j} z^{(i-j + L -1)} \\ &= \sum_{k=0}^{2L-2} d_{k} z^{k} \end{align*}$

$F(z)$ は次数が $0$ から $2L-2$ までの、負べきを含まない通常の多項式となっており、まさに知りたいのは $z^{k}$ の項の係数 $d_k$ を、 $C_{i, j}$ から計算する方法である。この係数 $d_k$ は、 $z$ の指数部分 $(i-j + L-1)$ が $k$ となるような全ての $C_{i,j}$ を足し合わせたものに等しくなるはずである（当たり前）。つまり

$\begin{align*} i-j + L-1 &= k \\ i-j &= k - (L-1) \end{align*}$

であり、この関係を満たす全ての $(i, j)$ のペア（ただし $0 \leq i, j \lt L$ ）についての $C_{i,j}$ の和が $d_k$ となる。 $L=3$ の場合を例に取ると、

$d_0 = C_{0,2}$
$d_1 = C_{0,1} + C_{1,2}$
$d_2 = C_{0,0} + C_{1,1} + C_{2,2}$
$d_3 = C_{1,0} + C_{2,1}$
$d_4 = C_{2,0}$

となっている。行列のインデックスを0スタートにしていることに注意する。

最小二乗法による振幅・位相の推定

MUSIC 法で周波数 $f_k$ が推定できれば、信号モデル $x(n) = \sum_{k=1}^{M} c_{k} \exp(j2\pi f_{k} n/F_s)$ は、未知数である複素振幅 $c_k = A_k \exp(j( 2\pi \phi_k))$ に関する線形方程式と見なすことができる。線形方程式を行列で表現すると

$\begin{align*} \mathbf{x} = \mathbf{V} \mathbf{c} \end{align*}$

となる。ここで $\mathbf{x}$ は観測信号ベクトル、 $\mathbf{V}$ は Vandermonde 行列、 $\mathbf{c}$ は求めるべき複素振幅ベクトルである。 $\mathbf{V}$ の各列は推定周波数 $f_k$ に対応するステアリングベクトル $\mathbf{a}(f_k)$ から成る。

この方程式は、 $\mathbf{V}$ の疑似逆行列 $\mathbf{V}^{+}$ を用いて、最小二乗解 $\hat{\mathbf{c}} = \mathbf{V}^{+} \mathbf{x}$ として解くことができる。得られた複素振幅 $\hat{\mathbf{c}}$ から、振幅 $A_k$ と位相 $\phi_k$ を抽出する。

$\begin{align*} A_k &= |\hat{c}_k|\\ \phi_k &= \arg \hat{c}_k \end{align*}$

分析対象が実信号の場合、振幅を2倍し、以下のように求める。絶対値を2倍する理由は下記の補足を参照してほしい。

$\begin{align*} A_k &= 2|\hat{c}_k|\\ \phi_k &= \arg \hat{c}_k \end{align*}$

補足：複素振幅ベクトルの絶対値を2倍する理由

以下，表記を簡単にするために連続時間信号 $t$ で考える。分析対象となる信号が、 $A \cos(2 \pi f t + \phi)$ という形をした複数のコサイン波の重ね合わせ（実信号）で書けるとする。オイラーの公式 $\exp(j\theta) = \cos(\theta) + j \sin(\theta)$ を使って、各コサイン波を複素正弦波の和に変換する：

$\begin{align*} A \cos(2\pi ft + \phi) = \frac{A}{2} \exp(j(2\pi ft + \phi)) + \frac{A}{2} \exp(-j(2\pi ft + \phi)) \end{align*}$

これは「周波数 $f$ のコサイン波は、実は正の周波数 $+f$ を持つ複素正弦波と、負の周波数 $-f$ を持つ複素正弦波のペアでできている」ということを意味する。しかし、MUSIC 法では通常、正の周波数 ( $f > 0$ ) のみを推定する。したがって、計算を簡略化するために信号モデルの正の周波数部分だけを使う。

$\begin{align*} y(t) = \frac{A}{2} \exp(j(2\pi ft + \phi)) \end{align*}$

少し変形して、振幅と位相を1つの複素振幅 $c$ にまとめる。

$\begin{align*} y(t) = \left\{ \frac{A}{2} \exp(j\phi) \right\} \exp(j2\pi ft) = c \exp(j2\pi ft) \end{align*}$

この $c$ が見つかれば、振幅 $A$ は 2 $|c|$ で求まり、位相 $\phi$ は $\arg c$ で求まるわけである。信号モデルと実信号の表現の違いに気をつける必要がある。

実装

以下に置いた。Enjoy!

github.com

推定器クラスの階層構造を図 1 に示す。

以下、これらのクラスについて説明する。

SinusoidParameters

周波数、振幅、位相のパラメータ群を保持する dataclass として SinusoidParameters を作成した。frozen=True とすることで、一度生成されたパラメータセットが不変であることを保証する。

@dataclass(frozen=True)
class SinusoidParameters:
    """A class to store the parameters of multiple sinusoids."""

    frequencies: npt.NDArray[np.float64]
    amplitudes: npt.NDArray[np.float64]
    phases: npt.NDArray[np.float64]

実装場所：utils/data_models.py
- このファイルにはサンプリング周波数や SNR など実験全体の設定を保持するための ExperimentConfig も定義してある。

基底クラス AnalyzerBase

パラメータ推定のうち共通の処理を、 AnalyzerBase クラス（基底クラス）によって実装する。SinusoidParametersクラスとの関係は、階層図において composition (has-a relationship) で示されている。

__init__: 解析に必要なパラメタ（サンプリング周波数、信号数など）で初期化する。
fit: 信号を受け取り、パラメタ推定の全プロセスの骨格に相当する。推定結果はインスタンス変数 self.est_params（型は SinusoidParameters）に格納される。
プロパティ (@property): fitの実行後、analyzer.frequencies のように直感的な形で推定結果にアクセスするための getter を提供する。これらは SinusoidParameters のメンバと一対一に対応する。
get_params: モデル固有のハイパパラメタ（MUSIC 法ならば周波数グリッド点の数など）の設定値を参照するためのメソッドである。階層図には書かれていないが、ハイパパラメタは AnalyzerParameters という TypedDict にも格納される。get_params は AnalyzerParameters に対するインタフェースを提供している。
プライベートメソッド (_で始まる): 派生クラスで必須となるメソッドを実装する。
- _build_covariance_matrix: 共分散行列の計算を行う。
- _estimate_amplitudes_phases : 周波数推定後の振幅・位相の推定を行う。
- _estimate_frequencies: 周波数推定のプライベートメソッドであるが、この時点では抽象メソッドとしておき、具体的な実装（オーバーライド）は派生クラスで行う。
実装場所：analyzers/base.py

基底クラス：MusicAnalyzerBase

次に、MUSIC 法に特有かつ共通の処理を、 AnalyzerBase クラスを継承した MusicAnalyzerBase クラス（基底クラス）によって実装する。基底クラスを2段階で用意した理由は、ESPRIT 法の実装時に AnalyzerBase クラスを継承した基底クラス（EspritAnalyzerBase）を用意するからであり、このほうが設計的に美しい。

_estimate_noise_subspace : 雑音部分空間の推定を実装する。
実装場所：analyzers/music/base.py

Mixin クラス：ForwardBackwardMixin

共分散行列の推定精度を高めるための Forward-Backward Averaging [5] を実装した。基底クラス MusicAnalyzerBase は共分散行列を推定するために _build_covariance_matrix というプライベートメソッドを用いるが、派生クラスが ForwardBackwardMixin を Mixin の形で多重継承することによって、このメソッドをオーバーライドする。

実装場所：mixins/covariance.py

派生クラス：SpectralMusicAnalyzer

SpectralMusicAnalyzer は MusicAnalyzerBase クラスを継承した派生クラスである。_estimate_frequencies メソッドを、スペクトル探索のアプローチで具体的に実装（オーバーライド）する。このクラス固有のプライベートメソッドは

_calculate_music_spectrum: MUSICスペクトルを計算するメソッド

である。

実装場所：analyzers/music/spectral.py

派生クラス：SpectralMusicAnalyzerFB

SpectralMusicAnalyzerFB は ForwardBackwardMixin と SpectralMusicAnalyzer を Mixin の形で多重継承したクラスである。

class SpectralMusicAnalyzerFB(ForwardBackwardMixin, SpectralMusicAnalyzer):
    """Spectral MUSIC analyzer enhanced with Forward-Backward averaging.

    Inherits from ForwardBackwardMixin to override the covariance matrix calculation.
    """

派生クラス：RootMusicAnalyzer

RootMusicAnalyzer もまた MusicAnalyzerBase クラスを継承した派生クラスである。 _estimate_frequencies メソッドを、多項式求根のアプローチで具体的に実装（オーバーライド）する。このクラス固有のプライベートメソッドは

_calculate_polynomial_coefficients: 多項式の係数を計算するメソッド

である。求根処理と後処理は、内部的に用意した別の関数を呼び出す形で実装する（numpy の polyroots をラップ）。

実装場所：analyzers/music/root.py

派生クラス：RootMusicAnalyzerFB

RootMusicAnalyzerFB は ForwardBackwardMixin と RootMusicAnalyzer を Mixin の形で多重継承したクラスである。

class RootMusicAnalyzerFB(ForwardBackwardMixin, RootMusicAnalyzer):
    """Root MUSIC analyzer enhanced with Forward-Backward averaging.

    Inherits from ForwardBackwardMixin to override the covariance matrix calculation.
    """

実験的評価

ここでは MUSIC 法の性能を簡単な実験を通じて評価する。

デモスクリプトをコマンドラインから動かして実験を行う上で、コマンドラインオプションを用意してある。例えば標本化周波数 44100 Hz、信号長 0.1 s、SNR 40 dB、正弦波の周波数が440 Hz, 445 Hz, 450 Hz、振幅の範囲を 0.5 〜1.5、周波数探索時におけるグリッド点数を 16384 にしたい場合は

python3 examples/run_comparison.py --fs 44100 --duration 0.1 --snr_db 40 \
              --freqs_true "440 445 450" --amp_range "0.5 1.5" --n_grids 16384

のようにする。リポジトリに置いてある main.py には、既に実装済である ESPRIT 法（次回の記事で紹介）も実行できるようにしてある。ゆえに以下に示す結果は実行結果の抜粋となる。

実験その1

MUSIC 法の性能を評価するため、以下の設定でシミュレーションを実施する。SNR は 40 dB とやや高めに設定した。

信号: 3つの近接した周波数 (440 Hz, 445 Hz, 450 Hz) を持つ正弦波の和
サンプリング周波数: 44100 Hz
信号長: 0.1 sec (FFT分解能: 10 Hz)
各正弦波の振幅値: 0.5〜1.5の範囲でランダム
SNR: 40 dB
グリッド点数（周波数探索時の分解能）: 16384

実験結果を以下に示す。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                40.0 dB
True Frequencies:   [440. 445. 450.] Hz
True Amplitudes:    [1.30039923 0.98442246 0.87698843]
True Phases:        [-0.94952831  1.01071672  1.09965033] rad

--- Running Spectral MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  N Grids: 16384

--- Estimation Results ---
Est Frequencies: [440.11170115 444.14942318 449.53305255] Hz
Est Amplitudes:  [1.47830862 1.07005009 0.96662503]
Est Phases:      [-1.00671587  1.09819513  1.14682031] rad

--- Estimation Errors ---
Freq Errors:  [ 0.11170115 -0.85057682 -0.46694745] Hz
Amp Errors:   [0.17790939 0.08562763 0.0896366 ]
Phase Errors: [-0.05718755  0.08747841  0.04716998] rad


--- Running Root MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333

--- Estimation Results ---
Est Frequencies: [440.00599141 444.83948865 450.03890814] Hz
Est Amplitudes:  [1.31482601 0.97266724 0.87505361]
Est Phases:      [-0.95566543  1.05456075  1.05716061] rad

--- Estimation Errors ---
Freq Errors:  [ 0.00599141 -0.16051135  0.03890814] Hz
Amp Errors:   [ 0.01442678 -0.01175522 -0.00193482]
Phase Errors: [-0.00613712  0.04384403 -0.04248972] rad

実験結果から、MUSIC 法が FFT の分解能限界 (10 Hz) よりもはるかに狭い 5 Hz間隔の周波数を分離・推定できていることが確認できる。これは、MUSIC 法が持つ超解像性能を示すものである。また、root-MUSIC 法による周波数の推定精度が spectral-MUSIC 法の推定精度を大きく上回ったことも確認できる。振幅と位相も問題なく推定できている。

実験その2

SNRを20 dB まで小さくし，他の条件は実験その1と揃えた状態で実験を行う。

信号: 3つの近接した周波数 (440 Hz, 445 Hz, 450 Hz) を持つ正弦波の和
サンプリング周波数: 44100 Hz
信号長: 0.1 sec (FFT分解能: 10 Hz)
SNR: 20 dB
各正弦波の振幅値: 0.5〜1.5の範囲でランダム
グリッド点数: 16384

実験結果を示す。ノイズの影響を受け、推定に失敗している。周波数が近接し、かつ MUSIC スペクトルがノイズに埋もれてしまった結果、ありえない周波数を誤検出したと考えられる。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                20.0 dB
True Frequencies:   [440. 445. 450.] Hz
True Amplitudes:    [1.16982233 0.8748777  0.56716902]
True Phases:        [-0.13502604  0.65792432 -1.94954471] rad


--- Running Spectral MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  N Grids: 16384

--- Estimation Results ---
Est Frequencies: [  440.11170115   448.18714521 18164.36550082] Hz
Est Amplitudes:  [0.857592   1.05505536 0.00715926]
Est Phases:      [ 0.03896099 -0.81346753 -1.68694601] rad

--- Estimation Errors ---
Freq Errors:  [1.11701154e-01 3.18714521e+00 1.77143655e+04] Hz
Amp Errors:   [-0.31223034  0.18017766 -0.56000976]
Phase Errors: [ 0.17398703 -1.47139185  0.2625987 ] rad


--- Running Root MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333

--- Estimation Results ---
Est Frequencies: [  439.54740547   448.47439399 18164.77455859] Hz
Est Amplitudes:  [0.85499289 1.04188646 0.00730656]
Est Phases:      [ 0.25313974 -0.96671257 -1.81650774] rad

--- Estimation Errors ---
Freq Errors:  [-4.52594526e-01  3.47439399e+00  1.77147746e+04] Hz
Amp Errors:   [-0.31482944  0.16700876 -0.55986246]
Phase Errors: [ 0.38816578 -1.62463689  0.13303697] rad

実験その3

ここでは SNR を実験その2と同様の 20 dB にして、かつ正弦波の周波数を離した状態（10 Hz 間隔）で推定させる。

信号: 3つの近接した周波数 (440 Hz, 450 Hz, 460 Hz) を持つ正弦波の和
サンプリング周波数: 44100 Hz
信号長: 0.1 sec (FFT分解能: 10 Hz)
各正弦波の振幅値: 0.5〜1.5の範囲でランダム
SNR: 20 dB
グリッド点数（周波数探索時の分解能）: 16384

実験結果を示す。多少の誤差はあるものの、3つの周波数の推定に成功している。SNR が小さくとも、周波数を十分に離せば推定できることが分かった。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                20.0 dB
True Frequencies:   [440. 450. 460.] Hz
True Amplitudes:    [0.53970806 1.48169218 1.27532951]
True Phases:        [ 0.97787732 -2.11985235 -1.02808099] rad


--- Running Spectral MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  N Grids: 16384

--- Estimation Results ---
Est Frequencies: [440.11170115 449.53305255 460.3003113 ] Hz
Est Amplitudes:  [0.47391103 1.46436949 1.24914052]
Est Phases:      [ 0.97046465 -1.99518855 -1.07836629] rad

--- Estimation Errors ---
Freq Errors:  [ 0.11170115 -0.46694745  0.3003113 ] Hz
Amp Errors:   [-0.06579703 -0.01732269 -0.02618899]
Phase Errors: [-0.00741267  0.1246638  -0.0502853 ] rad


--- Running Root MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333

--- Estimation Results ---
Est Frequencies: [439.60582965 449.47013942 460.17414509] Hz
Est Amplitudes:  [0.46274031 1.49516184 1.24895291]
Est Phases:      [ 1.10880204 -1.96690338 -1.03914951] rad

--- Estimation Errors ---
Freq Errors:  [-0.39417035 -0.52986058  0.17414509] Hz
Amp Errors:   [-0.07696775  0.01346966 -0.0263766 ]
Phase Errors: [ 0.13092472  0.15294897 -0.01106852] rad

実験その4

ここでは周波数探索のグリッド点を、実験その1やその3から減らしてみる。spectral-MUSIC 法のみを対象とする。

信号: 3つの近接した周波数 (440 Hz, 450 Hz, 460 Hz) を持つ正弦波の和
サンプリング周波数: 44100 Hz
信号長: 0.1 sec (FFT分解能: 10 Hz)
各正弦波の振幅値: 0.5〜1.5の範囲でランダム
SNR: 20 dB
グリッド点数（周波数探索時の分解能）: 4096

実験結果を示す。3つめの正弦波（460 Hz）の周波数の推定には失敗していることが分かる。それゆえ、グリッド点数は十分に確保しておく必要がある。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                20.0 dB
True Frequencies:   [440. 450. 460.] Hz
True Amplitudes:    [0.67367712 0.68575736 1.4302683 ]
True Phases:        [-2.77798553 -2.23526847 -1.76869682] rad

--- Running Spectral MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  N Grids: 4096
/home/***/work/esprit/music_param_class_demo.py:95: UserWarning: Expected 3 components, but found 2.
  warnings.warn(
Estimation incomplete or failed. Found 2 components.
Est Frequencies: [441.53846154 452.30769231] Hz

実験その5

spectral-MUSIC 法と root-MUSIC 法の実行時間を比較する。main 関数にて以下のコードを使って計測した。使用 OS とマシンは Ubuntu 24.04 (Intel Core i9-9900K) と MacBook Air (M1) である。サンプリング周波数は 44100 Hz 、信号長は 100 ms 、正弦波の個数は 3 とした。また Python のバージョンは Ubuntu が 3.12 で MacBook Air が 3.13 である。Numpy のバージョンは 2.3.2、Scipy のバージョンは 1.16.1 である。

import time 
...
start = time.perf_counter()
analyzer.fit(noisy_signal.astype(np.complex128))
end = time.perf_counter()
print(end - start)

それぞれの実行時間を以下に示す。Ubuntu と MacBook Air のどちらにおいても、spectral-MUSIC 法の実行時間が圧倒的に短い結果となった。root-MUSIC 法のボトルネックは多項式の求解にある。今回の実装では、求解に numpy.polynomial.polynomial.polyroots を採用したが、内部的にコンパニオン行列の固有値分解（計算量は多項式の次数の3乗に比例）を利用しているため、実行時間で大きく差が開いた。MUSIC 法と同様の枠組みに沿いつつ計算量を削減したい場合には、Min-Norm 法 [6] がある。先のリポジトリには既に実装済である。Ubuntu 24.04 上で動かしたところ，root Min-Norm の実行時間は root-MUSIC 法の概ね 30 % （約 2.5 秒）となり、大幅な計算量削減を達成している。

手法	Ubuntu 24.04 (Intel Core i9-9900K)	MacBook Air (M1)
spectral-MUSIC 法（グリッド点数 4096）	0.4924 秒	1.0125 秒
spectral-MUSIC 法（グリッド点数 8192）	0.5888 秒	1.1289 秒
spectral-MUSIC 法（グリッド点数 16384）	0.7389 秒	1.4161 秒
spectral-MUSIC 法（グリッド点数 32768）	1.1509 秒	1.9947 秒
root-MUSIC 法	8.8604 秒	15.09276 秒

実験その6

共分散行列の推定精度を改善するための Forward-Backward Averaging の導入前後で推定精度がどのように変化するかを検証する。実験条件は実験その 1 よりも SNR を小さくし、信号長を短くしている。

信号: 3つの近接した周波数 (440 Hz, 445 Hz, 450 Hz) を持つ正弦波の和
サンプリング周波数: 44100 Hz
信号長: 0.04 (FFT分解能: 10 Hz)
各正弦波の振幅値: 0.5 〜1.5 の範囲でランダム
SNR: 30 dB
グリッド点数（周波数探索時の分解能）: 16384

実験結果を以下に示す。Forward-Backward Averaging によって推定精度が改善したことが分かる。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    40 ms
SNR:                20.0 dB
True Frequencies:   [440. 460. 480.] Hz
True Amplitudes:    [0.75806191 1.06189213 0.56268528]
True Phases:        [-3.01681205 -1.1851941  -1.16639541] rad

--- Running Spectral MUSIC ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  N Grids: 16384

--- Estimation Results ---
Est Frequencies: [441.4576085  458.95440396 479.1430141 ] Hz
Est Amplitudes:  [0.83273335 1.14350704 0.5834617 ]
Est Phases:      [-3.12138567 -1.13723312 -1.11144532] rad

--- Estimation Errors ---
Freq Errors:  [ 1.4576085  -1.04559604 -0.8569859 ] Hz
Amp Errors:   [0.07467143 0.08161491 0.02077642]
Phase Errors: [-0.10457362  0.04796098  0.05495009] rad


--- Running Spectral MUSIC (Forward-Backward) ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  N Grids: 16384

--- Estimation Results ---
Est Frequencies: [440.11170115 460.3003113  479.1430141 ] Hz
Est Amplitudes:  [0.77239811 1.10557283 0.60856771]
Est Phases:      [-3.02878033 -1.24644535 -1.02350208] rad

--- Estimation Errors ---
Freq Errors:  [ 0.11170115  0.3003113  -0.8569859 ] Hz
Amp Errors:   [0.0143362  0.04368069 0.04588243]
Phase Errors: [-0.01196827 -0.06125125  0.14289333] rad

考察

実験結果は、近接した周波数成分を含む信号のパラメータ（周波数、振幅、位相）を MUSIC 法が高精度に推定できることを示している。特に、MUSIC 法が FFT の分解能限界を超える近接周波数（5 Hz 間隔）を正確に分離した上で、後続の最小二乗法が正確な振幅・位相を算出している点が重要である（実験その 1）。

実験その 3 の結果から、10 Hz 間隔の周波数の検出は SCDE でもできなかったので、大きなメリットである。しかし、周波数が近接した場合はノイズの影響を受けやすく（実験その 2）、また十分なグリッド点数を確保しなければならない（実験その 4）。グリッド点数を多くすれば、それだけ多くの計算時間を要するので、推定性能と計算時間の間にはトレードオフが存在する（実験その 4 および実験その 5）。

おわりに

本記事では、MUSIC 法の原理を解説し、Python による実装例と簡単な結果を示した。MUSIC 法は、信号の数が既知であるという強い仮定を置く一方で、高い周波数分解能を実現できる強力な手法である。次回の記事では、MUSIC 法と同じく部分空間法に分類されるが、スペクトル探索を必要としない ESPRIT 法について紹介する予定である。

MUSIC 法の計算コストは root-MUSIC はもとより、spectral-MUSIC も決して低いものではない。MUSIC 法の派生で、その計算コストを大幅に削減する Min-Norm 法 [6] はリポジトリに実装済であり、いずれ紹介記事を書く予定である。

参考文献

[1] Schmidt, R. O. (1979). “Multiple emitter location and signal parameter estimation,” in Proc. RADC, Spectral Estimation Workshop, Rome, NY, pp. 243–258.

[2] Bienvenu, G. (1979). “Influence of the spatial coherence of the background noise on high resolution passive methods,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing, Washington, DC, pp. 306–309.

[3] R.O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas and Propagat., vol. AP-34, no. 3, pp. 276-280, 1986.

[4] A. Barabell, "Improving the resolution performance of eigenstructure-based direction-finding algorithms," ICASSP '83. IEEE International Conference on Acoustics, Speech, and Signal Processing, Boston, MA, USA, 1983, pp. 336-339, doi: 10.1109/ICASSP.1983.1172124.

[5] Petre Stoica and Randolph L. Moses, "Spectral Analysis of Signals," Pearson Prentice Hall, 2005.

[6] R. Kumaresan and D. W. Tufts, "Estimating the Angles of Arrival of Multiple Plane Waves," in IEEE Transactions on Aerospace and Electronic Systems, vol. AES-19, no. 1, pp. 134-139, 1983.

[7] 唐沢好男，"高分解能到来方向推定の仕組み　これでわかる MUSIC 法と ESPRIT 法" http://www.radio3.ee.uec.ac.jp/ronbun/TR_YK_049_DOA_Estimation.pdf

[8] 浅野太，"音のアレイ信号処理," コロナ社，2011. https://www.coronasha.co.jp/np/isbn/9784339011166/