【Rで瞬時周波数推定】意外と単純 Teager–Kaiser Energy Method：ヒルベルト変換じゃないよ

振動する時系列信号を解析するとき、私たちはしばしば「この信号は、いまどれくらいの強さで振動しているのか」「振動の速さは時間とともにどう変わっているのか」といった問いを立てます。Teager–Kaiser energy operator（TKEO）は、まさにそのような問いに対して、答えを与える方法（非線形演算子）の一つです。
　TKEO の特徴は、フーリエ変換やヒルベルト変換のような「大域的な解析」を行わずに、ごく近傍のサンプルだけを使って、信号のエネルギー的な振る舞いを捉えようとする点にあります。この発想は、単なる数値処理の工夫ではなく、信号がどのように生成されているかという物理的・生理学的な視点から生まれました。
　その源流は、1980年代初頭の音声研究にさかのぼります。音声生成の理論を研究していた Huseyin M. Teager は、発声における気流や声帯振動が、線形モデルでは十分に説明できないことに注目しました。彼の研究では、「音声は単なる波形ではなく、エネルギーの流れとして理解すべきである」という考え方が強調されており、これが後の TKEO の思想的背景になります。　　　この考えを、信号処理の立場から明確な演算子として整理したのが James F. Kaiser です。Kaiser は、連続時間の振動信号に対して、振幅と瞬時周波数の積に比例する量が「物理的なエネルギー」を表すことに注目し、微分演算を用いた非線形エネルギー演算子を導入しました。
　さらに 1990年、Kaiser はこの考え方を離散時間信号に落とし込み、連続する3点だけから、

$\displaystyle{ \Psi[x[n]] = x^2[n] - x[n-1]\,x[n+1] }$

の形でエネルギーを計算する簡潔な差分アルゴリズムを提示します。これが、現在広く知られている Teager–Kaiser energy operator（TKEO）です。この差分形式の登場によって、TKEO は理論的な概念にとどまらず、実際の信号解析で使える実用的な手法として普及しました。
　その後、1990年代には TKEO を基礎として、瞬時振幅や瞬時周波数を分離推定する Energy Separation Algorithm（ESA）が提案され、音声信号だけでなく、生体信号や機械振動など、振幅変調 (AM)-周波数変調 (FM)型の振動を含む幅広い分野で応用されるようになります。
　今回はこのTKEOを使った瞬時周波数、瞬時振幅の推定を考えてみます。

1. 基本アイデア

　TKEO を理解するために、まず、基本的なコサイン波信号を考えます。

$\displaystyle{ x(t) = A \cos(2\pi f\, t) }$

この信号に対し、Teager–Kaiser energy operator（TKEO）は次で定義されます。

$\displaystyle{ \Psi[x(t)] = \left(\frac{dx}{dt}\right)^2 - x(t)\frac{d^2x}{dt^2} }$

この式に

$\displaystyle{ x(t) = A \cos(2\pi f\, t) }$

を代入すると、

$\displaystyle{ \Psi[x(t)] = A^2 (2\pi f)^2= 4\pi^2 A^2 f^2 }$

となります。この式を $f$ について解くと、

$\displaystyle{ f = \frac{1}{2\pi}\sqrt{\frac{\Psi[x(t)]}{A^2}} }$

となり、振幅 $A$ が既知でれば、周波数を求めることができます。
　実際の時系列解析では $A$ も未知なので、 $A$ を消すアイデアが必要です。そこで、信号の時間微分

$\displaystyle{ y(t)=\frac{dx(t)}{dt}=-2\pi f A \sin(2\pi f t) }$

を考えます。
　これに対する TKEO は

$\displaystyle{ \Psi[y(t)] = 4\pi^4 A^2 f^4 }$

となります。これを用い元信号と微分信号の比をとれば、 $A$ が消えて、

$\displaystyle{ \frac{\Psi[y(t)]}{\Psi[x(t)]} = \frac{4\pi^4 A^2 f^4}{4\pi^2 A^2 f^2} = 4\pi^2 f^2 }$

となります。したがって、

