レーダーのドップラー解析や振動診断の現場では、「観測できるデータがほんの数十点しかないのに、すぐ近くにある2つの周波数成分を見分けたい」という場面がよくあります。たとえば回転機械の振動から、互いに0.02 Hzしか離れていない2つのピークを切り分けたい。ところが、FFTにそのまま窓をかけてパワーを計算する古典的なピリオドグラムでは、周波数分解能はデータ長 $N$ で決まる $1/N$ で頭打ちになり、短いデータでは2つのピークが1つの山に潰れてしまいます。
この限界を破る発想が パラメトリックスペクトル推定 です。「観測した信号は、ある単純なモデル（ここでは自己回帰=ARモデル）から生成された」と仮定し、そのモデルのパラメータをデータから推定します。すると、推定したモデルが暗黙のうちに「データの外側」を合理的に補外してくれるため、短いデータでも鋭いピークを再現できます。AR法は、わずか数十点のデータから近接した2つの正弦波を分離できる「高分解能スペクトル推定」の代表格です。
この技術は実用上いたるところで使われています。応用先を挙げると、(1) レーダー・ソナーのドップラースペクトル推定 — 短い観測時間で複数のターゲット速度を分離する、(2) 音声処理のLPC（線形予測符号化） — 声道の共鳴（フォルマント）をARの極として抽出し、携帯電話の音声圧縮に使う、といった具合です。本記事では、ARモデルのパワースペクトル密度（PSD）の式を導出し、自己相関からYule-Walker方程式を立て、それをLevinson-Durbin再帰で高速に解く手順を一切省略せずに追います。最後に、近接2正弦波に対してピリオドグラムとAR法の分解能の差をPythonで実装して確かめ、次数 $p$ の選び方（AIC）と過大次数で生じるスプリアスピークまで触れます。
本記事の内容
- AR過程の直感的な意味と、PSDの式 $S(f) = \sigma^2 / |1 + \sum a_k e^{-j2\pi f k}|^2$ の導出
- 自己相関から Yule-Walker方程式 を立てる手順（省略なし）
- それを $O(p^2)$ で解く Levinson-Durbin再帰 の導出と実装
- 近接2正弦波に対する、ピリオドグラムとAR法の分解能比較（Python）
- 次数 $p$ の選び方（AIC）と、過大次数で生じるスプリアスピーク
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ノンパラメトリック法の限界とパラメトリック法の発想
スペクトル推定には、大きく分けて2つの流派があります。1つは、信号にモデルを仮定せず、観測データを直接フーリエ変換してパワーを測る ノンパラメトリック法（ピリオドグラム、Welch法など）。もう1つが、信号を生成する確率モデルを仮定し、そのパラメータを推定する パラメトリック法 です。
ノンパラメトリック法の代表であるピリオドグラムは、長さ $N$ のデータ $x[0], \dots, x[N-1]$ から
$$ \hat{S}_{\text{per}}(f) = \frac{1}{N}\left| \sum_{n=0}^{N-1} x[n]\, e^{-j2\pi f n} \right|^2 $$
として計算します。これは直感的で計算も速いのですが、致命的な弱点があります。観測窓が長さ $N$ の矩形窓であるため、スペクトルは幅 $\sim 1/N$ の主ローブを持つ「窓のスペクトル」とのたたみ込みになります。つまり、真のスペクトルがどんなに鋭いピークでも、ピリオドグラムでは幅 $1/N$ にぼやける のです。2つのピークの間隔が $1/N$ より狭ければ、両者は1つに溶けてしまいます。
ここで発想を変えます。「観測したデータは、ごく少数のパラメータを持つモデルが生み出した」と仮定したらどうでしょう。たとえば「この信号は、過去の自分の値の線形和に白色雑音を足したもの（=AR過程）だ」と仮定する。そうすれば、推定すべきはわずか $p+1$ 個のパラメータだけになり、データが短くてもパラメータは安定して決まります。そして決まったモデルは、観測区間の外まで信号がどう続くかを暗黙に予言してくれるので、実効的に「無限に長いデータを見た」かのような鋭いスペクトルが得られます。これがパラメトリック法の威力です。
では、その「単純なモデル」として最もよく使われるAR過程とは何か、次の節で具体的に見ていきます。
AR過程とは
AR（AutoRegressive、自己回帰）過程を一言で言えば、「今の値は、過去の自分の値たちの重み付き和に、小さなランダムなショックを足したもの」 です。たとえば今日の気温は、昨日・一昨日の気温にだいたい引きずられつつ、その日固有の偶然のゆらぎが乗っている、というイメージです。過去を引きずる「慣性」がモデルに組み込まれているのがポイントで、この慣性こそがスペクトルに鋭いピーク（共鳴）を生みます。
数式で書くと、次数 $p$ のAR過程 AR($p$) は
$$ \begin{equation} x[n] = -\sum_{k=1}^{p} a_k\, x[n-k] + w[n] \end{equation} $$
と表されます。ここで $a_1, \dots, a_p$ は AR係数、$w[n]$ は分散 $\sigma^2$ の白色雑音（平均0、互いに無相関）です。右辺第1項が「過去の自分への依存」、$w[n]$ が「その時刻の新しいショック（イノベーション）」です。係数の前にマイナス符号を付けているのは、後で出てくる伝達関数の分母を $A(z) = 1 + \sum a_k z^{-k}$ という素直な形にするための慣習で、深い意味はありません。
この式を移項すると、白色雑音 $w[n]$ から $x[n]$ を作る「フィルタ」の形がはっきりします。
$$ x[n] + \sum_{k=1}^{p} a_k\, x[n-k] = w[n] $$
左辺は、係数 $1, a_1, \dots, a_p$ を持つ 全極型（all-pole）IIRフィルタ の差分方程式そのものです。つまりAR過程は、白色雑音を全極フィルタに通したもの、と言い換えられます。フィルタの伝達関数や差分方程式の基礎はディジタルフィルタの基礎（FIR/IIR）を解説で扱っているので、不安な方は先に目を通しておくとよいでしょう。
白色雑音はあらゆる周波数を均等に含む「真っ白な」入力です。それを全極フィルタに通すと、フィルタの極（分母の零点）に近い周波数だけが強く増幅され、スペクトルに鋭いピークが立ちます。AR過程が共鳴を表現できるのは、この全極構造のおかげです。
この模式図がAR法の全体像です。平坦なスペクトルを持つ白色雑音 $w[n]$ を全極フィルタ $H(z)=1/A(z)$ に通すと、極が単位円に近い周波数だけが強く増幅され、出力 $x[n]$ のスペクトルに鋭いピークが現れます。つまりAR法は「どんな極を持つフィルタなら観測データを作れるか」を逆算する作業だと言えます。では、入力スペクトルがフィルタを通るとどう変わるのか、その関係を使ってAR過程のPSDを導きます。
ARパワースペクトル密度の導出
ゴールは、AR($p$) 過程のパワースペクトル密度（PSD）が
$$ \begin{equation} S(f) = \frac{\sigma^2}{\left| 1 + \sum_{k=1}^{p} a_k\, e^{-j2\pi f k} \right|^2} \end{equation} $$
という形になることを示すことです。導出の道具は、線形時不変（LTI）システムを通したときのPSDの変換則です。
LTIシステムを通したPSDの変換則
入力 $w[n]$、出力 $x[n]$、伝達関数 $H(z)$ のLTIフィルタを考えます。入力のPSDを $S_w(f)$、出力のPSDを $S_x(f)$ とすると、両者は周波数応答 $H(e^{j2\pi f})$ の振幅2乗を介して
$$ \begin{equation} S_x(f) = \left| H(e^{j2\pi f}) \right|^2 S_w(f) \end{equation} $$
で結ばれます。これはPSDの基本定理で、パワースペクトル密度（PSD）とはでも扱っています。直感的には「フィルタが各周波数の振幅を $|H|$ 倍するなら、パワー（振幅の2乗）は $|H|^2$ 倍される」ということです。
AR過程の伝達関数
AR過程の差分方程式 $x[n] + \sum_{k=1}^p a_k x[n-k] = w[n]$ の両辺をz変換します。時間シフト $x[n-k]$ のz変換が $z^{-k} X(z)$ になることを使うと
$$ X(z)\left( 1 + \sum_{k=1}^{p} a_k z^{-k} \right) = W(z) $$
となります。ここで分母多項式を
$$ A(z) = 1 + \sum_{k=1}^{p} a_k z^{-k} $$
と置くと、白色雑音 $w$ から信号 $x$ への伝達関数 $H(z) = X(z)/W(z)$ は
$$ H(z) = \frac{1}{A(z)} = \frac{1}{1 + \sum_{k=1}^{p} a_k z^{-k}} $$
です。分子が定数1で分母だけが多項式なので、これは確かに全極型フィルタです。
PSDの式を組み立てる
周波数応答は $H(z)$ に $z = e^{j2\pi f}$ を代入したものです（単位円上で評価する）。
$$ H(e^{j2\pi f}) = \frac{1}{1 + \sum_{k=1}^{p} a_k\, e^{-j2\pi f k}} $$
一方、白色雑音 $w[n]$ のPSDは全周波数で一定で、その値は分散 $\sigma^2$ に等しくなります（白色＝平坦なスペクトル）。
$$ S_w(f) = \sigma^2 $$
これらを変換則 $S_x(f) = |H(e^{j2\pi f})|^2 S_w(f)$ に代入すると
$$ S_x(f) = \left| \frac{1}{1 + \sum_{k=1}^{p} a_k e^{-j2\pi f k}} \right|^2 \sigma^2 = \frac{\sigma^2}{\left| 1 + \sum_{k=1}^{p} a_k e^{-j2\pi f k} \right|^2} $$
となり、目標の式 (2) が得られました。
この式の意味を噛みしめておきましょう。分母 $A(e^{j2\pi f}) = 1 + \sum a_k e^{-j2\pi f k}$ がある周波数 $f$ でゼロに近づくと、$S(f)$ は急激に大きくなります。$A(z) = 0$ の根が極であり、その極が単位円のすぐ内側にあると、対応する周波数で鋭いピークが立つわけです。AR法のスペクトルは「分母の谷」を「ピークの山」に反転させて作る、と覚えておくと見通しがよくなります。
後で手で追うAR(2)（自己相関 $r=[1.0, 0.6, 0.1]$ から係数 $a_1=-0.84375,\ a_2=0.40625$）の場合に、PSDと極配置が実際にどうなるかを先取りして見てみましょう。
左のPSDは $f\approx0.1$ あたりに1つの緩やかな山を持ち、右の極配置図では2つの共役な極が単位円の内側（半径約0.64）に位置しています。極の絶対値が1未満なのでこのフィルタは安定で、極が単位円に近づくほどピークは鋭くなる、という対応関係がはっきり読み取れます。
PSDの式 (2) は、AR係数 $a_k$ と雑音分散 $\sigma^2$ さえ分かれば計算できます。問題は、観測データからこの $a_k$ と $\sigma^2$ をどう推定するか。その鍵を握るのが、次に導くYule-Walker方程式です。
Yule-Walker方程式の導出
AR係数を求める最も基本的な方法が Yule-Walker（ユール・ウォーカー）方程式 です。これは、AR過程の自己相関関数が満たすべき関係式から導かれる連立1次方程式です。導出のアイデアは「差分方程式の両辺に $x[n-m]$ を掛けて期待値を取る」だけ。やってみましょう。
自己相関関数の定義
定常過程 $x[n]$ の自己相関関数を、ラグ $m$ に対して
$$ r[m] = E\big[ x[n]\, x[n-m] \big] $$
と定義します（実数信号を仮定、平均は0とします）。自己相関は「信号を $m$ だけずらした自分自身とどれだけ似ているか」を測る量で、$r[m] = r[-m]$（偶関数）という対称性を持ちます。
差分方程式に $x[n-m]$ を掛けて期待値を取る
AR過程の式 $x[n] + \sum_{k=1}^p a_k x[n-k] = w[n]$ の両辺に $x[n-m]$（$m \geq 0$）を掛けます。
$$ x[n]\,x[n-m] + \sum_{k=1}^{p} a_k\, x[n-k]\,x[n-m] = w[n]\,x[n-m] $$
両辺の期待値を取ります。期待値は線形なので和の外に出せて
$$ E\big[x[n]x[n-m]\big] + \sum_{k=1}^{p} a_k\, E\big[x[n-k]x[n-m]\big] = E\big[w[n]x[n-m]\big] $$
左辺の各期待値は自己相関に書き換えられます。$E[x[n]x[n-m]] = r[m]$、$E[x[n-k]x[n-m]] = r[m-k]$ なので
$$ \begin{equation} r[m] + \sum_{k=1}^{p} a_k\, r[m-k] = E\big[w[n]x[n-m]\big] \end{equation} $$
右辺の評価がカギ
右辺 $E[w[n]x[n-m]]$ を、$m$ の値で場合分けして評価します。ここがYule-Walker方程式の心臓部です。
$x[n-m]$ は、時刻 $n-m$ までの雑音 $w[n-m], w[n-m-1], \dots$ の重ね合わせで作られています（因果的なフィルタなので、未来の雑音は含まれません）。一方 $w[n]$ は時刻 $n$ のショックです。
- $m \geq 1$ のとき: $x[n-m]$ は $w[n]$ より過去の雑音だけからできているので、$w[n]$ とは無相関です。したがって $E[w[n]x[n-m]] = 0$。
- $m = 0$ のとき: $E[w[n]x[n]]$ を評価します。$x[n] = w[n] – \sum_k a_k x[n-k]$ を代入すると、$E[w[n]x[n]] = E[w[n]^2] – \sum_k a_k E[w[n]x[n-k]]$。第2項は上の理由で全てゼロ、第1項は雑音の分散なので $E[w[n]x[n]] = \sigma^2$。
この場合分けを式 (4) に反映させると、2つの式が得られます。
$m \geq 1$ の場合（Yule-Walker方程式）:
$$ \begin{equation} r[m] + \sum_{k=1}^{p} a_k\, r[m-k] = 0, \quad m = 1, 2, \dots, p \end{equation} $$
$m = 0$ の場合（分散の式）:
$$ \begin{equation} r[0] + \sum_{k=1}^{p} a_k\, r[k] = \sigma^2 \end{equation} $$
行列形式（正規方程式）
式 (5) を $m = 1, \dots, p$ について並べ、$r[-k] = r[k]$ の対称性を使って行列で書くと
$$ \begin{equation} \begin{bmatrix} r[0] & r[1] & \cdots & r[p-1] \\ r[1] & r[0] & \cdots & r[p-2] \\ \vdots & \vdots & \ddots & \vdots \\ r[p-1] & r[p-2] & \cdots & r[0] \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_p \end{bmatrix} = -\begin{bmatrix} r[1] \\ r[2] \\ \vdots \\ r[p] \end{bmatrix} \end{equation} $$
となります。左辺の行列 $\bm{R}$ は、各対角線上の要素がすべて等しい テプリッツ（Toeplitz）行列 であり、さらに対称行列です。この $p \times p$ の連立方程式を解けば $a_1, \dots, a_p$ が求まり、続いて式 (6) で $\sigma^2$ が求まります。
実際には真の自己相関 $r[m]$ は分からないので、観測データから推定した標本自己相関
$$ \hat{r}[m] = \frac{1}{N} \sum_{n=0}^{N-1-m} x[n]\, x[n+m], \quad m = 0, 1, \dots, p $$
を代入します。ここで $1/N$ で割る（$1/(N-m)$ ではない）のが バイアス付き推定 で、これを使うと $\bm{R}$ が必ず半正定値テプリッツ行列になり、得られるARフィルタが安定（極が単位円内）になることが保証されます。これは実用上とても重要な性質です。
式 (7) はそのまま $\bm{a} = -\bm{R}^{-1}\bm{r}$ と逆行列で解けますが、$p$ が大きいと逆行列計算は $O(p^3)$ かかります。テプリッツ構造を活かせばもっと速く解けます。それが次のLevinson-Durbin再帰です。
Levinson-Durbin再帰
Yule-Walker方程式の係数行列はテプリッツ対称行列という強い構造を持っています。Levinson-Durbin（レビンソン・ダービン）再帰は、この構造を巧みに使って、次数を $1, 2, \dots, p$ と1つずつ上げながら解を更新していくアルゴリズムです。計算量は $O(p^2)$ で、逆行列法の $O(p^3)$ より圧倒的に速く、しかもメモリも $O(p)$ で済みます。音声コーデックがリアルタイムでLPC係数を計算できるのは、この再帰のおかげです。
アイデア: 次数を1つ上げる
次数 $m$ までの解（AR係数 $a_1^{(m)}, \dots, a_m^{(m)}$ と予測誤差 $E_m$）が手元にあるとします。これを使って、次数 $m+1$ の解を「少しの計算」で作るのが目標です。
直感的には、次数 $m$ のモデルで予測しきれなかった「取りこぼし」が、$m+1$ 番目のラグの自己相関とどれだけ食い違っているかを測り、その食い違いを打ち消すように係数を1つ追加・既存係数を微調整する、という操作になります。この食い違いを表す量が 反射係数（PARCOR係数） $k_{m+1}$ です。
再帰の式
アルゴリズムは次のように進みます。記号として、次数 $m$ の係数ベクトルを $\bm{a}^{(m)} = (a_1^{(m)}, \dots, a_m^{(m)})$、予測誤差（残差分散）を $E_m$ と書きます。
初期化（次数0）:
$$ E_0 = r[0] $$
$m = 0, 1, \dots, p-1$ について繰り返す:
まず反射係数 $k_{m+1}$ を計算します。
$$ \begin{equation} k_{m+1} = -\frac{r[m+1] + \sum_{i=1}^{m} a_i^{(m)}\, r[m+1-i]}{E_m} \end{equation} $$
分子は「次数 $m$ のモデルが予測した自己相関」と「実際の $r[m+1]$」のズレに相当します。これを誤差 $E_m$ で正規化したものが反射係数です。
次に、新しい最高次係数を反射係数そのものとし、
$$ a_{m+1}^{(m+1)} = k_{m+1} $$
既存の係数を反射係数を使って更新します（$i = 1, \dots, m$）。
$$ \begin{equation} a_i^{(m+1)} = a_i^{(m)} + k_{m+1}\, a_{m+1-i}^{(m)} \end{equation} $$
この更新は「前向き予測の係数に、後ろ向き予測の係数（係数列を逆順にしたもの）を反射係数倍して足す」という形をしています。テプリッツ行列の対称性から、前向きと後ろ向きの予測問題が同じ解を共有することを利用しているのです。
最後に予測誤差を更新します。
$$ \begin{equation} E_{m+1} = E_m\, \left(1 – k_{m+1}^2\right) \end{equation} $$
反射係数が安定性を保証する
式 (11) から重要な性質が読み取れます。$E_{m+1} = E_m(1 – k_{m+1}^2)$ で、予測誤差は分散なので必ず $E_{m+1} \geq 0$。これが成り立つには
$$ |k_{m+1}| \leq 1 $$
が必要です。バイアス付き自己相関を使うと実際に常に $|k_m| < 1$ となり、このとき得られるARフィルタの全ての極が単位円の内側に入る（=フィルタが安定する）ことが証明されています。次数を上げるたびに誤差 $E_m$ が単調に減っていく様子も式 (11) から明らかで、$1-k^2 \leq 1$ だからです。
繰り返しが終わると、$\bm{a}^{(p)} = (a_1, \dots, a_p)$ がYule-Walkerの解、$E_p$ が雑音分散の推定値 $\sigma^2$ になります。これらをPSDの式 (2) に入れれば、ARスペクトルが描けます。
ここまでで、データ → 自己相関 → Levinson-Durbin → AR係数と $\sigma^2$ → PSD という一連の流れが完成しました。次は手計算で小さな例を追って、各ステップが本当に動くことを確認します。
具体例: AR(2)過程を手で追う
理論を具体的な数字で確かめましょう。AR(2)過程を考え、自己相関が $r[0] = 1.0,\ r[1] = 0.6,\ r[2] = 0.1$ と与えられたとして、Levinson-Durbin再帰でAR係数と雑音分散を求めてみます。
初期化: $E_0 = r[0] = 1.0$。
$m = 0$ → 次数1: 反射係数は式 (8) で $m=0$ のとき和の項がないので
$$ k_1 = -\frac{r[1]}{E_0} = -\frac{0.6}{1.0} = -0.6 $$
係数は $a_1^{(1)} = k_1 = -0.6$。予測誤差は
$$ E_1 = E_0(1 – k_1^2) = 1.0 \times (1 – 0.36) = 0.64 $$
$m = 1$ → 次数2: 反射係数の分子は $r[2] + a_1^{(1)} r[1] = 0.1 + (-0.6)(0.6) = 0.1 – 0.36 = -0.26$。よって
$$ k_2 = -\frac{-0.26}{0.64} = 0.40625 $$
新しい最高次係数 $a_2^{(2)} = k_2 = 0.40625$。既存係数の更新は式 (9) で $i=1$:
$$ a_1^{(2)} = a_1^{(1)} + k_2\, a_1^{(1)} = -0.6 + 0.40625 \times (-0.6) = -0.84375 $$
予測誤差は
$$ E_2 = E_1(1 – k_2^2) = 0.64 \times (1 – 0.40625^2) = 0.64 \times 0.8350 \approx 0.5344 $$
したがってこのAR(2)過程は $a_1 = -0.84375,\ a_2 = 0.40625,\ \sigma^2 \approx 0.5344$ と推定されます。反射係数はどちらも $|k| < 1$ を満たしており、得られたフィルタは安定です。
数字が小さいので、PSDの分母 $A(e^{j2\pi f}) = 1 + a_1 e^{-j2\pi f} + a_2 e^{-j4\pi f}$ が最小になる周波数の近くにピークが立つ、という構造もイメージできます。この具体例をそのままコードで再現し、さらに本題の「近接2正弦波の分離」へ進みましょう。
Pythonでの実装
ここからは実装です。まずLevinson-Durbin再帰を関数として書き、上の手計算と一致することを確かめます。続いて、近接した2つの正弦波＋雑音という短いデータに対し、ピリオドグラムとAR法のスペクトルを比較します。
日本語ラベルのために、最初にフォントを設定しておきます。
import numpy as np
import matplotlib, matplotlib.pyplot as plt
# 日本語フォント設定（環境にあるものを順に試す）
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
Levinson-Durbin再帰の実装
Yule-Walker方程式を解く中核です。入力は自己相関の配列 $r[0], \dots, r[p]$、出力はAR係数（先頭の1を含む $A = [1, a_1, \dots, a_p]$ の形）と雑音分散です。
import numpy as np
def levinson_durbin(r, p):
"""自己相関 r[0..p] からAR係数とノイズ分散を再帰で求める。
返り値: A=[1, a_1, ..., a_p], sigma2(=E_p), 反射係数のリスト"""
A = np.array([1.0]) # 次数0の係数（先頭1のみ）
E = r[0] # 予測誤差 E_0
reflection = []
for m in range(1, p + 1):
# 反射係数の分子: r[m] + Σ a_i r[m-i] （i=1..m-1）
acc = r[m] + np.dot(A[1:m], r[m-1:0:-1]) if m > 1 else r[1]
k = -acc / E # 反射係数 k_m
A = np.concatenate([A, [0.0]]) # 末尾に1要素拡張
A = A + k * A[::-1] # a_i += k * a_{m-i}（後ろ向き予測を加算）
E = E * (1.0 - k * k) # 予測誤差の更新
reflection.append(k)
return A, E, reflection
このコードの肝は A = A + k * A[::-1] の1行です。配列を逆順にした A[::-1] が「後ろ向き予測の係数」に対応し、それを反射係数 $k$ 倍して足すことで、式 (9) の更新と新係数 $a_{m+1}=k$ の追加を一度に行っています。次に、上の手計算（AR(2)、$r=[1.0, 0.6, 0.1]$）と一致するか確認します。
import numpy as np
r_demo = np.array([1.0, 0.6, 0.1])
A, sigma2, refl = levinson_durbin(r_demo, 2)
print("AR係数 A =", np.round(A, 5)) # [1, a_1, a_2]
print("ノイズ分散 sigma^2 =", round(sigma2, 5))
print("反射係数 k =", np.round(refl, 5))
実行すると A = [1. -0.84375 0.40625]、sigma^2 = 0.5344、k = [-0.6 0.40625] が得られ、先ほどの手計算と完全に一致します。反射係数がいずれも絶対値1未満なので、このARフィルタは安定です。手で追った再帰がコードでも正しく動くことが確認できました。
テストデータ: 近接した2正弦波＋雑音
いよいよ本題です。$f_1 = 0.20$ と $f_2 = 0.22$ という、わずか0.02しか離れていない2つの正弦波に白色雑音を加え、たった $N = 64$ 点だけ観測します。FFTの周波数分解能 $1/N = 0.0156$ はこの間隔0.02より細かくないので、ピリオドグラムには厳しい設定です。
import numpy as np
np.random.seed(4)
fs = 1.0 # サンプリング周波数（正規化）
N = 64 # 短いデータ長
n = np.arange(N)
f1, f2 = 0.20, 0.22 # 近接する2周波数
x = (np.sin(2*np.pi*f1*n) + np.sin(2*np.pi*f2*n)
+ 0.1*np.random.randn(N)) # 雑音を加える
print("データ長 N =", N, " FFT分解能 1/N =", round(1/N, 4),
" 周波数間隔 =", round(f2-f1, 4))
出力は N = 64 1/N = 0.0156 周波数間隔 = 0.02 です。周波数間隔0.02はFFTの分解能0.0156と同程度（むしろ分解能のほうが粗い）なので、ピリオドグラムでは2つを分離するのが難しいと予想できます。生成した時系列そのものを見てみましょう。
64点の波形は、2つの近い周波数が重なって生じる「うなり（ビート）」のため、振幅が時間とともにゆっくり変動しています。ただし観測区間が短いため、うなりの周期を1つ分すら見切れておらず、目で見て2つの周波数を読み取るのは不可能です。この同じデータに対して、2つの方法でスペクトルを推定して比べます。
ARスペクトルとピリオドグラムを計算する関数
標本自己相関を計算する関数と、AR係数からPSDを評価する関数を用意します。
import numpy as np
def autocorr_biased(x, maxlag):
"""バイアス付き標本自己相関 r[0..maxlag]（1/Nで正規化）"""
N = len(x)
r = np.zeros(maxlag + 1)
for m in range(maxlag + 1):
r[m] = np.sum(x[:N-m] * x[m:]) / N
return r
def ar_psd(A, sigma2, freqs):
"""AR係数A=[1,a_1,...,a_p]とsigma2からPSD S(f)を計算"""
p = len(A) - 1
S = np.empty(len(freqs))
for i, f in enumerate(freqs):
e = np.exp(-1j * 2*np.pi * f * np.arange(p + 1)) # [1, e^{-j2πf}, ...]
S[i] = sigma2 / np.abs(np.dot(A, e))**2 # 式(2)
return S
ar_psd はPSDの式 (2) をそのまま実装したものです。分母 $|A(e^{j2\pi f})|^2$ が小さい周波数で $S(f)$ が大きくなり、ピークが立ちます。Yule-Walkerの入力となる標本自己相関 $\hat{r}[m]$ を、実際にラグ20まで計算して可視化してみます。
$\hat{r}[0]=1$ 付近から始まり、ラグが増えるにつれて符号を変えながら振動的に減衰しています。この振動の周期がおおよそ信号の周期（中心周波数約0.21の逆数＝約4.8サンプル）に対応しており、自己相関の段階ですでに周期情報が埋め込まれていることが分かります。AR法はこの $\hat{r}[m]$ だけを入力に係数を決めます。これらを使って、次数を変えながらARスペクトルを描きます。
分解能の比較プロット
ピリオドグラム（FFT）と、次数 $p$ を変えたAR法のスペクトルを重ねて比較します。
import numpy as np
import matplotlib.pyplot as plt
freqs = np.linspace(0, 0.5, 2000)
# ピリオドグラム
X = np.fft.rfft(x)
f_per = np.fft.rfftfreq(N, 1/fs)
per = np.abs(X)**2 / N
plt.figure(figsize=(9, 5.5))
# 正規化して重ねる（ピーク=0 dB基準）
plt.plot(f_per, 10*np.log10(per/per.max()+1e-12),
color="gray", marker="o", ms=3, label="ピリオドグラム (FFT)")
for p, c in [(4, "tab:green"), (14, "tab:blue")]:
r = autocorr_biased(x, p)
A, sigma2, _ = levinson_durbin(r, p)
S = ar_psd(A, sigma2, freqs)
plt.plot(freqs, 10*np.log10(S/S.max()+1e-12), color=c, lw=2,
label=f"AR法 (p={p})")
for f0 in (f1, f2):
plt.axvline(f0, color="red", ls="--", alpha=0.6)
plt.text(f1, 2, "真の2周波数", color="red", fontsize=9)
plt.xlabel("正規化周波数 $f$")
plt.ylabel("正規化パワー [dB]")
plt.title(f"近接2正弦波の分解能比較（N={N}, 間隔={f2-f1:.2f}）")
plt.xlim(0.10, 0.35); plt.ylim(-25, 5)
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("/tmp/ar_resolution.png", dpi=150)
plt.show()
このグラフから、AR法の威力がはっきり読み取れます。ピリオドグラム（灰）は $f=0.22$ あたりに幅広い1つの主ローブを作るだけで、生のFFTグリッド上では2つの真のピークが1つに溶けてしまっています（細かいギザギザは雑音による偽の凹凸です）。低次のAR法（$p=4$、緑）も同様に1つの滑らかな山にしかなりません。ところが $p=14$ のAR法（青）になると、真の周波数 $0.20$ と $0.22$（赤破線）の位置にはっきりと分かれた2つのピークが現れ、間に浅い谷ができています。同じ64点のデータから、AR法が2正弦波を分離できているのが分かります。
次数を上げていくとどうなるか
「では次数はいくらでも上げればいいのか?」という疑問が湧きます。次数 $p$ を $2, 6, 14, 30$ と変えて、ピークの分離具合とスプリアス（偽の）ピークの出方を見てみます。
import numpy as np
import matplotlib.pyplot as plt
freqs = np.linspace(0, 0.5, 2000)
orders = [2, 6, 14, 30]
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for ax, p in zip(axes.ravel(), orders):
r = autocorr_biased(x, p)
A, sigma2, _ = levinson_durbin(r, p)
S = ar_psd(A, sigma2, freqs)
ax.plot(freqs, 10*np.log10(S/S.max()+1e-12), color="tab:blue", lw=2)
for f0 in (f1, f2):
ax.axvline(f0, color="red", ls="--", alpha=0.6)
ax.set_title(f"AR次数 p = {p}")
ax.set_xlabel("正規化周波数 $f$"); ax.set_ylabel("正規化パワー [dB]")
ax.set_xlim(0, 0.5); ax.set_ylim(-30, 5); ax.grid(alpha=0.3)
plt.suptitle("AR次数とスペクトルの分解能・スプリアス", y=1.02)
plt.tight_layout()
plt.savefig("/tmp/ar_orders.png", dpi=150)
plt.show()
4枚のグラフを見比べると、次数選びのトレードオフが見えます。$p=2$ と $p=6$ では分解能が足りず、$f=0.21$ あたりに1つの山しか立ちません。$p=14$ になると2本のピークが $0.20,\ 0.22$ にきれいに分離します。ところが $p=30$ まで上げると、本来存在しない位置に小さな偽のピーク（スプリアス）が現れ、スペクトル全体がギザギザになってきます。次数が大きすぎると、モデルが雑音の偶然のゆらぎまで「共鳴」として拾ってしまうのです。分解能を上げたいが、上げすぎると偽ピークが出る。この板挟みに目安を与える指標が、次に説明するAICです。
次数選択: AIC
最適な次数は、「モデルがデータをよく説明する（残差分散 $E_p$ が小さい）」ことと「パラメータが少ない（過適合しない）」ことのバランスで決めます。これを定量化した代表的な指標が 赤池情報量規準（AIC） です。AR次数選択では
$$ \text{AIC}(p) = N \ln E_p + 2(p + 1) $$
の形がよく使われます。第1項 $N \ln E_p$ は予測誤差が小さいほど下がる「当てはまりの良さ」、第2項 $2(p+1)$ はパラメータ数に対する「罰則」です。AICを最小にする $p$ を選びます。
import numpy as np
import matplotlib.pyplot as plt
p_max = 30
aic = []
for p in range(1, p_max + 1):
r = autocorr_biased(x, p)
A, Ep, _ = levinson_durbin(r, p)
aic.append(N * np.log(Ep) + 2 * (p + 1))
aic = np.array(aic)
p_best = np.arange(1, p_max + 1)[np.argmin(aic)]
plt.figure(figsize=(8, 4.5))
plt.plot(range(1, p_max + 1), aic, "o-", color="tab:purple")
plt.axvline(p_best, color="red", ls="--", label=f"AIC最小 p={p_best}")
plt.xlabel("AR次数 $p$"); plt.ylabel("AIC")
plt.title("AICによるAR次数の選択")
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("/tmp/ar_aic.png", dpi=150)
plt.show()
print("AICが選んだ次数 p =", p_best)
AIC曲線は、次数を上げると最初は急に下がり（当てはまりが改善）、ある次数を過ぎると罰則項が効いて緩やかに上昇に転じます。その谷底が「ちょうどよい次数」です。実際にこのデータでAICを計算すると、最小値は $p=3$ で得られ、谷底がかなり低い次数に来ていることが図からも読み取れます。これは「正弦波＋雑音」という信号に対してAICが次数を過小評価しやすいという、よく知られた性質です。純粋な正弦波はARモデルにとって「無限に鋭い極」に対応するため、低次でも残差分散がそこそこ小さくなり、罰則項が早めに勝ってしまうのです。
したがって実務では、AICの値だけを鵜呑みにせず、「AICの谷の位置」を下限の目安にしつつ、目的（いくつのピークを分離したいか）に応じて少し高めの次数も試す、という使い方が現実的です。今回なら、AICが示す低次では1つの山にしかならないので、2正弦波の分離を狙って $p=14$ 程度まで上げる、という判断になります。AICは万能ではありませんが、青天井に次数を上げないための「歯止め」として有用です。次数選択には他にも、AICの罰則を $\ln N \cdot (p+1)$ に強めたBIC（MDL）など複数の規準があり、用途に応じて使い分けます。
スペクトルとピーク周波数の数値確認
最後に、選んだ次数のARスペクトルからピーク周波数を数値で取り出し、真値とどれだけ合っているか確かめます。
import numpy as np
p = 14
r = autocorr_biased(x, p)
A, sigma2, refl = levinson_durbin(r, p)
freqs = np.linspace(0, 0.5, 4000)
S = ar_psd(A, sigma2, freqs)
# バンド[0.15,0.30]内の局所最大（ピーク）を検出
band = (freqs > 0.15) & (freqs < 0.30)
fb, Sb = freqs[band], S[band]
peaks = [fb[i] for i in range(1, len(Sb)-1)
if Sb[i] > Sb[i-1] and Sb[i] > Sb[i+1]]
print("検出したピーク周波数 =", [round(pk, 4) for pk in peaks])
print("真の周波数 =", [f1, f2])
print("反射係数の最大絶対値 =", round(max(abs(np.array(refl))), 4),
"(<1なら安定)")
出力は 検出したピーク周波数 = [0.2004, 0.2197] のように2つになり、真の周波数 $0.20,\ 0.22$ を誤差0.001程度で当てています。図でもオレンジの検出点（0.2004, 0.2197）が赤破線の真値とほぼ重なっており、目視でも一致が確認できます。反射係数の最大絶対値も 0.9358 と1未満で、推定したARフィルタが安定であることが確認できます。短い64点のデータから、間隔0.02の2正弦波を周波数誤差わずかで分離できました。これがパラメトリックAR法の実力です。
Levinson-Durbinの内部量：誤差の単調減少と安定性
ここで、Levinson-Durbin再帰が回っている間の内部量を覗いてみましょう。次数を上げるごとに予測誤差 $E_m$ がどう減るか、各段の反射係数 $|k_m|$ がどこにあるかを実測でプロットします。
左図のとおり、予測誤差 $E_m$ は次数とともに単調に減少しています。これは式 (11) の $E_{m+1}=E_m(1-k_m^2)$ で $1-k_m^2\leq1$ だからで、理論どおりの振る舞いです。右図の反射係数はすべて $|k_m|<1$（最大でも約0.94）に収まっており、得られたARフィルタが安定（全極が単位円内）であることが保証されています。
スペクトルのピークは「極の偏角」に、ピークの鋭さは「極が単位円にどれだけ近いか」に対応します。$p=14$ で推定したフィルタの極を複素平面に描いてみましょう。
14個の極のうち、$f_1=0.20$ と $f_2=0.22$ に対応する角度（赤破線）の近くに、単位円のすぐ内側へ張り付いた極のペアが見えます。これらが2本の鋭いピークを生む「信号の極」です。一方、原点近くに散らばる極は雑音を表現するための極で、単位円から離れているためスペクトルにはほとんど影響しません。極の配置を見れば、AR法が「どの極を信号、どの極を雑音と解釈したか」が一目で分かります。
超分解能はどこまで効くか：間隔を狭めて試す
「AR法は分解能が高い」と言いますが、間隔をどこまで狭めると分離できなくなるのでしょうか。$N=64$ に固定したまま、2つの周波数の間隔 $\Delta f$ を変え、各間隔で200試行ずつ「2つのピークを分離できたか（間に2 dB以上の谷ができたか）」の成功率を測ります。ピリオドグラムは生のFFTグリッド（実分解能 $1/N$）で評価します。
import numpy as np
def resolves_two(spec_db, freqs, fa, fb):
i1 = np.argmin(np.abs(freqs-fa)); i2 = np.argmin(np.abs(freqs-fb))
lo, hi = min(i1,i2), max(i1,i2)
if hi-lo < 2: return False
valley = np.min(spec_db[lo:hi+1])
pk = min(np.max(spec_db[max(0,lo-3):lo+1]), np.max(spec_db[hi:hi+4]))
return (pk - valley) >= 2.0
dfs = [0.008, 0.010, 0.012, 0.0156, 0.020, 0.025, 0.030, 0.040]
N = 64; nn = np.arange(N); fc = 0.21
freqs = np.linspace(0, 0.5, 3000)
rng = np.random.default_rng(0)
for df in dfs:
fa, fb = fc-df/2, fc+df/2; ar_ok = per_ok = 0
for _ in range(200):
xt = np.sin(2*np.pi*fa*nn)+np.sin(2*np.pi*fb*nn)+0.1*rng.standard_normal(N)
A, s, _ = levinson_durbin(autocorr_biased(xt, 14), 14)
Sar = ar_psd(A, s, freqs)
if resolves_two(10*np.log10(Sar/Sar.max()+1e-12), freqs, fa, fb): ar_ok += 1
Xp = np.fft.rfft(xt); fp = np.fft.rfftfreq(N, 1.0); Sp = np.abs(Xp)**2
if resolves_two(10*np.log10(Sp/Sp.max()+1e-12), fp, fa, fb): per_ok += 1
print(f"Δf={df:.4f} AR={ar_ok/200:.2f} ピリオドグラム={per_ok/200:.2f}")
結果は明快です。レイリー限界 $1/N=0.0156$（赤点線）より広い $\Delta f=0.025$ では、AR法が約95%の確率で分離できるのに対し、生のFFTグリッドのピリオドグラムは $\Delta f=0.03$ まで広げないと分離できません（成功率0%→100%が急に切り替わる）。AR法の曲線がピリオドグラムより左側、つまりより狭い間隔まで分離できる領域に張り出していることが、超分解能の正体です。一方で $\Delta f \lesssim 0.0156$ になると雑音の中ではAR法でも安定して分離するのは難しく、超分解能にも限界があることが分かります。
雑音への頑健性
最後に、AR法のピーク周波数推定が雑音にどれだけ強いかを調べます。$N=64$、$p=14$ のまま、加える白色雑音の標準偏差を $0.02$ から $0.8$ まで変え、各レベルで200試行の平均周波数誤差を測ります。
import numpy as np
N = 64; nn = np.arange(N); f1, f2 = 0.20, 0.22
freqs = np.linspace(0, 0.5, 3000); rng = np.random.default_rng(7)
for nl in [0.02, 0.05, 0.1, 0.2, 0.4, 0.8]:
errs = []
for _ in range(200):
xt = np.sin(2*np.pi*f1*nn)+np.sin(2*np.pi*f2*nn)+nl*rng.standard_normal(N)
A, s, _ = levinson_durbin(autocorr_biased(xt, 14), 14)
Sar = ar_psd(A, s, freqs); band = (freqs>0.15)&(freqs<0.30)
fb, sb = freqs[band], Sar[band]
pk = sorted([fb[i] for i in range(1,len(sb)-1) if sb[i]>sb[i-1] and sb[i]>sb[i+1]],
key=lambda v: -ar_psd(A,s,np.array([v]))[0])[:2]
errs.append((abs(sorted(pk)[0]-f1)+abs(sorted(pk)[1]-f2))/2 if len(pk)==2 else np.nan)
print(f"雑音σ={nl} 平均誤差={np.nanmean(errs):.4f}")
雑音が $\sigma\lesssim0.1$ の範囲では平均誤差は0.001程度とほぼ一定で、AR法は本文の設定（$\sigma=0.1$）を含む低雑音域では非常に正確です。ところが $\sigma$ が0.2を超えると誤差も標準偏差（エラーバー）も急に増大し、$\sigma=0.8$ では平均誤差が0.027まで悪化します。雑音が大きいと「信号の極」と「雑音の極」の区別が曖昧になり、ピーク位置が揺らぐためです。AR法の高分解能はあくまで十分なSNRが前提だ、という実務上の注意がここから読み取れます。
まとめ
本記事では、AR法によるパラメトリックスペクトル推定を、PSDの導出からYule-Walker方程式、Levinson-Durbin再帰、そして近接2正弦波の分離実験まで一気通貫で解説しました。
- ARモデル: 信号を「過去の自分の線形和＋白色雑音」とみなす全極型モデル。白色雑音を全極フィルタに通したものに等しい
- ARのPSD: 伝達関数 $H(z)=1/A(z)$ と変換則 $S_x=|H|^2\sigma^2$ から、$S(f)=\sigma^2/|1+\sum a_k e^{-j2\pi f k}|^2$。分母の谷がスペクトルのピークになる
- Yule-Walker方程式: 差分方程式に $x[n-m]$ を掛けて期待値を取ると、自己相関を係数とするテプリッツ連立方程式 $\bm{R}\bm{a}=-\bm{r}$ が立つ
- Levinson-Durbin再帰: テプリッツ構造を使い $O(p^2)$ で解く。反射係数 $|k_m|<1$ が安定性を保証し、$E_{m+1}=E_m(1-k_m^2)$ で誤差が単調に減る
- 分解能: ピリオドグラムは $1/N$ で頭打ちだが、AR法は短いデータでも近接2正弦波を分離できる
- 次数選択: 低すぎると分解能不足、高すぎるとスプリアス。AIC $=N\ln E_p + 2(p+1)$ の最小で折り合いをつける
AR法は、レーダーのドップラー解析や音声のLPC、振動診断など「短いデータで鋭いスペクトルがほしい」場面で威力を発揮します。一方で、モデルが合っていないとき（例: 広帯域雑音や移動平均成分が強いとき）は誤った結果を出すので、ノンパラメトリック法と使い分ける感覚が大切です。
次のステップとして、以下の記事も参考にしてください。