前回記事で紹介した論文の続編。

アトラクタの安定性が最大化される領域を発見。その起源を考察し、重み行列の「スペクトル集中」にあることを導出。

その領域において、「力」のせめぎ合いによる、自己組織化現象を確認。その場所は不安定ぽく見えるのだけど、実は最も記憶に適した場所だったという。

※ 2026年4月10日追記：NOLTA誌にアクセプト

2025-12-02

高容量カーネルHopfieldネットワークの性能を定量評価したプレプリントを公開

ホップフィールドネットワーク

KLRによる学習を導入することで、記憶容量 $P/N > 4.0$ になることを確認した（ $P$ はパターン数， $N$ はニューロン数）。

さらに、ネットワークサイズ $N$ と最適なカーネル幅 $\gamma$ の間に明確なスケーリング則が存在することを発見した。

arxiv.org

※ 2026年4月10日追記：NOLTA誌にアクセプト

2025-11-07

日本語5母音の「声道断面積関数」をPythonで可視化してみた

音声 python

はじめに

とあることがきっかけで，日本語5母音の声道断面積関数をプロットしてみたくなったので，やってみたということ．

荒井先生の論文には5母音の声道断面の直径パラメタが載っているので（Table 1），これを参考にする．

www.jstage.jst.go.jp

直径パラメタから断面積を計算し，scipyの PchipInterpolator によって補完する．

実装

以下に置いた．Enjoy!

https://gist.github.com/tam17aki/1adf19685e0dbd367fab606ab8d4c6ab

結果

図1から図5に結果を示す．

おわりに

可視化は楽しい．

2025-10-11

Mac で GTAGS for Python & Emacs

emacs python

自分用メモ

brew install global

cp /opt/homebrew/share/gtags/gtags.conf ~/.globalrc

.globalrcに対して、

common:\ :skip= ...

の最後に

mypy_cache/,ruff_cache/,*.pyc,venv/:

など、タグづけにあたりスキップしたいディレクトリや拡張子を追記

pythonファイルにタグ付けするためにPygmentsをインストールしておく

pip install Pygments

タグ付け

gtags --gtagslabel=pygments

gtags-modeのために以下を評価

https://github.com/Ergus/gtags-mode