$\displaystyle{ f = \frac{1}{2\pi} \sqrt{\frac{\Psi[\dot{x}(t)]}{\Psi[x(t)]}} }$

を計算すれば、周波数 $f$ を求められます。これが Teager–Kaiser Energy Method による瞬時周波数推定の基本式です。
　さらに、瞬時振幅も同様に求めることができます。先ほど求めた関係式

$\displaystyle{ \Psi[x(t)] = 4\pi^2 A^2 f^2 }$

に、上で得た周波数 $f$ の表式を代入すると、

$\displaystyle{ A^2 = \frac{\Psi[x(t)]}{4\pi^2 f^2} }$

となります。ここで $f$ を消去すると、

$\displaystyle{ A = \frac{\Psi[x(t)]}{\sqrt{\Psi[\dot{x}(t)]}} }$

が得られます。この式により、TKEO を用いて瞬時周波数だけでなく、瞬時振幅も同時に推定できることが分かります。

2. 振幅変調 (AM)–周波数変調 (FM)信号への拡張

　ここまでは、振幅 $A$ 、周波数 $f$ が一定の正弦波を仮定してきました。しかし、現実の信号では、

振幅が時間とともに変化する（AM：amplitude modulation）
周波数も時間とともに変化する（FM：frequency modulation）

という状況が一般的です。そこで、次のような AM–FM 信号を考えます。

$\displaystyle{ x(t) = a(t)\cos\bigl(\phi(t)\bigr), \qquad f(t) = \frac{1}{2\pi}\frac{d\phi(t)}{dt} }$

ここで、

$a(t)$ ：ゆっくり変化する振幅
$\phi(t)$ ：位相
$f(t)$ ：瞬時周波数（Hz）

とします。

■ [補足説明] 周波数が変化する信号と「瞬時周波数」

　周波数が時間依存する信号について、上の式の書き方に違和感を感じる学生がいるかもしれませんので、補足の説明をしておきます。
　周波数が時間とともに変化する信号を記述するとき、次のように書いたら良いのではないかと考えるかもしれません。

$\displaystyle{ x(t) = a(t)\cos\bigl(2\pi f(t)\, t\bigr) }$

一見すると「周波数が $f(t)$ に従っているコサイン波」を表しているように見えます。しかし、この書き方では、瞬時周波数は $f(t)$ になりません。
　周波数 とは本来「位相がどれだけの速さで回転しているか」を表す量です。したがって、瞬時周波数は 位相の時間微分として定義されます。ところが、

$\displaystyle{ \phi(t) = 2\pi f(t)\, t }$

とおいて微分すると，

$\displaystyle{ \frac{d\phi(t)}{dt} = 2\pi\Bigl[f(t) + t\,\frac{df(t)}{dt}\Bigr] }$

となります。つまり、このときの瞬時周波数は

$\displaystyle{ f_{\mathrm{inst}}(t) = \frac{1}{2\pi}\frac{d\phi(t)}{dt} = f(t) + t\,\frac{df(t)}{dt} }$

であり、 $f(t)$ そのものではありません。 $f(t)$ が時間とともに変化する場合、 $t\,df(t)/dt$ の項が無視できなくなり、位相の進み方は直感と大きくずれてしまいます。
　この問題を避けるためには、考え方を逆にします。つまり、「周波数を先に決めるのではなく」、「位相を基本量として定義」します。
　そのため、周波数が変化する信号は次の形で書く必要があります。

$\displaystyle{ x(t) = a(t)\cos\bigl(\phi(t)\bigr) }$

そして、瞬時周波数は定義として

$\displaystyle{ f(t) = \frac{1}{2\pi}\frac{d\phi(t)}{dt} }$

で与えられます。

2.1 微分の計算

　話を本題に戻して、信号

$\displaystyle{ x(t)=a(t)\cos\phi(t) }$

の時間微分を計算します。積の微分より、

$\displaystyle{ \frac{dx(t)}{dt} = \dot a(t)\cos\phi(t) - 2\pi f(t)\,a(t)\sin\phi(t) }$

