サイレンの音は近づくと高く、遠ざかると低く聞こえます。鳥のさえずりは一瞬のうちに音程が滑らかに変化します。こうした信号は、どの瞬間を切り取っても「そのときの周波数」が違います。ところが、私たちが最初に学ぶフーリエ変換は、信号全体をひとつの周波数スペクトルに押し込めてしまいます。「いま何Hzで鳴っているか」という、時々刻々と変わる周波数を取り出したいときには、フーリエ変換だけでは足りません。
この「瞬間ごとの周波数」を数学的にきちんと定義するための道具が、解析信号（analytic signal） と 瞬時周波数（instantaneous frequency） です。実数の信号 $x(t)$ にヒルベルト変換を組み合わせて複素信号 $z(t)$ を作ると、その大きさが「いまの振幅」、偏角の傾きが「いまの周波数」を表すようになります。
瞬時周波数は幅広い場面で使われます。たとえばレーダーやソナーの チャープ信号（周波数を掃引したパルス）の設計と解析、通信の FM復調、機械の振動診断における回転数の追跡、地震波や生体信号（心電図・脳波）の時間周波数解析などです。EMD（経験的モード分解）やヒルベルト・ホァン変換の中核でもあります。本記事では、解析信号のスペクトルがなぜ正の周波数だけになるのかを丁寧に導出し、瞬時振幅・瞬時位相・瞬時周波数を定義したうえで、scipy.signal.hilbert を使って設計どおりの掃引周波数が復元できることを確かめます。
本記事の内容
- 解析信号 $z(t) = x(t) + j\,\mathcal{H}[x](t)$ の定義と直感
- 解析信号のスペクトルが正周波数のみになることの導出
- 瞬時振幅・瞬時位相・瞬時周波数の定義と意味
- 単一成分でないと瞬時周波数が破綻すること（Bedrosianの制約）
- Pythonでチャープ信号・AM+FM信号に適用し、掃引周波数を復元する実装
- 位相アンラップがなぜ必要かの実演
全体像を先に眺めておきましょう。測定できるのは実数の $x(t)=\cos$ だけで、その値はつねに $\pm1$ の間を振動するため「いまの振幅・周波数」を一目では読み取れません。そこでヒルベルト変換で虚部 $j\sin$ を補い、複素信号 $z=a\,e^{j\phi}$（回転フェーザ）へ”持ち上げる”のが解析信号の発想です。この矢印の長さが瞬時振幅、回る速さが瞬時周波数に対応します。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
解析信号とは
まず、なぜわざわざ複素信号を作るのか、という素朴な動機から入りましょう。実数の余弦波 $x(t) = \cos(2\pi f_0 t)$ を考えます。この信号の「振幅」を知りたいとき、瞬間値 $x(t)$ をそのまま見ても意味がありません。$x(t)$ は $+1$ と $-1$ の間を振動しているので、ある瞬間の値だけでは「振幅は1だ」とは言えないのです。振幅を知るには、少なくとも1周期ぶんを眺める必要があります。
ところがオイラーの公式を思い出すと、$\cos(2\pi f_0 t)$ は複素指数 $e^{j 2\pi f_0 t}$ の実部です。もし信号そのものが $e^{j 2\pi f_0 t}$ だったら、その大きさは常に $|e^{j 2\pi f_0 t}| = 1$ で、偏角は $2\pi f_0 t$ です。つまり複素表現にできれば、「大きさ＝振幅」「偏角の傾き＝周波数」が瞬間ごとに読み取れます。回転する矢印（フェーザ）を思い浮かべてください。矢印の長さがその瞬間の強さ、矢印が回る速さがその瞬間の周波数です。
問題は、私たちが測定できるのは実数の $x(t) = \cos(2\pi f_0 t)$ だけで、虚部 $\sin(2\pi f_0 t)$ は与えられていない点です。そこで、実部だけから「本来あるべき虚部」を作り出すのがヒルベルト変換の役目です。 余弦を正弦に、つまり位相を $90^\circ$ 遅らせた信号を虚部として補ってやれば、実信号を回転フェーザに”持ち上げる”ことができます。この持ち上げた複素信号を解析信号と呼びます。
解析信号 $z(t)$ は、実信号 $x(t)$ とそのヒルベルト変換 $\mathcal{H}[x](t)$ を使って
$$ z(t) = x(t) + j\,\mathcal{H}[x](t) $$
と定義されます。実部が元の信号、虚部がヒルベルト変換です。$x(t) = \cos(2\pi f_0 t)$ なら、後で見るように $\mathcal{H}[x](t) = \sin(2\pi f_0 t)$ となり、$z(t) = \cos(2\pi f_0 t) + j\sin(2\pi f_0 t) = e^{j 2\pi f_0 t}$ ときれいな回転フェーザになります。
ここで自然な疑問が生まれます。ヒルベルト変換を虚部に足すと、なぜちょうど回転フェーザ、すなわち「負の周波数が消えた信号」になるのでしょうか。次の節でヒルベルト変換の周波数特性を確認し、そこから解析信号のスペクトルを導きます。
ヒルベルト変換の復習
解析信号の性質を導くには、ヒルベルト変換が周波数領域でどう振る舞うかを押さえておく必要があります。ヒルベルト変換は時間領域では次の主値積分で定義されます。
$$ \mathcal{H}[x](t) = \frac{1}{\pi}\, \mathrm{p.v.}\!\int_{-\infty}^{\infty} \frac{x(\tau)}{t – \tau}\, d\tau $$
これは $x(t)$ とカーネル $1/(\pi t)$ の畳み込みです。畳み込みはフーリエ変換で掛け算になるので、$1/(\pi t)$ のフーリエ変換さえわかれば、周波数領域での作用が決まります。よく知られているように、$1/(\pi t)$ のフーリエ変換は $-j\,\mathrm{sgn}(f)$ です。すなわちヒルベルト変換のフーリエ領域での伝達関数は
$$ \widehat{\mathcal{H}}(f) = -j\,\mathrm{sgn}(f) = \begin{cases} -j & (f > 0) \\ 0 & (f = 0) \\ +j & (f < 0) \end{cases} $$
です。ここで $\mathrm{sgn}$ は符号関数です。$-j = e^{-j\pi/2}$ は「位相を $90^\circ$ 遅らせる」操作、$+j = e^{+j\pi/2}$ は「位相を $90^\circ$ 進める」操作に対応します。振幅は変えず（$|-j| = 1$）、位相だけを正の周波数では $-90^\circ$、負の周波数では $+90^\circ$ 回します。
図は伝達関数 $-j\,\mathrm{sgn}(f)$ の虚部を周波数に対して描いたものです。正の周波数では $-1$（$-j$ に相当し位相 $-90^\circ$）、負の周波数では $+1$（$+j$ に相当し位相 $+90^\circ$）と、原点をはさんで階段状に符号が反転します。$f=0$ の直流成分だけは $0$ で、ヒルベルト変換が平均値を通さないことがここに現れています。
この作用を余弦波で確認しましょう。$\cos(2\pi f_0 t)$ のスペクトルは $\pm f_0$ に立つ2本のインパルスです。正の周波数側 $+f_0$ の成分には $-j$ が掛かり、負の周波数側 $-f_0$ の成分には $+j$ が掛かります。時間領域に戻すと、この操作はちょうど余弦を正弦に変えます。
$$ \mathcal{H}[\cos(2\pi f_0 t)] = \cos\!\left(2\pi f_0 t – \frac{\pi}{2}\right) = \sin(2\pi f_0 t) $$
同様に $\mathcal{H}[\sin(2\pi f_0 t)] = -\cos(2\pi f_0 t)$ です。位相を $90^\circ$ 遅らせるという直感が、そのまま余弦→正弦の変換になっているわけです。
時間領域で見ると、変換後の波形（橙破線）は元の余弦（シアン）を横にちょうど $1/4$ 周期ぶんずらした正弦になっています。振幅は $1$ のまま変わらず、位相だけが $90^\circ$ 遅れているのが読み取れます。これがフーリエ領域で $-j$ を掛けた効果の時間波形での姿です。
なお $f = 0$（直流成分）では伝達関数が $0$ なので、ヒルベルト変換は直流を通しません。これは後で瞬時振幅を扱うとき、信号の平均値（オフセット）が結果を乱す要因になることと関係します。ヒルベルト変換の周波数特性が $-j\,\mathrm{sgn}(f)$ であることを押さえたので、次はこれを解析信号のスペクトルに直接使います。
解析信号のスペクトルは正周波数のみ
解析信号を作る本当のご利益は、スペクトルから負の周波数がきれいに消えることにあります。これを導出しましょう。目標は、$z(t) = x(t) + j\,\mathcal{H}[x](t)$ のフーリエ変換 $Z(f)$ が正の周波数だけを持つことを示すことです。
$x(t)$ のフーリエ変換を $X(f)$ とします。ヒルベルト変換は周波数領域で $-j\,\mathrm{sgn}(f)$ を掛ける操作なので、$\mathcal{H}[x](t)$ のフーリエ変換は $-j\,\mathrm{sgn}(f)\,X(f)$ です。したがって $z(t) = x(t) + j\,\mathcal{H}[x](t)$ のフーリエ変換は、フーリエ変換の線形性から
$$ Z(f) = X(f) + j\cdot\big(-j\,\mathrm{sgn}(f)\,X(f)\big) $$
となります。ここで $j \cdot (-j) = 1$ を使うと、括弧の中の係数が整理できます。
$$ Z(f) = X(f) + \mathrm{sgn}(f)\,X(f) = \big(1 + \mathrm{sgn}(f)\big)X(f) $$
あとは $1 + \mathrm{sgn}(f)$ を周波数の符号ごとに場合分けするだけです。$\mathrm{sgn}(f)$ は $f>0$ で $+1$、$f=0$ で $0$、$f<0$ で $-1$ なので、
$$ Z(f) = \begin{cases} 2X(f) & (f > 0) \\ X(0) & (f = 0) \\ 0 & (f < 0) \end{cases} $$
が得られます。これが解析信号の核心です。負の周波数成分が完全に消え、正の周波数成分は振幅が2倍になります。 直流成分だけはそのまま残ります。
左は実信号 $\cos$ のスペクトルで、$\pm f_0$ に振幅 $1/2$ の線が鏡像対称に2本立っています。右の解析信号では負の周波数側の線が消え、正の周波数側だけが振幅 $1$（2倍）の1本に折りたたまれています。2本を1本にまとめてもエネルギーが保たれるよう、残った成分の振幅が2倍になるのがポイントです。
なぜ「正周波数だけ」がうれしいのでしょうか。実信号のスペクトルは必ず共役対称、すなわち $X(-f) = X^*(f)$ です。負の周波数は正の周波数の”鏡像”にすぎず、独立な情報を持ちません。実数の余弦波 $\cos(2\pi f_0 t)$ が $\pm f_0$ の2本で表されるのも、この鏡像のせいです。負の周波数を捨てて正の周波数を2倍にすることは、鏡像を折りたたんで1本のフェーザにまとめる操作にほかなりません。こうして $\cos(2\pi f_0 t)$ は $2\times$（$+f_0$ の1本）となり、時間領域では $e^{j 2\pi f_0 t}$ という素直な回転フェーザに戻ります。
因子2が入る理由も、この視点で腑に落ちます。エネルギーを保つために、$\pm f_0$ の2本（各振幅 $1/2$）を1本に折りたたむと振幅が $1$ になる、というつじつま合わせです。実際 $\cos(2\pi f_0 t) = \tfrac12 e^{j2\pi f_0 t} + \tfrac12 e^{-j2\pi f_0 t}$ で、負側を消して正側を2倍にすると $e^{j2\pi f_0 t}$ そのものになります。
スペクトルが片側だけになったことで、$z(t)$ は「一方向に回るフェーザ」になりました。この回転フェーザの大きさと偏角を読み取れば、瞬間ごとの振幅と位相が得られます。次の節でそれらを正式に定義します。
瞬時振幅・瞬時位相・瞬時周波数の定義
解析信号 $z(t)$ は複素数なので、極形式で書けます。大きさと偏角に分ければ、
$$ z(t) = a(t)\, e^{j\phi(t)} $$
と表せます。ここで登場する2つの量が、瞬時振幅と瞬時位相です。
複素平面で見ると、解析信号は原点から伸びる矢印（フェーザ）が時間とともに回転する様子として描けます。矢印の長さ $a(t)$ がその瞬間の振幅、矢印の向き $\phi(t)$ が瞬時位相です。単純な余弦の解析信号 $e^{j2\pi f_0 t}$ では長さが $1$ のまま一定速度で反時計回りに回り、その回転の速さがそのまま周波数に対応します。
瞬時振幅（instantaneous amplitude） は解析信号の絶対値です。
$$ a(t) = |z(t)| = \sqrt{x(t)^2 + \mathcal{H}[x](t)^2} $$
これは信号の”包絡線（エンベロープ）”を表します。$x = \cos$、$\mathcal{H}[x] = \sin$ なら $a(t) = \sqrt{\cos^2 + \sin^2} = 1$ となり、余弦波の振幅がそのまま1本の水平線として得られます。実部だけを見ていたら $-1$ から $+1$ まで振動していた値が、虚部を補うことで一定値の”強さ”に変換されるのがポイントです。AMラジオの復調が包絡線検波であることを思い出すと、$a(t)$ は変調された振幅そのものを取り出す量だとわかります。
瞬時位相（instantaneous phase） は解析信号の偏角です。
$$ \phi(t) = \arg z(t) = \arctan\!\frac{\mathcal{H}[x](t)}{x(t)} $$
回転フェーザがいまどの角度を向いているかを表します。$e^{j2\pi f_0 t}$ なら $\phi(t) = 2\pi f_0 t$ と時間に比例して増え続けます。
そして本記事の主役、瞬時周波数（instantaneous frequency） は、瞬時位相の時間微分（を $2\pi$ で割ったもの）です。
$$ f(t) = \frac{1}{2\pi}\frac{d\phi(t)}{dt} $$
なぜ位相の傾きが周波数なのでしょうか。定常な余弦波 $\cos(2\pi f_0 t)$ を考えると、位相は $\phi(t) = 2\pi f_0 t$ で、その傾きは $d\phi/dt = 2\pi f_0$ です。これを $2\pi$ で割れば $f_0$ に戻ります。つまり普通の周波数の定義と完全に一致します。瞬時周波数はこれを一般化し、「位相が時々刻々どれだけ速く進んでいるか」をそのまま周波数と呼ぶ、という自然な拡張です。角周波数で言えば $\omega(t) = d\phi/dt$ が瞬時角周波数です。
この定義の威力は、位相が時間に比例しない場合、つまり周波数が変化する信号で発揮されます。たとえば周波数を線形に掃引するチャープ信号 $x(t) = \cos\!\big(2\pi(f_0 t + \tfrac{k}{2}t^2)\big)$ では、瞬時位相が $\phi(t) = 2\pi(f_0 t + \tfrac{k}{2}t^2)$ なので、微分すると
$$ f(t) = \frac{1}{2\pi}\frac{d}{dt}\Big[2\pi\big(f_0 t + \tfrac{k}{2}t^2\big)\Big] = f_0 + k t $$
と、時間とともに直線的に増える周波数が得られます。$k$ が掃引率（1秒あたりの周波数増加）です。フーリエ変換ならこの信号は $f_0$ から $f_0 + kT$ まで広がった1つの帯としてしか見えませんが、瞬時周波数なら「時刻 $t$ には $f_0 + kt$ で鳴っている」と時間ごとに読み取れます。
ここまでの定義は魔法のように見えますが、実は成立する条件があります。どんな信号にでも $f(t)$ を計算してよいわけではないのです。次の節で、瞬時周波数が意味を持つための条件（単一成分であること、Bedrosianの制約）を見ていきます。
瞬時周波数が意味を持つ条件（Bedrosianの制約）
瞬時周波数は強力ですが、乱用すると無意味な数値を吐き出します。落とし穴を2つ押さえておきましょう。
単一成分でないと破綻する
最大の注意点は、瞬時周波数は「その瞬間にひとつの周波数しか鳴っていない」単一成分（mono-component）信号にしか意味を持たない ことです。試しに2つの余弦波の和 $x(t) = \cos(2\pi f_1 t) + \cos(2\pi f_2 t)$ を考えます。この信号の解析信号は $z(t) = e^{j2\pi f_1 t} + e^{j2\pi f_2 t}$ で、これは2本のフェーザの合成です。2本のフェーザを足すと、合成ベクトルの長さは時間とともに脈打ち（うなり＝ビート）、向きも一定の速さでは回りません。
具体的に合成フェーザの偏角を微分すると、瞬時周波数は $f_1$ と $f_2$ の間を激しく揺れ動き、しかも合成ベクトルが原点近くを通る瞬間には 負の周波数 や発散的なスパイクさえ生じます。物理的には $f_1$ と $f_2$ の2つの音が同時に鳴っているのに、$f(t)$ は「1つの周波数」を無理やり答えようとして破綻するのです。「いまの周波数」という概念は、そもそも単一成分でなければ定義できません。多成分信号を扱いたいなら、まず帯域通過フィルタやEMDで単一成分に分解してから各成分に瞬時周波数を適用する、という手順が必要です。
Bedrosianの定理
もうひとつ、振幅変調された信号 $x(t) = a(t)\cos(\phi(t))$ の包絡線と位相を正しく取り出せる条件が Bedrosianの定理 です。私たちは「振幅 $a(t)$ が瞬時振幅、位相 $\phi(t)$ が瞬時位相として素直に取り出せる」と暗に期待していますが、それが成り立つのはヒルベルト変換が
$$ \mathcal{H}[a(t)\cos\phi(t)] = a(t)\,\mathcal{H}[\cos\phi(t)] = a(t)\sin\phi(t) $$
のように、ゆっくり変わる振幅 $a(t)$ を”素通し”できるときだけです。Bedrosianの定理は、この分離が成り立つための十分条件を周波数領域で与えます。振幅 $a(t)$ のスペクトルが低周波側に、搬送波 $\cos\phi(t)$ のスペクトルが高周波側にあって、両者の帯域が重ならない（$a(t)$ の最高周波数 < $\cos\phi(t)$ の最低周波数）とき、上の分離が厳密に成り立ちます。
図はBedrosian条件が成り立つ状況の模式的なスペクトルです。包絡線 $a(t)$ の帯域（シアン）は原点近くの低周波に、搬送波 $\cos\phi(t)$ の帯域（橙）はずっと高い周波数に位置し、両者の間に周波数の隙間があります。この「帯域が重ならない」条件が満たされると、ヒルベルト変換は包絡線に触れず搬送波だけに $90^\circ$ の位相シフトを与え、$a(t)$ と位相をきれいに分離できます。
直感的には、包絡線が搬送波よりずっとゆっくり変化していれば、ヒルベルト変換は速く振動する搬送波だけに $90^\circ$ の位相シフトを与え、ゆっくりした包絡線には触れない、ということです。AM信号なら「変調信号の帯域が搬送波周波数よりずっと低い」という通常の設計条件がそのままBedrosian条件になっています。逆に包絡線が搬送波と同じくらい速く動くと帯域が重なり、瞬時振幅・瞬時位相の分離が崩れて、包絡線に搬送波のさざなみが漏れ込みます。
この2つの制約は、後のPython実験で実際に破綻の様子を見ると腑に落ちます。まずは制約が満たされる素直なチャープ信号で、瞬時周波数がきちんと復元できることを確かめましょう。
Pythonでの実装
ここからは実装です。scipy.signal.hilbert は名前に反して「ヒルベルト変換」そのものではなく、解析信号 $z(t) = x(t) + j\,\mathcal{H}[x](t)$ を返す 関数である点に注意してください（実部が元信号、虚部がヒルベルト変換）。まずは共通の準備として、日本語フォントとヘルパ関数を用意します。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.signal import hilbert
# 日本語フォント設定
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
def instantaneous_frequency(x, fs):
"""実信号 x から解析信号を作り、瞬時振幅・瞬時位相・瞬時周波数を返す"""
z = hilbert(x) # 解析信号 z(t) = x + j H[x]
amp = np.abs(z) # 瞬時振幅（包絡線）
phase = np.unwrap(np.angle(z)) # 瞬時位相（アンラップ必須）
inst_freq = np.diff(phase) / (2*np.pi) * fs # 瞬時周波数 = (1/2pi) dphi/dt
return z, amp, phase, inst_freq
np.angle は偏角を $(-\pi, \pi]$ に折り返して返すため、位相が $\pi$ を超えるたびに $2\pi$ の飛びが入ります。この飛びをそのまま微分すると巨大なスパイクになるので、np.unwrap で連続な位相に直してから微分します。この「位相アンラップ」の重要性は後で図で確認します。瞬時周波数は位相の差分を $2\pi$ で割り、サンプリング周波数 $f_s$ を掛けて Hz に直しています。
単一成分の余弦波で解析信号を確認する
まず、もっとも単純な余弦波で解析信号の振る舞いを確かめます。狙いは、実部が元の余弦・虚部が正弦（$90^\circ$ 遅れ）になり、包絡線が一定になることを目で見ることです。
fs = 1000.0 # サンプリング周波数 [Hz]
t = np.arange(0, 1.0, 1/fs) # 1秒間
f0 = 10.0 # 10 Hz の余弦波
x = np.cos(2*np.pi*f0*t)
z, amp, phase, inst_freq = instantaneous_frequency(x, fs)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(t, z.real, label="実部 x(t)=cos")
ax[0].plot(t, z.imag, "--", label="虚部 H[x](t)=sin")
ax[0].plot(t, amp, "r", lw=2, label="瞬時振幅 |z(t)|")
ax[0].set_xlim(0, 0.3); ax[0].set_xlabel("時間 [s]"); ax[0].set_ylabel("振幅")
ax[0].set_title("解析信号の実部・虚部・包絡線"); ax[0].legend(loc="upper right")
ax[1].plot(t[:-1], inst_freq, "g", lw=2)
ax[1].axhline(f0, color="k", ls=":", label="設計値 10 Hz")
ax[1].set_ylim(0, 20); ax[1].set_xlabel("時間 [s]"); ax[1].set_ylabel("瞬時周波数 [Hz]")
ax[1].set_title("瞬時周波数"); ax[1].legend()
plt.tight_layout(); plt.show()
左の図から、虚部（破線）が実部（実線）よりちょうど $90^\circ$ 遅れた正弦波になっていること、そして赤い瞬時振幅が振動せず $1$ の水平線になっていることが読み取れます。実部だけを見ていたら $\pm1$ に振動していた”強さ”が、解析信号にすると一定値として取り出せるわけです。右の図では、両端を除いて瞬時周波数が $10$ Hz にぴたりと張り付いており、定常な余弦波の瞬時周波数がその周波数そのものになるという理論どおりの結果です。両端がわずかに乱れるのは、ヒルベルト変換がFFTベースで信号を周期的とみなすための端点効果（エッジ効果）で、実務では両端を数サンプル捨てて解釈します。
線形チャープ信号で掃引周波数を復元する
次が本命です。周波数を時間とともに直線的に掃引するチャープ信号を作り、瞬時周波数が設計どおりの直線を復元するかを確かめます。$f(t) = f_0 + kt$ となるよう位相を $\phi(t) = 2\pi(f_0 t + \tfrac{k}{2}t^2)$ に設定します。
fs = 2000.0
T = 2.0
t = np.arange(0, T, 1/fs)
f_start, f_end = 20.0, 200.0 # 20 Hz -> 200 Hz へ掃引
k = (f_end - f_start) / T # 掃引率 [Hz/s]
phase_true = 2*np.pi*(f_start*t + 0.5*k*t**2)
x = np.cos(phase_true)
z, amp, phase, inst_freq = instantaneous_frequency(x, fs)
f_design = f_start + k*t # 設計上の瞬時周波数（直線）
plt.figure(figsize=(11, 4))
plt.subplot(1, 2, 1)
plt.plot(t, x, lw=0.6); plt.plot(t, amp, "r", lw=2, label="瞬時振幅")
plt.xlim(0, 0.3); plt.xlabel("時間 [s]"); plt.ylabel("振幅")
plt.title("チャープ信号と包絡線"); plt.legend()
plt.subplot(1, 2, 2)
plt.plot(t[:-1], inst_freq, "g", lw=1.5, label="推定した瞬時周波数")
plt.plot(t, f_design, "k--", lw=1.5, label="設計した掃引 f0+kt")
plt.xlabel("時間 [s]"); plt.ylabel("周波数 [Hz]")
plt.title("瞬時周波数 vs 設計掃引"); plt.legend()
plt.tight_layout(); plt.show()
err = np.abs(inst_freq - f_design[:-1])
print(f"瞬時周波数の平均誤差（両端5%除外）: "
f"{np.mean(err[len(err)//20 : -len(err)//20]):.3f} Hz")
右の図で、推定した瞬時周波数（緑）が設計した直線 $f_0 + kt$（黒破線）とほぼ完全に重なります。$20$ Hz から $200$ Hz へと $2$ 秒かけて滑らかに上昇する様子が、時間ごとの周波数として正しく取り出せています。フーリエ変換ではこの信号は $20$〜$200$ Hz の帯としか見えませんが、解析信号を使えば「いつ何Hzか」が復元できるわけです。出力される平均誤差も端点を除けば1Hz以下に収まり、掃引の中央部では推定が極めて正確であることがわかります。左の図の包絡線が(端点を除き)ほぼ一定なのは、チャープが振幅変調を含まない純粋な周波数変調だからです。
AM+FM信号：包絡線と瞬時周波数を同時に取り出す
現実の信号は、振幅も周波数も同時に変化することがよくあります。ここでは振幅変調（AM）と周波数変調（FM）を同時にかけた信号で、包絡線 $a(t)$ と瞬時周波数 $f(t)$ を別々に取り出せることを確かめます。Bedrosian条件を満たすよう、包絡線はゆっくり（$2$ Hz）、搬送波は速く（$60$ Hz 中心）に設定します。
fs = 2000.0
t = np.arange(0, 2.0, 1/fs)
# ゆっくりした振幅変調（2 Hz）: Bedrosian条件を満たすため搬送波よりずっと低周波
a_true = 1.0 + 0.5*np.cos(2*np.pi*2.0*t)
# 正弦波状に周波数変調された搬送波（中心 60 Hz、偏移 ±15 Hz、変調 3 Hz）
fc, f_dev, fm = 60.0, 15.0, 3.0
phase_true = 2*np.pi*(fc*t - (f_dev/(2*np.pi*fm))*np.cos(2*np.pi*fm*t))
f_true = fc + f_dev*np.sin(2*np.pi*fm*t) # 設計上の瞬時周波数
x = a_true * np.cos(phase_true)
z, amp, phase, inst_freq = instantaneous_frequency(x, fs)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(t, x, lw=0.5, color="0.6", label="AM+FM 信号")
ax[0].plot(t, amp, "r", lw=2, label="推定包絡線 |z|")
ax[0].plot(t, a_true, "b--", lw=1.5, label="設計包絡線 a(t)")
ax[0].set_xlabel("時間 [s]"); ax[0].set_ylabel("振幅")
ax[0].set_title("瞬時振幅（包絡線）の復元"); ax[0].legend(loc="upper right")
ax[1].plot(t[:-1], inst_freq, "g", lw=1.0, label="推定した瞬時周波数")
ax[1].plot(t, f_true, "k--", lw=1.5, label="設計した瞬時周波数")
ax[1].set_ylim(30, 90); ax[1].set_xlabel("時間 [s]"); ax[1].set_ylabel("周波数 [Hz]")
ax[1].set_title("瞬時周波数の復元"); ax[1].legend(loc="upper right")
plt.tight_layout(); plt.show()
左の図では、赤い推定包絡線が青破線の設計包絡線 $a(t) = 1 + 0.5\cos(2\pi\cdot2 t)$ にぴたりと重なり、$2$ Hz でゆっくり脈打つ振幅がきれいに取り出せています。右の図でも、推定した瞬時周波数が $60 \pm 15$ Hz を $3$ Hz で正弦振動する設計値とよく一致します。振幅の情報（包絡線）と周波数の情報（瞬時周波数）が、ひとつの解析信号から分離して取り出せている のが最大のポイントです。これはBedrosian条件（包絡線が搬送波よりずっと低周波）が満たされているおかげで、ヒルベルト変換がゆっくりした振幅を素通しし、速い搬送波だけに位相シフトを与えているからです。
位相アンラップがなぜ必要か
先ほどのコードで np.unwrap を使いました。これを外すとどうなるかを見て、アンラップの必要性を実感しましょう。np.angle の生の出力は $(-\pi, \pi]$ に折り返されるため、位相が増え続ける信号では $\pi$ を超えるたびに $-\pi$ へジャンプします。
fs = 1000.0
t = np.arange(0, 0.5, 1/fs)
x = np.cos(2*np.pi*30*t) # 30 Hz
z = hilbert(x)
phase_wrapped = np.angle(z) # ラップされたまま（-pi, pi]
phase_unwrapped = np.unwrap(phase_wrapped)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(t, phase_wrapped, "r")
ax[0].set_xlabel("時間 [s]"); ax[0].set_ylabel("位相 [rad]")
ax[0].set_title("アンラップ前：2πの飛びが残る")
ax[1].plot(t, phase_unwrapped, "b")
ax[1].set_xlabel("時間 [s]"); ax[1].set_ylabel("位相 [rad]")
ax[1].set_title("アンラップ後：直線的に増加")
plt.tight_layout(); plt.show()
# 飛びを含んだまま微分すると瞬時周波数が破綻する
if_wrapped = np.diff(phase_wrapped)/(2*np.pi)*fs
if_unwrapped = np.diff(phase_unwrapped)/(2*np.pi)*fs
print(f"アンラップ前の瞬時周波数の最大値: {np.max(np.abs(if_wrapped)):.0f} Hz（本来30 Hz）")
print(f"アンラップ後の瞬時周波数の最大値: {np.max(np.abs(if_unwrapped)):.1f} Hz")
左の図はのこぎり歯状で、位相が $\pi$ に達するたびに $-\pi$ へ急落しています。この不連続を微分すると、飛びの箇所で巨大な負のスパイク（数百Hz級）が現れ、瞬時周波数がまったく使い物になりません。右のアンラップ後の位相はきれいな直線で、その傾きを取れば一定の $30$ Hz が得られます。出力の数値でも、アンラップ前は本来 $30$ Hz のはずが飛びのせいで極端な値になり、アンラップ後は $30$ Hz 付近に収まることが確認できます。瞬時周波数は位相の微分で定義される以上、微分する前に位相を連続にしておくアンラップが不可欠 なのです。
多成分信号では瞬時周波数が破綻する
最後に、単一成分でないと瞬時周波数が意味を失うことを実演します。$40$ Hz と $80$ Hz の2つの余弦波を足した信号に、そのまま瞬時周波数を計算してみます。
fs = 2000.0
t = np.arange(0, 0.5, 1/fs)
x = np.cos(2*np.pi*40*t) + np.cos(2*np.pi*80*t) # 2成分の和
z, amp, phase, inst_freq = instantaneous_frequency(x, fs)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(t, amp, "r", lw=1.5)
ax[0].set_xlabel("時間 [s]"); ax[0].set_ylabel("|z(t)|")
ax[0].set_title("包絡線がうなる（ビート）")
ax[1].plot(t[:-1], inst_freq, "g", lw=1.0)
ax[1].axhline(40, color="k", ls=":", label="40 Hz")
ax[1].axhline(80, color="b", ls=":", label="80 Hz")
ax[1].set_ylim(-100, 200); ax[1].set_xlabel("時間 [s]"); ax[1].set_ylabel("瞬時周波数 [Hz]")
ax[1].set_title("瞬時周波数が激しく振れ、負にもなる"); ax[1].legend()
plt.tight_layout(); plt.show()
左の包絡線は一定にならず、2成分の差 $80-40=40$ Hz でうなり（ビート）を打っています。右の瞬時周波数はもっと深刻で、$40$ Hz と $80$ Hz の平均あたりを中心に激しく振動し、包絡線が $0$ に近づく瞬間には 負の周波数 にまで振れています。これは「いまの周波数」を無理やり1つの数で答えようとした結果の破綻です。2つの音が同時に鳴っている信号に「瞬間の周波数」を問うこと自体が無意味なので、瞬時周波数を使う前にはまず帯域通過フィルタやEMDで信号を単一成分に分解しておく必要があります。この図は、瞬時周波数を安易に適用してはいけないという教訓を端的に示しています。
まとめ
本記事では、解析信号と瞬時周波数について、定義・導出・実装を通して解説しました。
- 解析信号 $z(t) = x(t) + j\,\mathcal{H}[x](t)$ は、実信号にヒルベルト変換（$90^\circ$ 位相シフト）を虚部として足した複素信号で、実信号を回転フェーザに”持ち上げる”
- ヒルベルト変換の周波数特性 $-j\,\mathrm{sgn}(f)$ を使うと、解析信号のスペクトルは $Z(f) = (1+\mathrm{sgn}(f))X(f)$ となり、負の周波数が消え、正の周波数が2倍 になる
- 極形式 $z(t) = a(t)e^{j\phi(t)}$ から、瞬時振幅 $a(t)=|z(t)|$（包絡線）、瞬時位相 $\phi(t)=\arg z(t)$、瞬時周波数 $f(t)=\frac{1}{2\pi}\frac{d\phi}{dt}$ が定義される
- 線形チャープでは瞬時周波数が設計どおりの直線 $f_0+kt$ を、AM+FM信号では包絡線と瞬時周波数を分離して復元できる
- 瞬時周波数は 単一成分 信号にしか意味を持たず、多成分だと激しく振れて負にもなる。振幅・位相の正しい分離には Bedrosianの制約（包絡線が搬送波よりずっと低周波）が必要
- 位相の微分で周波数を得るため、微分の前に 位相アンラップ が不可欠
解析信号と瞬時周波数は、時間周波数解析への入り口です。ここから先は、単一成分への分解（EMD／ヒルベルト・ホァン変換）、時間周波数分布（短時間フーリエ変換やウェーブレット）へと広がっていきます。次のステップとして、以下の記事も参考にしてください。