(with-eval-after-load 'gtags-mode
  (setq gtags-mode-update-args "--gtagslabel=pygments")
  (defun gtags-mode--imenu-advice ()
    "Make imenu use Global."
    (when (and buffer-file-name (gtags-mode--local-plist default-directory))
      (gtags-mode--filter-find-symbol
       '("--file") (buffer-file-name)
       (lambda (name _code _file line)
         (list name line #'gtags-mode--imenu-goto-function))))))

2025-10-02

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

信号処理 ESPRIT python

はじめに

前回記事から引き続き、本記事も「正弦波のパラメタ推定シリーズ」である。

今回は、高分解能な周波数推定手法のひとつであるESPRIT (Estimation of Signal Parameters via Rotational Invariance Techniques) アルゴリズムについて解説し、Python によるクラスベースの実装を紹介する。ESPRIT 法は、前回記事で紹介した MUSIC 法と同じく部分空間法に分類されるが、スペクトル探索を必要とせず、信号部分空間の回転不変性を利用して直接パラメタを計算する点で大きく異なる。前回記事と同様、周波数の推定後に、正弦波の振幅と位相も最小二乗法を用いて推定する。それゆえ記事のタイトルに再び「複数正弦波のパラメタ推定」が現れている。

はじめに
ESPRIT 法の理論的背景
最小二乗法による振幅・位相の推定
実装
実験
考察
おわりに
参考文献

ESPRIT 法の理論的背景

問題設定（信号モデル）と部分空間

ここは前回記事の MUSIC 法と同様の定式化であり、記号の確認程度の目的でさらっと済ませる。

解析対象の信号が $M$ 個の複素正弦波と白色雑音の和で構成されると仮定する。計算機で扱うために、この信号を標本化周波数 $F_s$ [Hz] でサンプリングした離散信号 $x[n]$ を考える。

$\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*}$

ここで、 $A_k$ は振幅、 $\phi_k$ は初期位相、 $f_k$ は物理的な周波数 [Hz] であり、 $w[n]$ は離散時間の白色雑音である。正規化角周波数 $\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*}$

ここで、 $c_k = A_k \exp(j\phi_k)$ は複素振幅である。 $\omega_k$ は、「サンプルが1つ進むごとに位相がどれだけ変化するか」を表す量であり、単位は [rad/sample] である。ESPRIT 法が直接的に推定するのは、この正規化角周波数 $\omega_k$ である。

上記の目的を達成するため、信号の共分散行列 $\mathbf{R}_x$ を解析する。 $w[n]$ が白色雑音であるという仮定の下では、 $\mathbf{R}_x$ の固有値は2つのグループに分かれるのであった：

信号部分空間 $\mathbf{E}_s$ ： $M$ 個の大きな固有値に対応する固有ベクトルが張る空間。信号と雑音の両方のエネルギーを含む。
雑音部分空間 $\mathbf{E}_n$ ：残りの小さな固有値に対応する固有ベクトルが張る空間。雑音のエネルギーのみを含む。

分析対象となる正弦波が実信号の場合、 $\mathbf{E}_s$ は $2M$ 個の固有値（大きさ順）に対応した $2M$ 本の固有ベクトルによって張られる。 $2M$ 個の固有値を考える理由は、 $M$ 個の正の周波数成分に加えて、 $M$ 個の負の周波数成分（複素共役）を含ませるためである。

ESPRIT 法は、この信号部分空間 $\mathbf{E}_s$ を、雑音のない理想的な信号が張る部分空間の「推定値」として利用する。

回転不変性について

信号部分空間が持つ「回転不変性」を理解するため、雑音のない単一の複素正弦波 $s[n]$ を考える。

$\begin{align*} s[n] = A \exp(j(\omega n + \phi)) \end{align*}$

この信号を1サンプルだけ時間シフトさせると、 $s[n+1]$ は以下のようになる。

$\begin{align*} s[n+1] &= A \exp(j(\omega (n+1) + \phi))\\ &= A \exp(j(\omega n + \phi + \omega))\\ &= { A \exp(j(\omega n + \phi)) } \exp(j\omega) \\ & = s[n] \exp(j\omega) \end{align*}$

この関係式は、「信号を1サンプル進める操作は、複素平面上で $\exp(j\omega)$ という係数を掛けてベクトルを回転させる操作と等価である」ことを意味している。

部分空間への拡張

上記の関係式は、信号部分空間全体に対しても拡張できる。信号部分空間を張る基底ベクトルは、信号 $s[n]$ の線形結合で表現される。簡単のため、信号部分空間を $L$ 次元のベクトルで表現し、それを $\mathbf{v}$ とする。

$\begin{align*} \mathbf{v} = [s[0], s[1], ..., s[L-1]]^{\top} \end{align*}$

いま、

最終要素を除いたベクトル $\mathbf{v}_{1} = [s[0], s[1], ..., s[L-2]]^{\top}$
先頭要素を除いたベクトル $\mathbf{v}_{2} = [s[1], s[2], ..., s[L-1]]^{\top}$

を考えると、各要素の関係 $s[n+1] = s[n] \exp(j\omega)$ から、ベクトル全体としても以下の関係が成り立つ。

$\begin{align*} \mathbf{v}_2 =\mathbf{v}_1 \exp(j\omega) \end{align*}$

信号が $M$ 個の正弦波の和である場合、この関係は行列形式に拡張され、さらには信号部分空間全体に対しても関係は拡張できる。推定された信号部分空間 $\mathbf{E}_s$ から、

最終行を除いた部分行列 $\mathbf{E}_1$
先頭行を除いた部分行列 $\mathbf{E}_2$

を生成すると、 $M$ 個の回転因子 $\exp(j\omega_k)$ を対角成分に持つ対角行列 $\boldsymbol{\Phi}$ を用いて、以下の関係式が成り立つ。

$\begin{align*} \mathbf{E}_2 \approx \mathbf{E}_1 \boldsymbol{\Phi} \end{align*}$

近似的に等号が成り立つとしているのは、実際には雑音が加わったうえでの「推定」だからである。雑音が一切含まれない場合、信号部分空間 $\mathbf{E}_s$ は、雑音のない真の信号が張る部分空間と完全に一致する。それゆえ関係式は等号で成り立つ。なお等号で成り立つことの厳密な証明は、観測ベクトルを Vandermonde 行列 $\mathbf{V}$ を使って表現し、 $\mathbf{V}$ に対する回転不変性に基づいて示される。しかし長くなるので本記事では割愛し、記事を改めて紹介したい。

補足：最小二乗法と総最小二乗法（Total Least Squares; TLS）について

最小二乗法

$\mathbf{E}_2 \approx \mathbf{E}_1 \boldsymbol{\Phi}$ という近似的な関係式があった。最小二乗法は、この関係を以下のようにモデル化する。

$\begin{align*} \mathbf{E}_2 = \mathbf{E}_1 \boldsymbol{\Phi} + \boldsymbol{\Delta} \end{align*}$

$\mathbf{E}_1$ （説明変数）は正確で、誤差は全く含まれておらず、すべての誤差は $\mathbf{E}_2$ （従属変数）側に残差 $\boldsymbol{\Delta}$ として押し付けられるとしたものである。この残差 $\boldsymbol{\Delta}$ の大きさ（フロベニウスノルム）を最小化するような $\boldsymbol{\Psi}$ を見つけるのが目的である。

総最小二乗法（TLS）

TLSでは、雑音は $\mathbf{E}_s$ 全体に影響を与えているのだから、 $\mathbf{E}_1$ も $\mathbf{E}_2$ も両方とも誤差を含んでいるはずだ、と仮定する。

$\begin{align*} \mathbf{E}_2 + \boldsymbol{\Delta}_2 = (\mathbf{E}_1 + \boldsymbol{\Delta}_1) \boldsymbol{\Phi} \end{align*}$

ゆえに最適化の目的が、 $[\boldsymbol{\Delta}_1 \mid \boldsymbol{\Delta}_2]$ という結合された誤差行列の大きさを最小化するような $\boldsymbol{\Psi}$ を見つけることへと変わる。その解は $[\mathbf{E}_1 \mid \mathbf{E}_2]$ という拡張行列の特異値分解を通じて得られる。

一般に、SNR が十分に高い場合は最小二乗法と TLS の結果に大きな差はないとされている。しかし、SNR が低い環境では，雑音の影響をより正確にモデル化している TLS の方が、頑健で精度の高い推定結果を与えるとされている。

以下、周波数推定の説明は最小二乗法を採用したものとして進める。

回転不変性を利用した周波数推定

上記の関係式 $\mathbf{E}_2 \approx \mathbf{E}_1 \boldsymbol{\Phi}$ を、 $\boldsymbol{\Phi}$ の推定値である回転作用素 $\boldsymbol{\Psi}$ に関する方程式とみなし、 $\boldsymbol{\Psi}$ について近似的に解く（ $\mathbf{E}_2 \approx \mathbf{E}_1 \boldsymbol{\Psi}$ ； $\hat{\boldsymbol{\Phi}} = \boldsymbol{\Psi}$ ）。その解は疑似逆行列 $\mathbf{E}_{1}^{+}$ を用いて以下のように書ける。

$\begin{align*} \boldsymbol{\Psi} = \mathbf{E}_1^{+} \mathbf{E}_2 \end{align*}$

$\boldsymbol{\Phi}$ は $\exp(j\omega_k)$ を対角成分に持つ対角行列であった。したがって、推定された回転作用素 $\boldsymbol{\Psi}$ の固有値 $\lambda_k$ が、 $\exp(j\omega_k)$ の推定値となる。

$\begin{align*} \lambda_k \approx \exp(j\omega_k) \end{align*}$

よって $\boldsymbol{\Psi}$ の固有値を計算し、その偏角を取ることで、正規化角周波数 $\omega_k$ の推定値が直接的に求まる。

$\begin{align*} \widehat{\omega}_k = \arg \lambda_k \end{align*}$

最後に、推定された正規化角周波数 $\widehat{\omega}_k$ [rad/sample] を、標本化周波数 $F_s$ [sample/sec] を用いて物理的な周波数 $\widehat{f}_k$ [Hz] に変換する。

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

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

MUSIC 法のときと同様に、推定済の周波数から、振幅と位相を最小二乗法に基づいて推定する。この方法論はすでに前回記事で紹介したので、本記事では割愛する。

実装

以下に置いた。前回と同様のリポジトリである。

github.com

今回の記事に関連した推定器クラスの階層構造を図 1 に示す。

なおパラメータを格納するための dataclass である SinusoidParameters は MUSIC 法の実装と同様なので説明を割愛する。また全ての推定器の基底クラスである AnalyzerBase クラスについても説明を割愛する。

基底クラス：EspritAnalyzerBase

まず MUSIC 法と共通の基底クラス AnalyzerBase を継承した EspritAnalyzerBase によって、ESPRIT 法に共通の処理を実装する。これは ESPRIT 法の派生手法である Unitary ESPRIT 法や FFT-ESPRIT法（次回以降の記事で紹介）を既に実装しており、その実装における処理の共通化を図ったということである。具体的には推定された正規化角周波数 [rad/sample] を変換・ソートして物理的な周波数 [Hz] へと変換する処理が実装されている。

実装場所：base.py

基底クラス：EVDBasedEspritAnalyzer

ESPRIT 法の派生手法のうち，固有値分解および特異値分解に基づく手法を実装するために作成された基底クラスである。信号部分空間 $\mathbf{E}_s$ を推定するための _estimate_signal_subspace というプライベートメソッドが用意されている。

派生クラス：StandardEspritAnalyzer

EVDBasedEspritAnalyzerを継承し、 StandardEspritAnalyzer クラスとしてカプセル化する。_estimate_signal_subspace は具体的にこのクラスで実装される。

実装場所：standard.py

周波数を推定した後、MUSIC 法の前回記事で紹介したものと全く同じように、最小二乗法に基づいて振幅と位相を算出する。

派生クラス：StandardEspritAnalyzerFB

MUSIC 法と同様に、ESPRIT 法でも共分散行列の推定に前方後方平均化が利用できる。そのため、ForwardBackwardMixin クラスをMixinの形で多重継承することで StandardEspritAnalyzerFB クラスを定義している。

ソルバクラス：LSEspritSolver, TLSEspritSolver

LS と TLS による周波数推定機能は別クラスに分離した。

実装場所：solvers.py

Unitary ESPRIT 法の実装でも LS および TLS に基づく推定を行っており、共通の処理をソルバとして分離したわけである。ソルバは StandardEspritAnalyzer のインスタンス変数として保持する。

実験

実験その1

まずは MUSIC 法の前回記事・実験その1と条件を揃えて、以下の設定でシミュレーションを実施する。

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

実験結果を以下に示す。5 Hz 間隔の周波数は問題なく推定できている。推定精度も MUSIC 法から大きく改善した。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                40.0 dB
True Frequencies:   [440. 445. 450.] Hz
True Amplitudes:    [0.70398135 0.8183088  1.27613856]
True Phases:        [ 1.81840048 -2.17014818  2.63627968] rad


--- Results Summary ---
Method          | Time (s) | Freq RMSE (Hz) | Amp RMSE | Phase RMSE (rad)
----------------|----------|----------------|----------|-----------------
Spectral MUSIC  | 0.614614 | 1.496398       | 0.543398 | 0.729801         
ESPRIT (LS)     | 0.264415 | 0.402776       | 0.074715 | 0.077601

実験その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の範囲でランダム

ESPRIT 法といえども、前方平均化の共分散行列（通常の共分散行列）だけでは、5 Hz 間隔の周波数を正確に推定することはできなかった。実際、440 Hz と 450 Hz は概ね推定できているものの、その間に挟まれた 445 Hz は推定できていない。前方後方平均化では、たとえ 20 dB でも概ね周波数は推定できており（とはいえ 445 Hz の推定誤差は大きいが）、推定精度が改善された。すなわち耐雑音性能の向上が確認できた。もちろん、真の振幅と位相の値によっては推定に失敗することもあるが、その失敗頻度は前方平均化の場合よりも小さい（繰り返し実験すると体感ではっきりと分かるぐらい）。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                20.0 dB
True Frequencies:   [440. 445. 450.] Hz
True Amplitudes:    [1.13594625 0.7704213  0.56046685]
True Phases:        [-1.41625345 -1.01210672  1.42849417] rad


--- Running ESPRIT (LS) ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  Solver: LSEspritSolver
Elapsed Time: 0.2862 seconds

--- Estimation Results ---
Est Frequencies: [440.65161003 442.32703136 450.59822304] Hz
Est Amplitudes:  [1.93169584 1.54954816 0.46822731]
Est Phases:      [-1.99577992  0.0358281   1.50267909] rad

--- Estimation Errors ---
Freq Errors:  [ 0.65161003 -2.67296864  0.59822304] Hz
Amp Errors:   [ 0.79574959  0.77912686 -0.09223954]
Phase Errors: [-0.57952647  1.04793482  0.07418492] rad


--- Running ESPRIT FB (LS) ---
Analyzer Parameters:
  Subspace Ratio: 0.3333
  Solver: LSEspritSolver
Elapsed Time: 0.3153 seconds

--- Estimation Results ---
Est Frequencies: [439.81928947 444.61215622 450.26492115] Hz
Est Amplitudes:  [1.10013963 0.7397889  0.50147595]
Est Phases:      [-1.38948566 -1.01215115  1.37926545] rad

--- Estimation Errors ---
Freq Errors:  [-0.18071053 -0.38784378  0.26492115] Hz
Amp Errors:   [-0.03580662 -0.0306324  -0.0589909 ]
Phase Errors: [ 2.67677971e-02 -4.44275827e-05 -4.92287190e-02] rad


--- Results Summary ---
Method               | Time (s) | Freq RMSE (Hz) | Amp RMSE | Phase RMSE (rad)
---------------------|----------|----------------|----------|-----------------
ESPRIT (LS)          | 0.286208 | 1.625549       | 0.645178 | 0.692705        
ESPRIT FB (LS)       | 0.315255 | 0.290553       | 0.043590 | 0.032352

実験その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

実験結果を以下に示す。問題なく推定できている。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                20.0 dB
True Frequencies:   [440. 450. 460.] Hz
True Amplitudes:    [1.03063261 1.20461868 0.64799218]
True Phases:        [-1.24657106 -2.09448393  0.61981189] rad

--- Results Summary ---
Method               | Time (s) | Freq RMSE (Hz) | Amp RMSE | Phase RMSE (rad)
---------------------|----------|----------------|----------|-----------------
ESPRIT (LS)          | 0.279197 | 0.320150       | 0.031997 | 0.116930        
ESPRIT FB (LS)       | 0.313936 | 0.163312       | 0.016757 | 0.059470

実験その4

前方後方平均化の効果をさらに確認するため、SNR を10 dB まで下げ、他の条件は実験その3と揃えて実験する。

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

実験結果を示す。前方後方平均化により、周波数の推定精度に改善が見られ、雑音に対する頑健性が示された。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    100 ms
SNR:                10.0 dB
True Frequencies:   [440. 450. 460.] Hz
True Amplitudes:    [1.18763578 0.69405842 0.53693511]
True Phases:        [1.96869805 2.40032009 1.04519625] rad


--- Results Summary ---
Method               | Time (s) | Freq RMSE (Hz) | Amp RMSE | Phase RMSE (rad)
---------------------|----------|----------------|----------|-----------------
ESPRIT (LS)          | 0.269526 | 0.337318       | 0.011295 | 0.143783        
ESPRIT FB (LS)       | 0.317221 | 0.175467       | 0.014797 | 0.064985

実験その5

前方後方平均化の効果をさらに確認するため、信号長を 0.08 sec まで短くし、他の条件は実験その4と揃えて実験する。

信号: 3つの近接した周波数 (440 Hz, 450 Hz, 460 Hz) を持つ正弦波の和
サンプリング周波数: 44100 Hz
信号長: 0.08 sec
各正弦波の振幅値: 0.5〜1.5の範囲でランダム
SNR: 10 dB

実験結果を以下に示す。前方後方平均化による推定精度の改善はもはや明らかで、短時間信号での優位性も確認できた。

--- Experiment Setup ---
Sampling Frequency: 44100.0 Hz
Signal Duration:    80 ms
SNR:                10.0 dB
True Frequencies:   [440. 450. 460.] Hz
True Amplitudes:    [1.4964774  0.88632881 1.14531102]
True Phases:        [ 1.8977685  -1.83208434 -2.07535027] rad


--- Results Summary ---
Method               | Time (s) | Freq RMSE (Hz) | Amp RMSE | Phase RMSE (rad)
---------------------|----------|----------------|----------|-----------------
ESPRIT (LS)          | 0.163163 | 4.680445       | 0.244014 | 0.726905        
ESPRIT FB (LS)       | 0.181734 | 0.393764       | 0.040396 | 0.105671

おまけ

ESPRIT 法と MUSIC 法の実行時間を比較する。以下のコード片（抜粋）を使って計測した。サンプリング周波数は 44100 Hz 、信号長は 100 ms 、正弦波の個数は 3 とした。また Python のバージョンは Ubuntu が 3.12 で MacBook Air が 3.13 である。Numpy のバージョンは 2.3.2、Scipy のバージョンは 1.16.1 である。ちなみに Numpy/Scipy にはデフォルトの Open BLAS を使用している。MUSIC 法は spectral-MUSIC を比較対象にする。

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

実行時間は以下の通りである（MUSIC 法は前回記事の結果を引用）。Ubuntu と MacBook Air のどちらにおいても、ESPRIT 法の実行時間がより短い結果となった。

Ubuntu 24.04 (Intel Core i9-9900K)
- ESPRIT法：0.3838 秒
- MUSIC法（グリッド点数 4096）：0.49242 秒
- MUSIC法（グリッド点数 8192）：0.5888 秒
- MUSIC法（グリッド点数 16384）：0.73899 秒
- MUSIC法（グリッド点数 32768）：0.73899 秒
MacBook Air (M1)
- ESPRIT法：0.9017 秒
- MUSIC法（グリッド点数 4096）：1.0125 秒
- MUSIC法（グリッド点数 8192）：1.1289 秒
- MUSIC法（グリッド点数 16384）：1.4161 秒
- MUSIC法（グリッド点数 16384）：1.9947 秒

考察

実験結果から、ESPRIT 法がMUSIC 法と同様に、FFT の分解能限界を超える超解像性能を持つことが確認できた（実験その1, その3）。共分散行列の推定における前方後方平均化の導入により、周波数の推定精度が大きく改善されることも確認できた（実験その2, その4, その5）。 ESPRIT 法と MUSIC 法は、同じ部分空間法に分類されるが、以下の点で対照的である。

計算方法: MUSIC 法はスペクトル探索によるピーク検出を行うのに対し、ESPRIT 法は回転作用素の固有値計算によって直接パラメタを求める。一般に、探索グリッドの密度に依存しない ESPRIT 法の方が計算コスト面で有利とされることが多い。実際、「おまけ」の実験において、ESPRIT 法の実行時間は MUSIC 法よりも短時間であった．
拡張性: ESPRIT 法は、固有値の絶対値を評価することで、減衰率も同時に推定する減衰正弦波モデルへと自然に拡張できる。これについては、今後の記事で詳しく紹介する予定である。

一方で、両手法に共通する課題として、信号に含まれる正弦波の数（モデル次数）を事前に知る必要がある点が挙げられる。

おわりに

今回の実装と実験を通じて、ESPRIT 法の素晴らしさを実感した。次回以降の記事では（次回とは言ってない）、この ESPRIT 法を拡張し、減衰正弦波のパラメタ推定へと応用する。ESPRIT 法における計算量のボトルネックは固有値分解および特異値分解であり、その計算量を削減するために多くの手法が提案されている。本記事で名前しか紹介していない Unitary ESPRIT 法は、推定対象が複素正弦波信号の場合に ESPRIT 法の計算量がさらに削減される手法である。そのほか高速な ESPRIT 法実現のための Nyström-based ESPRIT や FFT-ESPRIT は全てリポジトリに実装済なので、早く動かしたいひとはそちらをどうぞ。

参考文献

[1] R. Roy and T. Kailath, "ESPRIT-estimation of signal parameters via rotational invariance techniques," in IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 7, pp. 984-995, 1989. https://ieeexplore.ieee.org/document/32276

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

[3] Roland Badeau, Remy Boyer, and Bertrand David, "EDS parametric modeling and tracking of audio signals," Proc. of the 5th International Conference on Digital Audio Effects (DAFx), 2002. DAFx Paper Archive - EDS Parametric Modeling and Tracking of Audio Signals

[4] Roland Badeau, Bertrand David, and Gael Richard, "A new perturbation analysis for signal enumeration in rotational invariance techniques," in IEEE Transactions on Signal Processing, vol. 54, no. 2, pp. 450–458, 2006. https://dl.acm.org/doi/abs/10.1109/TSP.2005.861899

[5] George P. Kafentzis, "On the Parameter Estimation of Sinusoidal Models for Speech and Audio Signals, " https://arxiv.org/abs/2401.01255

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

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

2025-09-29

Intel MKL の Numpy / Scipy で高速化したいとき（pip install）

python

以下のリポジトリより，

github.com

このコマンドでインストールする．

pip install numpy scipy --extra-index-url https://urob.github.io/numpy-mkl

体感で倍速になった気がする（固有値分解，特異値分解）．

2025-09-20

pandocを使って数式混じりのmarkdownをpdfに変換する方法（Mac）

LaTeX

Mac でのログ。変換対象のmarkdownはlatexによる数式が混じった文章である（github）。

pandoc はbrewでインストールする
Eisvogel のページからテンプレートをダウンロード https://github.com/Wandmalfarbe/pandoc-latex-template/releases
テンプレートを解凍したフォルダにある eisvogel.latex を、~/.pandoc/templates フォルダ（ない場合は作る）に置く
pdflatexが必要なので brew install basictex をする
ここを参考に、sudo tlmgr install 経由でたくさんパッケージをインストールする
以下のコマンドでmarkdown（例ではyour_text.md）をpdfに変換する（例では-o your_text.pdf）

pandoc your_text.md \
    -o your_text.pdf \
    --from markdown+tex_math_dollars \
    --template eisvogel \
    --syntax-highlighting=idiomatic

説明

--from markdown+tex_math_dollars: 入力形式を指定している。markdown に加えて +tex_math_dollars とすることで、$...$ (インライン) と $$...$$ (ディスプレイ) 形式の数式を正しく解釈するように Pandoc に指示
--template eisvogel: 見た目をいい感じにするテンプレート Eisvogel を使う
--syntax-highlighting=idiomatic: コードブロックを綺麗にハイライト表示する効果

注意

$$...$$ 形式で囲まれた箇所では、\begin{align} 〜 \end{align} や \begin{align*} 〜 \end{align*} だとpandocによる変換時にエラーが出る。\begin{aligned} 〜 \end{aligned} ならエラーは出ない。

備忘録

高容量カーネルHopfieldネットワークの学習メカニズムに関するプレプリントを公開

高容量カーネルHopfieldネットワークの性能を定量評価したプレプリントを公開

日本語5母音の「声道断面積関数」をPythonで可視化してみた

はじめに

実装

結果

おわりに

Mac で GTAGS for Python & Emacs

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

はじめに

ESPRIT 法の理論的背景

問題設定（信号モデル）と部分空間

回転不変性について

部分空間への拡張

回転不変性を利用した周波数推定

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

実装

実験

考察

おわりに

参考文献

Intel MKL の Numpy / Scipy で高速化したいとき（pip install）

pandocを使って数式混じりのmarkdownをpdfに変換する方法（Mac）