ここで、 $\dot a(t)$ は時間 $t$ での1階微分、 $\ddot a(t)$ は2階微分を表すとします。この表記を用いれば 2 階微分は、

$\displaystyle{ \frac{d^2x(t)}{dt^2} = \ddot a(t)\cos\phi(t) - 4\pi f(t)\dot a(t)\sin\phi(t) \\ \qquad - 2\pi a(t)\dot f(t)\sin\phi(t) - (2\pi f(t))^2 a(t)\cos\phi(t) }$

となります。
　ここで AM–FM モデルの本質的な仮定を用います。

$a(t)$ は位相 $\phi(t)$ に比べて 十分ゆっくり変化
すなわち $\dot a(t),\ \ddot a(t),\ \dot f(t)$ は $2\pi f(t), \ a(t)$ に比べて小さい

このとき、支配的な項だけを残すと、

$\displaystyle{ \frac{dx(t)}{dt} \approx -2\pi f(t)\,a(t)\sin\phi(t) }$

$\displaystyle{ \frac{d^2x(t)}{dt^2} \approx -(2\pi f(t))^2 a(t)\cos\phi(t) }$

となります。

2.3 TKEO の計算

　上の計算結果を、Teager–Kaiser エネルギー演算子

$\displaystyle{ \Psi[x(t)] = \left(\frac{dx}{dt}\right)^2 - x(t)\frac{d^2x}{dt^2} }$

に代入します。
　まず第1項は、

$\displaystyle{ \left(\frac{dx}{dt}\right)^2 \approx (2\pi)^2 a^2(t)\, f^2(t) \, \sin^2\phi(t) }$

第2項は、

$\displaystyle{ \begin{aligned} x(t)\frac{d^2x}{dt^2} &\approx a(t)\cos\phi(t)\, \bigl[-(2\pi f(t))^2 a(t)\cos\phi(t)\bigr] \\ &= -(2\pi)^2 a^2(t)\, f^2(t)\, \cos^2\phi(t) \end{aligned} }$

したがって、

$\displaystyle{ \Psi[x(t)] \approx a^2(t)\, \dot \phi^2(t) \bigl[\sin^2\phi(t)+\cos^2\phi(t)\bigr] }$

ここで、

$\displaystyle{ \sin^2\phi(t)+\cos^2\phi(t)=1 }$

より、

$\displaystyle{ \Psi[x(t)] \approx 4\pi^2\,a^2(t)\,f^2(t) }$

が得られます。

2.4 振幅 $a(t)$ を消して瞬時周波数 $f(t)$ を推定する

　ここまでで、AM–FM 信号

$\displaystyle{ x(t)=a(t)\, \cos\phi(t), \qquad f(t)=\frac{1}{2\pi}\frac{d\phi(t)}{dt} }$

に対して、TKEO が近似的に

$\displaystyle{ \Psi[x(t)] \approx 4\pi^2\,a^2(t)\,f^2(t) }$

を与えることが分かりました。しかしこの式だけでは、 $a(t)$ も未知なので $f(t)$ を直接求められません。そこで、微分信号を併用して $a(t)$ を消します。

(1) 微分信号を作る

　信号の時間微分を

$\displaystyle{ y(t)=\frac{dx(t)}{dt} }$

とします。AM–FM の「ゆっくり変化する」仮定の下では、

$\displaystyle{ y(t)=\frac{dx(t)}{dt}\approx -2\pi f(t)\,a(t)\sin\phi(t) }$

となりました。

(2) 微分信号の TKEO を計算する

　上式は「振幅が $2\pi \, f(t)\, a(t)$ 、位相が $\phi(t)-\pi/2$ のコサイン波」とみなせます。したがって、先ほどと同じ議論をそのまま適用でき、

$\displaystyle{ \Psi[y(t)] \approx 4\pi^2\,\{2\pi \, f(t)\, a(t)\}^2\,f^2(t) }$

となります。右辺を整理すると、

$\displaystyle{ \Psi[y(t)] \approx 16\pi^4\,a^2(t)\,f^4(t) }$

が得られます。

(3) 比をとって $a(t)$ を消す

　ここで、 $\Psi[y(t)]$ を $\Psi[x(t)]$ で割ります。

$\displaystyle{ \frac{\Psi[y(t)]}{\Psi[x(t)]} \approx \frac{16\pi^4\,a^2(t)\,f^4(t)}{4\pi^2\,a^2(t)\,f^2(t)} }$

分子・分母の共通因子 $a^ 2(t)$ と $f^ 2(t)$ が消えて、

$\displaystyle{ \frac{\Psi[y(t)]}{\Psi[x(t)]} \approx 4\pi^2\,f^2(t) }$

となります。

(4) 瞬時周波数 $f(t)$ の推定

　上式を $f(t)$ について解けば、

$\displaystyle{ f(t)\approx \frac{1}{2\pi}\sqrt{\frac{\Psi[y(t)]}{\Psi[x(t)]}} }$

が得られます。これが Teager–Kaiser Energy Method による瞬時周波数推定の基本式です。

(5) 瞬時振幅 $a(t)$ の推定

　瞬時周波数 $f(t)$ が求まれば、瞬時振幅も TKEO を用いて推定できます。AM–FM 信号に対しては、前節で示したように、

$\displaystyle{ \Psi[x(t)] \approx 4\pi^2\,a^2(t)\,f^2(t) }$

が成り立ちます。したがって、この式を $a(t)$ について解くと、

$\displaystyle{ a(t) \approx \frac{1}{2\pi} \frac{\sqrt{\Psi[x(t)]}}{f(t)} }$

が得られます。さらに、前節で導いた瞬時周波数の推定式

$\displaystyle{ f(t)\approx \frac{1}{2\pi}\sqrt{\frac{\Psi[y(t)]}{\Psi[x(t)]}} }$

を用いて $f(t)$ を消去すると、

$\displaystyle{ a(t) \approx \frac{\Psi[x(t)]}{\sqrt{\Psi[y(t)]}} }$

となります。この結果から、瞬時振幅も推定できます。

2.5 離散時系列（サンプル列）への適用

　実際の時系列解析では、信号は離散時間サンプリングされた列

$\displaystyle{ x[n] = x(n\Delta t) }$

として与えられます。ここで $\Delta t$ はサンプリング間隔であり、サンプリング周波数は $f_s=1/\Delta t$ です。このとき、瞬時周波数も Hz（1秒あたり）で推定したいので、式の中に $\Delta t$ を明示的に入れておきます。

(1) TKEO の離散版

　連続時間の TKEO は

$\displaystyle{ \Psi[x(t)] = \left(\frac{dx}{dt}\right)^2 - x(t)\frac{d^2x}{dt^2} }$

でした。この式の離散時間版は

$\displaystyle{ \Psi[x[n]] = x^2[n] - x[n-1]x[n+1] }$

です。この結果は、連続時間版の単純な差分化では導けません。以下では、なぜこの式が良いのかについて簡単に説明します。
　まず、周波数 $f$ （Hz）の正弦波を、サンプリング間隔 $\Delta t$ で観測した離散信号

$\displaystyle{ x[n] = A \cos\bigl(2\pi f\, n\Delta t\bigr) }$

を考えます。これは、連続時間信号 $x(t)=A\cos(2\pi f t)$ をサンプリングしたものです。
　この信号を、離散時間 TKEO の式として導入された

$\displaystyle{ \Psi[x[n]] = x^2[n] - x[n-1]x[n+1] }$

に代入してみます。
　計算に必要な部分を整理していくと、

$\displaystyle{ x[n-1] = A\cos\bigl(2\pi f (n-1)\Delta t\bigr) }$

$\displaystyle{ x[n+1] = A\cos\bigl(2\pi f (n+1)\Delta t\bigr) }$

なので、

$\displaystyle{ x[n-1]x[n+1] = A^2 \cos\bigl(2\pi f (n-1)\Delta t\bigr) \cos\bigl(2\pi f (n+1)\Delta t\bigr) }$

となります。ここで、三角関数の積和の公式

$\displaystyle{ \cos(\alpha-\beta)\cos(\alpha+\beta) = \cos^2\alpha - \sin^2\beta }$

を用いると、

$\displaystyle{ x[n-1]x[n+1] = A^2 \bigl\{ \cos^2(2\pi f n\Delta t) - \sin^2(2\pi f\Delta t) \bigr\} }$

が得られます。
　さらに、

$\displaystyle{ x^2[n] = A^2\cos^2(2\pi f n\Delta t) }$

との差をとると

$\displaystyle{ \Psi[x[n]] = x^2[n] - x[n-1]x[n+1] = A^2 \sin^2(2\pi f\Delta t) }$

となります。
　この結果から分かるように、離散時間の TKEO は

振幅の二乗 $A^ 2$
周波数 $f$
サンプリング間隔 $\Delta t$

のみに依存し、時間インデックス $n$ には依存しません。言い換えると、素直な離散化で差分の表現をえると、 $n$ 依存性が残ってしまうのです。

(2) 離散時系列の「微分信号」の計算

　連続時間では $y(t)=\dot x(t)$ を用いましたが、離散時間では差分を使います。最も簡単には、

$\displaystyle{ y[n] = x[n] - x[n-1] }$

と定義します（中心差分を使うことも可能ですが、ここでは実装が単純な形にします）。　差分信号にも同様に TKEO を適用し、

$\displaystyle{ \Psi[y[n]] = y^2[n] - y[n-1]\,y[n+1] }$

とします。

(3) 瞬時周波数 $f$ の推定

　連続時間で得た基本式

$\displaystyle{ f(t) \approx \frac{1}{2\pi} \sqrt{\frac{\Psi[y(t)]}{\Psi[x(t)]}} }$

を、離散時間で Hz のまま使える形にすると、サンプリング間隔 $\Delta t$ によってスケールが変わるため、

$\displaystyle{ f[n] \approx \frac{1}{2\pi\Delta t}\, \sqrt{\frac{\Psi[y[n]]}{\Psi[x[n]]}} }$

となります。サンプリング周波数 $f_s=1/\Delta t$ を用いれば、

$\displaystyle{ f[n] \approx \frac{f_s}{2\pi}\, \sqrt{\frac{\Psi[y[n]]}{\Psi[x[n]]}} }$

となります。
　この式は実装が簡単で広く用いられていますが、離散時間における低周波近似（すなわち $\sin\Omega \approx \Omega$ 、 $\Omega$ は rad/sample）に基づいている点に注意が必要です。
　理想的な正弦波

$\displaystyle{ x[n]=A\cos(\Omega n) }$

に対しては、離散 TKEO により次の厳密な関係が成り立ちます。

$\displaystyle{ \Psi[x[n]] = A^2\sin^2\Omega }$

また、差分信号 $\displaystyle{ y=x[n]-x[n-1] }$ については、

$\displaystyle{ \frac{\Psi[y[n]]}{\Psi[x[n]]} = 4\sin^2\!\left(\frac{\Omega}{2}\right) }$

が得られます。したがって、これを逆に解くことで、1 サンプルあたりの角周波数 $\Omega[n]$ は

$\displaystyle{ \Omega[n] = 2\arcsin\!\left( \frac{1}{2} \sqrt{\frac{\Psi[y[n]]}{\Psi[x[n]]}} \right) }$

と厳密に求められます。これを Hz に変換すると、

$\displaystyle{ f[n] = \frac{f_s}{2\pi}\, \Omega[n] = \frac{f_s}{\pi}\, \arcsin\!\left( \frac{1}{2} \sqrt{\frac{\Psi[y[n]]}{\Psi[x[n]]}} \right) }$

が、離散時間に整合した厳密な推定式です。

3. Rでやってみる

　実際に Teager–Kaiser Energy operator（TKEO）法を R で試してみます。ここでは、振幅も周波数も時間変化する AM–FM 信号を人工的に作り、TKEO から

瞬時周波数 $f(t)$ （Hz）
瞬時振幅 $a(t)$

を同時に推定できることを確認します。まずは、理論どおりに動くかを見たいので、ノイズなしのデモから始めます。ヒルベルト変換と同様に、解析対象の信号は0を中心に上下にほぼ対象に振動している必要があることに注意してください。

3.1 スクリプトのポイント

　このデモスクリプトは、大きく次の 4 ステップです。

AM–FM 信号を合成 ゆっくり変化する振幅 $a(t)$ と、ゆっくり変化する瞬時周波数 $f(t)$ を設定し、位相

$\displaystyle{ \phi(t)=2\pi\int_0^t f(\tau)\,d\tau }$

を数値積分（累積和）で作って、
$\displaystyle{ x(t)=a(t)\cos\phi(t) }$
を生成します。これで真値 $a _ {\rm true}(t)$ 、 $f _ {\rm true}(t)$ も同時に得られます。
離散 TKEO を計算 離散時間サンプル列では、

$\displaystyle{ \Psi[x[n]]=x^2[n]-x[n-1]x[n+1] }$

を用います。
微分の代わりに差分

$\displaystyle{ y[n]=x[n]-x[n-1] }$

を作り、 $\Psi[y[n]]$ も同様に計算します。

周波数 $f$ を推定 理想的な正弦波では、離散 TKEO により

$\displaystyle{ \frac{\Psi[y[n]]}{\Psi[x[n]]}=4\sin^2\!\left(\frac{\Omega[n]}{2}\right) }$

が成り立つので、これを逆に解いて

$\displaystyle{ \Omega[n]=2\arcsin\!\left(\frac{1}{2}\sqrt{\frac{\Psi[y[n]]}{\Psi[x[n]]}}\right) }$

を得て、Hz に変換します（ $f_s$ はサンプリング周波数）。

$\displaystyle{ f[n]=\frac{f_s}{2\pi}\,\Omega[n] }$

振幅 $a$ を推定 離散 TKEO の関係

$\displaystyle{ \Psi[x[n]]=a^2[n]\sin^2\Omega[n] }$

から、

$\displaystyle{ a[n]=\frac{\sqrt{\Psi[x[n]]}}{|\sin\Omega[n]|} }$

により振幅を推定します。

最後に、(i) 信号 $x(t)$ 、(ii) 周波数の真値と推定値の重ね描き、(iii) 振幅の真値と推定値の重ね描きを 3 段の図として表示します。

3.2 実装上の工夫

$\displaystyle{\Psi[x[n]]}$ や比が数値的に不安定になるのを防ぐため、ごく小さな $\varepsilon$ を用いて下限を設けています。
理論上、 $\displaystyle{\Psi[y]/\Psi[x]}$ は $[0, 4$ ] に入るので、推定が破綻しないようにこの範囲にクリップしています。
端点は $n-1$ と $n+1$ が必要なため、最初と最後では $\Psi$ が計算できず、一部に NA が出ます（これは正常です）。

############################################################
# TKEO demo in R
#  - Estimate instantaneous frequency f[n] and amplitude a[n]
#  - Discrete-time formulas consistent with Psi[x[n]]=A^2 sin^2(Omega)
############################################################

# -----------------------------
# 0) Helper: discrete TKEO
# -----------------------------
tkeo <- function(x) {
  n <- length(x)
  psi <- rep(NA_real_, n)
  if (n >= 3L) {
    psi[2:(n - 1L)] <- x[2:(n - 1L)]^2 - x[1:(n - 2L)] * x[3:n]
  }
  psi
}

# -----------------------------
# 1) Synthetic AM–FM signal (noise-free)
# -----------------------------
fs <- 200
dt <- 1 / fs
Tsec <- 10
n <- Tsec * fs
t <- (0:(n - 1)) * dt

# True amplitude (AM)
a_true <- 1.0 + 0.4 * sin(2 * pi * 0.2 * t)

# True instantaneous frequency (FM) [Hz]
f_true <- 8 + 3 * sin(2 * pi * 0.1 * t)

# Phase: phi(t) = 2*pi * integral f(t) dt
phi <- 2 * pi * cumsum(f_true) * dt

# Signal
x <- a_true * cos(phi)

# -----------------------------
# 2) TKEO-based estimation (discrete, correct formulas)
# -----------------------------
# Difference signal
y <- c(NA_real_, diff(x))   # y[n] = x[n] - x[n-1]

psi_x <- tkeo(x)
psi_y <- tkeo(y)

# Safety
eps <- 1e-12
psi_x_pos <- pmax(psi_x, eps)
psi_y_pos <- pmax(psi_y, 0)

# Ratio r = Psi[y]/Psi[x] = 4 sin^2(Omega/2)  (ideal sinusoid)
r <- psi_y_pos / psi_x_pos
r <- pmin(pmax(r, 0), 4)  # theoretical range [0,4]

# Omega_hat (rad/sample), internal only (not shown in the blog text if you prefer)
Omega_hat <- 2 * asin(0.5 * sqrt(r))

# Frequency estimate in Hz
f_hat <- (fs / (2 * pi)) * Omega_hat

# Amplitude estimate (consistent with Psi[x]=a^2 sin^2(Omega))
a_hat <- sqrt(psi_x_pos) / pmax(abs(sin(Omega_hat)), eps)

# -----------------------------
# 3) Plot (requested colors)
# -----------------------------
op <- par(no.readonly = TRUE)
par(mfrow = c(3, 1), mar = c(5, 5, 3, 1), las = 1)

# (Top) signal col=4
plot(t, x, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "x(t)",
     main = "AM–FM signal")

# (Middle) frequency overlay: true col=4, estimated col=2
plot(t, f_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Frequency [Hz]",
     main = "Instantaneous frequency: true vs TKEO")
lines(t, f_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True f(t)", "TKEO estimate"),
       col = c(4, 2), lwd = 2, bty = "n")

# (Bottom) amplitude overlay: true col=4, estimated col=2
plot(t, a_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Amplitude",
     main = "Amplitude: true vs TKEO estimate")
lines(t, a_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True a(t)", "TKEO estimate"),
       col = c(4, 2), lwd = 2, bty = "n")

par(op)

　この R スクリプトを実行すると、下図のような結果が得られます。上段には生成した AM-FM 信号、中段には瞬時周波数の真値と TKEO による推定結果、下段には瞬時振幅の真値と 推定結果を重ねて表示しています。ノイズを含まない理想条件では、瞬時周波数・瞬時振幅のいずれについても、推定結果が真値とほぼ完全に一致していることが分かります。
　とくに重要なのは、振幅が時間とともに変化しているにもかかわらず、瞬時周波数の推定がその影響をほとんど受けていない点です。これは、TKEO 法では 元信号と微分信号の比を取ることで振幅成分が理論的に消去されていることを、数値的にも裏づけています。また、振幅推定についても、離散 TKEO の関係式を正しく用いることで、時間変化する包絡を高い精度で復元できていることが確認できます。
　この結果は、AM–FM 仮定が満たされる範囲において、TKEO 法が瞬時周波数と瞬時振幅を同一の枠組みのもとで同時に推定できることを示しています。次の節では、この理想的な一致が、わずかなノイズの付加によってどのように崩れるのかを確認します。

3.3 TKEO法はノイズに弱い

　上の例では、ノイズを含まない理想的な AM–FM 信号を用いていました。その場合、TKEO による瞬時周波数・振幅推定は理論どおり非常にうまくいくことが確認できました。しかし、TKEO 法には明確な弱点があります。それが「ノイズに対する脆弱性」です。
　このことを確認するために、ごく小さなノイズを加えた場合に、推定結果がどのように変化するかを見てみます。　先ほどの R スクリプトでは、信号生成部分は次のようになっていました。

# Signal
x <- a_true * cos(phi)

ここに、平均 0、標準偏差 0.005 のガウス雑音を加えます。

# Signal
x <- a_true * cos(phi) + rnorm(n,0,0.005)

この修正以外は一切変更せず、まったく同じ手順で TKEO による解析を行います。　結果は下図です。図を見ればわかるように、ノイズを含むTKEO法は、まったく役に立ちません。

　このままでは実用上、役に立ちそうにない TKEO 法ですが、どのように工夫すれば実データにも耐える手法へと改善できるのか、その具体的な解決策を次の記事で示します。

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。