AM ラジオ放送の電波をオシロスコープで眺めると、搬送波の振幅がゆっくり変動している様子が見えます。この「振幅の変動」——つまり包絡線——を取り出せれば、もとの音声信号が復元できるはずです。しかし、数学的に「信号の振幅だけを取り出す」とはどういう操作でしょうか？ 信号 $x(t)$ をそのまま見ても、正と負に振動する値が得られるだけで、ゆるやかな振幅の輪郭は直接は読み取れません。
ここで登場するのがヒルベルト変換（Hilbert Transform）です。ヒルベルト変換を使えば、実数信号 $x(t)$ に対して「位相が 90° ずれた相棒」$\hat{x}(t)$ を作ることができ、この 2 つを組み合わせて解析信号（analytic signal）という複素信号を構成できます。解析信号の絶対値が包絡線、偏角の微分が瞬時周波数——こうして信号の振幅情報と周波数情報を綺麗に分離できるのです。
ヒルベルト変換を理解すると、次のような幅広い分野が見通しよくなります。
- 通信工学 — SSB（単側波帯）変調の生成、包絡線検波、瞬時周波数の推定はすべてヒルベルト変換に帰着します
- 振動解析・機械診断 — 回転機械の振動信号から包絡線を抽出し、ベアリング故障の兆候を検出するのにヒルベルト変換が使われています
- 地震工学・音響信号処理 — 地震波や音響信号の瞬時振幅・瞬時周波数を追跡し、時間-周波数解析の基礎となります
- レーダー・ソナー — 受信信号の I/Q 成分生成（直交検波）はヒルベルト変換と密接に結びついています
本記事の内容
- ヒルベルト変換の直感的理解 — 「信号の位相を 90° 回す変換」
- コーシーの主値積分による数学的定義
- 周波数領域での解釈 — $-j\,\mathrm{sgn}(f)$ との乗算
- 解析信号 $z(t) = x(t) + j\hat{x}(t)$ の構成と性質
- 包絡線の抽出（AM 信号への応用）
- 瞬時周波数の抽出（FM 信号への応用）
- SSB 変調への応用 — 位相法による単側波帯生成
- 離散ヒルベルト変換の実装方法
- Python による包絡線検出・瞬時周波数計算・SSB 信号生成
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ヒルベルト変換の直感 — 信号の位相を 90° 回す
ヒルベルト変換を一言で要約すると、信号に含まれるすべての周波数成分の位相を一律に 90°（$\pi/2$）ずらす操作です。
日常のアナロジーで考えてみましょう。ブランコを漕いでいるとき、ブランコの位置（変位）は $\cos(\omega t)$ のように変化します。一方、ブランコの速度は位置を微分した $-\omega\sin(\omega t)$ であり、位置が最も高い点では速度がゼロ、位置がゼロを通過する瞬間に速度が最大になります。つまり、速度は位置に対して位相が 90° 進んでいます。ヒルベルト変換は、この「位相を 90° ずらす」操作を、特定の周波数だけでなく信号に含まれるあらゆる周波数成分に対して同時に行うのです。しかも微分とは違い、振幅は変えません。$\cos(\omega t)$ をヒルベルト変換すると $\sin(\omega t)$ になり、振幅は 1 のまま保たれます。
なぜ「90° ずらした信号」が欲しいのでしょうか。それは、元の信号 $x(t) = A(t)\cos\theta(t)$ と、90° ずれた信号 $\hat{x}(t) = A(t)\sin\theta(t)$ の 2 つがあれば、
$$ A(t) = \sqrt{x(t)^2 + \hat{x}(t)^2} $$
と計算するだけで包絡線（振幅 $A(t)$）が取り出せるからです。$\cos^2\theta + \sin^2\theta = 1$ というピタゴラスの定理が効いているわけです。同様に、
$$ \theta(t) = \arctan\frac{\hat{x}(t)}{x(t)} $$
から瞬時位相を取り出し、その時間微分から瞬時周波数も得られます。
直感を掴んだところで、次にヒルベルト変換の厳密な数学的定義を見ていきましょう。
ヒルベルト変換の数学的定義 — コーシーの主値積分
ヒルベルト変換は、実数信号 $x(t)$ を入力として別の実数信号 $\hat{x}(t)$ を出力する線形変換です。「すべての周波数成分の位相を 90° ずらす」という操作を、時間領域の畳み込み積分として表現すると、次の定義式が得られます。
$$ \hat{x}(t) = \mathcal{H}\{x(t)\} = \frac{1}{\pi} \, \mathrm{p.v.}\!\int_{-\infty}^{\infty} \frac{x(\tau)}{t – \tau} \, d\tau $$
ここで $\mathrm{p.v.}$ はコーシーの主値（Cauchy principal value）を表します。各記号の意味を確認しましょう。
- $x(\tau)$: 入力信号を積分変数 $\tau$ で表したもの
- $t – \tau$: 現在時刻 $t$ と積分変数 $\tau$ の差。$\tau = t$ で分母がゼロになる
- $\mathrm{p.v.}$: $\tau = t$ の特異点を避けて積分する操作
なぜコーシーの主値が必要なのか
被積分関数の分母 $t – \tau$ は $\tau = t$ でゼロになるため、通常の意味では積分が発散してしまいます。コーシーの主値は、特異点の両側から対称的に近づいて積分を定義する方法です。
$$ \mathrm{p.v.}\!\int_{-\infty}^{\infty} \frac{x(\tau)}{t – \tau} \, d\tau = \lim_{\varepsilon \to 0^+} \left[ \int_{-\infty}^{t-\varepsilon} \frac{x(\tau)}{t – \tau} \, d\tau + \int_{t+\varepsilon}^{\infty} \frac{x(\tau)}{t – \tau} \, d\tau \right] $$
特異点 $\tau = t$ の左側と右側で $1/(t – \tau)$ は逆符号になるため、$x(\tau)$ が $\tau = t$ の近傍で十分滑らかであれば、特異点近傍の寄与は互いに打ち消し合い、主値積分は有限の値に収束します。
畳み込みとしての表現
ヒルベルト変換は、入力信号 $x(t)$ とインパルス応答 $h(t) = 1/(\pi t)$ の畳み込みとして書くことができます。
$$ \hat{x}(t) = x(t) * h(t) = x(t) * \frac{1}{\pi t} $$
この $h(t) = 1/(\pi t)$ がヒルベルト変換フィルタのインパルス応答です。$h(t)$ は $t = 0$ に特異点を持つ奇関数であり、時間的に減衰しないため非因果的（non-causal）な理想フィルタであることがわかります。言い換えれば、ヒルベルト変換は物理的にリアルタイム実現できない理想的な操作ですが、信号処理の理論ツールとしては非常に強力です。
畳み込みの表現が得られたということは、フーリエ変換を使えば周波数領域での解釈が可能になるはずです。次に、この変換が周波数領域でどのように見えるかを調べてみましょう。
周波数領域での解釈 — $-j\,\mathrm{sgn}(f)$ との乗算
畳み込み定理により、時間領域での畳み込みは周波数領域での乗算に対応します。ヒルベルト変換フィルタ $h(t) = 1/(\pi t)$ のフーリエ変換 $H(f)$ を求めましょう。
$1/(\pi t)$ のフーリエ変換の導出
$h(t) = 1/(\pi t)$ のフーリエ変換を直接計算するために、まず符号関数（signum function）のフーリエ変換対を利用します。
符号関数は次のように定義されます。
$$ \mathrm{sgn}(t) = \begin{cases} +1 & (t > 0) \\ 0 & (t = 0) \\ -1 & (t < 0) \end{cases} $$
符号関数のフーリエ変換は、$\mathrm{sgn}(t) = 2u(t) – 1$（$u(t)$ は単位ステップ関数）の関係と超関数論の結果を組み合わせると、次のように導かれます。
$$ \mathcal{F}\{\mathrm{sgn}(t)\} = \frac{1}{j\pi f} $$
これをフーリエ変換の双対性（duality）で読み替えます。双対性とは「$g(t)$ のフーリエ変換が $G(f)$ なら、$G(t)$ のフーリエ変換は $g(-f)$」という性質です。上の結果に双対性を適用すると、
$$ \mathcal{F}\left\{\frac{1}{j\pi t}\right\} = \mathrm{sgn}(-f) = -\mathrm{sgn}(f) $$
$\mathrm{sgn}(f)$ は奇関数なので $\mathrm{sgn}(-f) = -\mathrm{sgn}(f)$ を用いました。両辺に $-j$ を掛けて $j$ を分母から外すと、
$$ \mathcal{F}\left\{\frac{1}{\pi t}\right\} = -j\,\mathrm{sgn}(f) $$
したがって、ヒルベルト変換フィルタの伝達関数は次の通りです。
$$ \boxed{H(f) = -j\,\mathrm{sgn}(f)} $$
伝達関数の意味
この伝達関数が何を表しているか、具体的に整理しましょう。
$$ H(f) = -j\,\mathrm{sgn}(f) = \begin{cases} -j & (f > 0) \\ 0 & (f = 0) \\ +j & (f < 0) \end{cases} $$
$-j = e^{-j\pi/2}$ なので、正の周波数成分には位相 $-90°$（$-\pi/2$）の回転を施し、$+j = e^{+j\pi/2}$ なので、負の周波数成分には位相 $+90°$（$+\pi/2$）の回転を施します。一方、$|H(f)| = 1$（$f \neq 0$）なので振幅は一切変えません。
まさに「全帯域位相シフタ」——すべての周波数成分の位相を $90°$ ずらすフィルタです。正と負の周波数で回転方向が逆なのは、実数信号の場合にヒルベルト変換後も実数信号であり続けるための整合条件です。実信号のフーリエ変換は $X(-f) = X^*(f)$（エルミート対称性）を満たし、ヒルベルト変換後もこの性質が保たれる必要があります。
周波数領域でのヒルベルト変換
以上をまとめると、ヒルベルト変換は周波数領域で次のように表現されます。
$$ \hat{X}(f) = H(f) \cdot X(f) = -j\,\mathrm{sgn}(f) \cdot X(f) $$
ここで $X(f) = \mathcal{F}\{x(t)\}$、$\hat{X}(f) = \mathcal{F}\{\hat{x}(t)\}$ です。この関係を使えば、入力信号をフーリエ変換し、$-j\,\mathrm{sgn}(f)$ を乗じ、逆フーリエ変換するだけでヒルベルト変換が計算できます。離散信号の場合には FFT を使ってこの操作を効率的に実装でき、これが後で紹介する Python の scipy.signal.hilbert の内部動作原理でもあります。
周波数領域の解釈がわかったところで、次にヒルベルト変換の最も重要な応用——解析信号の構成に進みましょう。
解析信号の構成と性質
解析信号とは
ヒルベルト変換の真価は、実数信号 $x(t)$ から解析信号（analytic signal）$z(t)$ を構成できることにあります。解析信号は、元の信号を実部、ヒルベルト変換を虚部とする複素信号です。
$$ \boxed{z(t) = x(t) + j\hat{x}(t)} $$
「解析」という名前は、この信号が複素解析の意味で特別な性質を持つことに由来します。具体的には、解析信号のフーリエ変換は負の周波数成分がゼロになります。
解析信号の周波数スペクトルの導出
$z(t) = x(t) + j\hat{x}(t)$ のフーリエ変換を計算しましょう。
$$ Z(f) = X(f) + j\hat{X}(f) $$
先ほど導出した $\hat{X}(f) = -j\,\mathrm{sgn}(f) \cdot X(f)$ を代入します。
$$ Z(f) = X(f) + j \cdot \left[-j\,\mathrm{sgn}(f)\right] \cdot X(f) $$
$j \cdot (-j) = -j^2 = -(-1) = 1$ なので、
$$ Z(f) = X(f) + \mathrm{sgn}(f) \cdot X(f) = \left[1 + \mathrm{sgn}(f)\right] X(f) $$
ここで $1 + \mathrm{sgn}(f)$ の値を場合分けします。
$$ 1 + \mathrm{sgn}(f) = \begin{cases} 2 & (f > 0) \\ 1 & (f = 0) \\ 0 & (f < 0) \end{cases} $$
これは単位ステップ関数の 2 倍に相当します。つまり、
$$ \boxed{Z(f) = 2u(f) \cdot X(f)} $$
ここで $u(f)$ はステップ関数（$f > 0$ で 1、$f < 0$ で 0）です。
この結果の意味は明快です。解析信号は元の信号の正の周波数成分だけを 2 倍にして保持し、負の周波数成分を完全に除去した信号です。実信号のフーリエ変換は正と負で共役対称 $X(-f) = X^*(f)$ なので、正の周波数成分だけで信号の全情報を担っています。2 倍にするのは、負の周波数成分を捨てた分のエネルギーを補償するためです。
解析信号の極表示 — 包絡線と瞬時位相
複素信号 $z(t)$ を極形式で書くと、
$$ z(t) = A(t) \, e^{j\theta(t)} $$
ここで、
$$ A(t) = |z(t)| = \sqrt{x(t)^2 + \hat{x}(t)^2} $$
また、$z(t)$ の偏角を取ると瞬時位相が得られます:
$$ \theta(t) = \arg z(t) = \arctan\frac{\hat{x}(t)}{x(t)} $$
$A(t)$ が包絡線（envelope）、$\theta(t)$ が瞬時位相（instantaneous phase）です。そして瞬時位相を時間微分すれば瞬時角周波数が得られます。
$$ \omega_i(t) = \frac{d\theta(t)}{dt} $$
瞬時周波数 $f_i(t)$ は $\omega_i(t)/(2\pi)$ です。
$$ f_i(t) = \frac{1}{2\pi}\frac{d\theta(t)}{dt} $$
これで信号の振幅情報（包絡線）と周波数情報（瞬時周波数）を分離するための数学的道具が揃いました。次に、これらの概念を具体的な信号に適用して理解を深めましょう。
ヒルベルト変換の基本的な計算例
例1: 余弦関数のヒルベルト変換
最も基本的な例として、$x(t) = \cos(2\pi f_0 t)$ のヒルベルト変換を計算しましょう。
周波数領域で考えます。$\cos(2\pi f_0 t)$ のフーリエ変換は、
$$ X(f) = \frac{1}{2}\left[\delta(f – f_0) + \delta(f + f_0)\right] $$
ヒルベルト変換の周波数表現を適用すると、
$$ \hat{X}(f) = -j\,\mathrm{sgn}(f) \cdot X(f) = -j\,\mathrm{sgn}(f) \cdot \frac{1}{2}\left[\delta(f – f_0) + \delta(f + f_0)\right] $$
$f_0 > 0$ のとき、$\delta(f – f_0)$ は $f = f_0 > 0$ にのみ値を持つので $\mathrm{sgn}(f) = 1$、$\delta(f + f_0)$ は $f = -f_0 < 0$ にのみ値を持つので $\mathrm{sgn}(f) = -1$ です。したがって、
$$ \hat{X}(f) = \frac{1}{2}\left[-j \cdot \delta(f – f_0) + j \cdot \delta(f + f_0)\right] $$
$-j = e^{-j\pi/2}$ であることを思い出し、逆フーリエ変換を計算すると、
$$ \hat{x}(t) = \frac{1}{2}\left[-j \cdot e^{j2\pi f_0 t} + j \cdot e^{-j2\pi f_0 t}\right] = \frac{e^{j2\pi f_0 t} – e^{-j2\pi f_0 t}}{2j} = \sin(2\pi f_0 t) $$
途中でオイラーの公式 $\sin\alpha = (e^{j\alpha} – e^{-j\alpha})/(2j)$ を使いました。
$$ \boxed{\mathcal{H}\{\cos(2\pi f_0 t)\} = \sin(2\pi f_0 t)} $$
まさに「位相が 90° 遅れた信号」が得られました。コサインをヒルベルト変換するとサインになる——これは冒頭で述べた直感そのものです。
例2: 正弦関数のヒルベルト変換
同様に、$x(t) = \sin(2\pi f_0 t)$ をヒルベルト変換すると、
$$ \mathcal{H}\{\sin(2\pi f_0 t)\} = -\cos(2\pi f_0 t) $$
となります。サインの 90° 遅れはマイナスコサインです。
例3: 解析信号の構成
$x(t) = \cos(2\pi f_0 t)$ の解析信号を構成すると、
$$ z(t) = \cos(2\pi f_0 t) + j\sin(2\pi f_0 t) = e^{j2\pi f_0 t} $$
見事に正の周波数の複素指数関数だけが残りました。この結果は $Z(f) = 2u(f) \cdot X(f)$ の公式から直接得られる $Z(f) = \delta(f – f_0)$ と整合しています。
基本的な計算例で理論が確認できました。次に、ヒルベルト変換の重要な応用である包絡線の抽出について見ていきましょう。
包絡線の抽出 — AM 信号への応用
AM 変調信号の復習
振幅変調（AM）信号は次の形をしています。
$$ x(t) = A_c \left[1 + m \cdot s(t)\right] \cos(2\pi f_c t) $$
ここで $A_c$ は搬送波振幅、$m$ は変調指数、$s(t)$ は正規化されたベースバンド信号（$|s(t)| \leq 1$）、$f_c$ は搬送波周波数です。
この信号の「包絡線」は $A(t) = A_c[1 + m \cdot s(t)]$ であり、これがベースバンド信号 $s(t)$ の情報を担っています。
ヒルベルト変換による包絡線の導出
AM 信号が狭帯域信号（搬送波周波数 $f_c$ がベースバンド信号の帯域幅よりも十分に高い）であるとき、ヒルベルト変換は搬送波の位相を 90° ずらす操作に帰着します。
$x(t) = A(t)\cos(2\pi f_c t)$ のヒルベルト変換は、$A(t)$ がゆっくり変化する包絡線で $\cos(2\pi f_c t)$ が高周波搬送波であるとき、Bedrosian の定理により近似的に次のようになります。
$$ \hat{x}(t) \approx A(t)\sin(2\pi f_c t) $$
Bedrosian の定理は、低周波成分 $A(t)$ と高周波成分 $\cos(2\pi f_c t)$ のスペクトルが重ならない場合に、積のヒルベルト変換において低周波部分はそのまま保持され、高周波部分だけが位相シフトされることを主張します。つまり、包絡線 $A(t)$ はヒルベルト変換を「素通り」し、搬送波だけがコサインからサインに変わるのです。
したがって解析信号は、
$$ z(t) = A(t)\cos(2\pi f_c t) + jA(t)\sin(2\pi f_c t) = A(t)e^{j2\pi f_c t} $$
包絡線は解析信号の絶対値として取り出せます。
$$ \boxed{A(t) = |z(t)| = \sqrt{x(t)^2 + \hat{x}(t)^2}} $$
この関係は変調指数 $m$ が 1 を超える過変調の場合でも正しく動作します。従来の整流＋ローパスフィルタによる包絡線検波は過変調で歪みますが、ヒルベルト変換に基づく方法は $A(t)$ が負の値をとっても $|z(t)|$ として正しい振幅を返します。
ここまでで包絡線（信号の振幅情報）の抽出方法がわかりました。次は、信号の周波数情報を取り出す「瞬時周波数」について詳しく見ていきましょう。
瞬時周波数の抽出 — FM 信号への応用
瞬時周波数の定義
解析信号 $z(t) = A(t)e^{j\theta(t)}$ の瞬時位相 $\theta(t)$ を時間微分して $2\pi$ で割ったものが瞬時周波数です。
$$ f_i(t) = \frac{1}{2\pi}\frac{d\theta(t)}{dt} $$
瞬時周波数は「信号がその瞬間にどの周波数で振動しているか」を表す量です。単一周波数の正弦波 $\cos(2\pi f_0 t)$ であれば瞬時周波数は定数 $f_0$ ですが、FM 信号のように周波数が時間変化する信号では $f_i(t)$ も時間とともに変動します。
FM 信号への適用
FM（周波数変調）信号は次の形をしています。
$$ x_{\mathrm{FM}}(t) = A_c \cos\!\left[2\pi f_c t + 2\pi k_f \int_0^t s(\tau)\,d\tau\right] $$
ここで $k_f$ は周波数感度 [Hz/V]、$s(t)$ はベースバンド信号です。瞬時位相は、
$$ \theta(t) = 2\pi f_c t + 2\pi k_f \int_0^t s(\tau)\,d\tau $$
瞬時角周波数を求めるため、この瞬時位相を時間で微分します。積分の微分は被積分関数そのものになるので、次式を得ます。
$$ \frac{d\theta}{dt} = 2\pi f_c + 2\pi k_f \, s(t) $$
したがって瞬時周波数は、
$$ f_i(t) = f_c + k_f \, s(t) $$
つまり、FM 信号の瞬時周波数から搬送波周波数 $f_c$ を引けば、ベースバンド信号 $s(t)$ を $k_f$ 倍したものが得られます。これがヒルベルト変換を用いた FM 復調の原理です。
瞬時周波数の数値計算
実際の数値計算では、瞬時位相 $\theta(t) = \arctan[\hat{x}(t)/x(t)]$ の微分を求める必要があります。ここでいくつかの注意点があります。
位相のアンラッピング: $\arctan$ の値域は $(-\pi, \pi]$ に限られるため、位相が $\pm\pi$ を越えると不連続な飛びが生じます。これを補正する操作が位相アンラッピング（phase unwrapping）です。NumPy では numpy.unwrap 関数がこの処理を行います。
数値微分: 位相をアンラッピングした後、中心差分や勾配関数（numpy.gradient）で数値微分します。
$$ f_i[n] \approx \frac{\theta_{\text{unwrap}}[n+1] – \theta_{\text{unwrap}}[n-1]}{2 \cdot 2\pi \Delta t} $$
ここで $\Delta t = 1/f_s$ はサンプリング間隔です。
包絡線と瞬時周波数という 2 つの基本応用を理解したところで、次にヒルベルト変換の通信工学における重要な応用である SSB 変調について解説します。
SSB 変調への応用 — 位相法による単側波帯生成
DSB と SSB の違い
通常の AM 変調（DSB: 両側波帯）では、ベースバンド信号のスペクトルが搬送波の上側と下側に対称的にコピーされます。しかし、上側波帯（USB: Upper Sideband）と下側波帯（LSB: Lower Sideband）は同じ情報を担っているため、片方だけ送れば帯域幅を半分に節約できます。これが SSB（Single Sideband）変調です。
SSB 変調を実現する方法にはいくつかありますが、そのうちの一つが位相法（phasing method）であり、ヒルベルト変換を直接利用します。
位相法の原理
位相法では、ベースバンド信号 $s(t)$ とそのヒルベルト変換 $\hat{s}(t)$ の両方を用いて SSB 信号を生成します。
USB（上側波帯） を残す場合:
$$ \boxed{x_{\mathrm{USB}}(t) = s(t)\cos(2\pi f_c t) – \hat{s}(t)\sin(2\pi f_c t)} $$
LSB（下側波帯） を残す場合:
$$ \boxed{x_{\mathrm{LSB}}(t) = s(t)\cos(2\pi f_c t) + \hat{s}(t)\sin(2\pi f_c t)} $$
位相法の導出
なぜこの式で SSB が得られるのか、フーリエ変換を用いて確認しましょう。
DSB-SC（搬送波抑圧両側波帯）信号は $x_{\mathrm{DSB}}(t) = s(t)\cos(2\pi f_c t)$ であり、そのフーリエ変換は、
$$ X_{\mathrm{DSB}}(f) = \frac{1}{2}\left[S(f – f_c) + S(f + f_c)\right] $$
これは $S(f)$ を $\pm f_c$ にシフトしたもので、上側波帯と下側波帯の両方を含みます。
一方、$\hat{s}(t)\sin(2\pi f_c t)$ のフーリエ変換を計算しましょう。$\hat{s}(t)$ のフーリエ変換は $\hat{S}(f) = -j\,\mathrm{sgn}(f) \cdot S(f)$ であり、$\sin(2\pi f_c t) = [e^{j2\pi f_c t} – e^{-j2\pi f_c t}]/(2j)$ のフーリエ変換はインパルスの差です。これらの畳み込みを計算すると、
$$ \mathcal{F}\{\hat{s}(t)\sin(2\pi f_c t)\} = \frac{1}{2j}\left[\hat{S}(f – f_c) – \hat{S}(f + f_c)\right] $$
$\hat{S}(f) = -j\,\mathrm{sgn}(f) \cdot S(f)$ を代入して整理します。$f > f_c$ の領域（上側波帯）では $\mathrm{sgn}(f – f_c) = 1$ なので $\hat{S}(f – f_c) = -jS(f – f_c)$ となり、$f < -f_c$ の領域（$f + f_c < 0$）では $\mathrm{sgn}(f + f_c) = -1$ なので $\hat{S}(f + f_c) = jS(f + f_c)$ となります。
これらを $x_{\mathrm{USB}}(t) = s(t)\cos(2\pi f_c t) – \hat{s}(t)\sin(2\pi f_c t)$ に代入して周波数領域で整理すると、上側波帯成分は残り、下側波帯成分はキャンセルされることが確認できます。
直感的に言えば、ヒルベルト変換による 90° 位相シフトと搬送波の 90° 位相シフト（$\cos \to \sin$）が組み合わさることで、一方の側波帯では位相が揃って強め合い、他方では逆位相になって打ち消し合うのです。
位相法の実用上の課題
位相法は原理的にはシンプルですが、実用上はベースバンド信号の全帯域にわたって正確に 90° の位相シフトを実現する必要があります。アナログ回路では広帯域の位相シフタを作るのが難しいため、かつてはフィルタ法（SSB フィルタを使う方式）が主流でした。しかし、ディジタル信号処理では FFT を用いてヒルベルト変換を高精度に実装できるため、位相法がソフトウェア無線（SDR）などで広く使われるようになっています。
SSB 変調の位相法を理解できたところで、次にディジタル信号処理でヒルベルト変換をどう実装するかを詳しく見ていきましょう。
離散ヒルベルト変換
連続から離散へ
離散信号 $x[n]$（$n = 0, 1, \dots, N-1$）に対してヒルベルト変換を行うには、連続時間の理論を離散周波数領域に翻訳します。基本的なアイデアは以下の通りです。
- $x[n]$ の DFT（離散フーリエ変換）$X[k]$ を計算する
- 負の周波数成分をゼロにし、正の周波数成分を 2 倍にする
- 逆 DFT で時間領域に戻す
この操作で得られるのは解析信号 $z[n]$ であり、その虚部が離散ヒルベルト変換 $\hat{x}[n]$ です。
具体的なアルゴリズム
$N$ 点の実数信号 $x[n]$ に対して、DFT $X[k]$（$k = 0, 1, \dots, N-1$）を計算した後、以下の重み関数 $h[k]$ を乗じます。
$N$ が偶数の場合:
$$ h[k] = \begin{cases} 1 & k = 0 \text{ (DC)} \\ 2 & k = 1, 2, \dots, N/2 – 1 \text{ (正の周波数)} \\ 1 & k = N/2 \text{ (ナイキスト周波数)} \\ 0 & k = N/2 + 1, \dots, N-1 \text{ (負の周波数)} \end{cases} $$
$N$ が奇数の場合:
$$ h[k] = \begin{cases} 1 & k = 0 \text{ (DC)} \\ 2 & k = 1, 2, \dots, (N-1)/2 \text{ (正の周波数)} \\ 0 & k = (N+1)/2, \dots, N-1 \text{ (負の周波数)} \end{cases} $$
DC 成分（$k = 0$）とナイキスト成分（$k = N/2$、$N$ 偶数の場合）は正負の周波数の境界上にあるため、係数を 1（そのまま保持）とします。
解析信号の DFT は $Z[k] = h[k] \cdot X[k]$ であり、逆 DFT で時間領域に戻すと複素信号 $z[n]$ が得られます。
$$ z[n] = \mathrm{IDFT}\{h[k] \cdot X[k]\} $$
この操作は SciPy の scipy.signal.hilbert 関数が内部で行っている処理そのものです。ただし注意が必要なのは、SciPy の hilbert 関数はヒルベルト変換 $\hat{x}[n]$ ではなく解析信号 $z[n]$ を返すという点です。ヒルベルト変換そのものが必要な場合は、返り値の虚部を取り出します。
窓関数と端点効果
離散ヒルベルト変換を FFT で実装する場合、信号の両端での打ち切り効果（端点効果）に注意が必要です。FFT は信号を周期的とみなすため、信号の始まりと終わりが不連続だと、スペクトルに漏れ（リーケージ）が生じ、ヒルベルト変換の精度が悪化します。
対策としては、信号の両端にゼロパディングを十分に施す方法や、窓関数（ハニング窓やハミング窓など）を適用してから変換する方法があります。
理論の解説が一通り終わりました。ここからは Python を使って、ヒルベルト変換の各応用を実際に実装し、理論の正しさを視覚的に確認していきましょう。
Python 実装 1: ヒルベルト変換の基本と包絡線検出
まずは最も基本的な応用として、AM 変調信号からヒルベルト変換を用いて包絡線を抽出するコードを実装します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
# --- パラメータ設定 ---
fs = 10000 # サンプリング周波数 [Hz]
T = 0.1 # 信号長 [秒]
t = np.arange(0, T, 1/fs)
fc = 500 # 搬送波周波数 [Hz]
fm = 30 # ベースバンド周波数 [Hz]
m = 0.7 # 変調指数
# --- AM 変調信号の生成 ---
envelope_true = 1 + m * np.cos(2 * np.pi * fm * t) # 真の包絡線
x_am = envelope_true * np.cos(2 * np.pi * fc * t) # AM 信号
# --- ヒルベルト変換による包絡線検出 ---
z = hilbert(x_am) # 解析信号 (scipy は解析信号を返す)
envelope_hilbert = np.abs(z) # 包絡線 = |z(t)|
# --- 可視化 ---
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
axes[0].plot(t * 1000, x_am, color='steelblue', alpha=0.7, label='AM signal $x(t)$')
axes[0].plot(t * 1000, envelope_true, 'r--', linewidth=2, label='True envelope')
axes[0].plot(t * 1000, envelope_hilbert, 'k-', linewidth=2, label='Hilbert envelope')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('AM Signal and Envelope Detection via Hilbert Transform')
axes[0].legend(loc='upper right')
axes[0].grid(True, alpha=0.3)
# 推定誤差
error = envelope_hilbert - envelope_true
axes[1].plot(t * 1000, error * 1000, color='darkorange')
axes[1].set_xlabel('Time [ms]')
axes[1].set_ylabel('Error [×10⁻³]')
axes[1].set_title('Envelope Estimation Error')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
上のグラフから、いくつかの重要な点が読み取れます。
- ヒルベルト変換による包絡線（黒の実線）が真の包絡線（赤の破線）にほぼ完全に一致している — 搬送波周波数 $f_c = 500\,\mathrm{Hz}$ がベースバンド周波数 $f_m = 30\,\mathrm{Hz}$ に対して十分高いため、Bedrosian の定理の条件が満たされ、高精度な包絡線推定が実現しています。
- 誤差は $10^{-3}$ のオーダーで極めて小さい — 推定誤差のグラフを見ると、端点付近でやや誤差が大きくなっていますが、信号の中央部では無視できるレベルです。端点効果は FFT ベースの処理に特有のものであり、ゼロパディングで緩和できます。
- AM 信号の正・負の振動をまたいで包絡線が滑らかに取り出されている — 単純な整流（絶対値をとってローパスフィルタ）と異なり、ヒルベルト変換は瞬間的に正確な包絡線を与えます。
次に、瞬時周波数の計算を FM 信号に適用してみましょう。
Python 実装 2: 瞬時周波数の計算
FM 変調信号に対してヒルベルト変換を適用し、瞬時周波数を正しく推定できるかを検証します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
# --- パラメータ設定 ---
fs = 10000 # サンプリング周波数 [Hz]
T = 0.2 # 信号長 [秒]
t = np.arange(0, T, 1/fs)
fc = 500 # 搬送波周波数 [Hz]
fm = 5 # ベースバンド周波数 [Hz]
kf = 150 # 周波数感度 [Hz/V]
# --- FM 変調信号の生成 ---
baseband = np.cos(2 * np.pi * fm * t) # ベースバンド信号
phase_integral = (kf / fm) * np.sin(2 * np.pi * fm * t) # kf * ∫s(τ)dτ
x_fm = np.cos(2 * np.pi * fc * t + 2 * np.pi * phase_integral) # FM 信号
# 真の瞬時周波数
f_inst_true = fc + kf * baseband
# --- ヒルベルト変換による瞬時周波数推定 ---
z = hilbert(x_fm)
inst_phase = np.unwrap(np.angle(z)) # 位相アンラッピング
f_inst_est = np.gradient(inst_phase, 1/fs) / (2 * np.pi) # 数値微分
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
# FM 波形
axes[0].plot(t * 1000, x_fm, color='steelblue', alpha=0.7)
axes[0].set_ylabel('Amplitude')
axes[0].set_title('FM Signal')
axes[0].grid(True, alpha=0.3)
# 瞬時周波数の比較
axes[1].plot(t * 1000, f_inst_true, 'r--', linewidth=2, label='True $f_i(t)$')
axes[1].plot(t * 1000, f_inst_est, 'k-', linewidth=1.5, label='Estimated $f_i(t)$')
axes[1].set_ylabel('Frequency [Hz]')
axes[1].set_title('Instantaneous Frequency')
axes[1].legend(loc='upper right')
axes[1].grid(True, alpha=0.3)
# 推定誤差
error_f = f_inst_est - f_inst_true
axes[2].plot(t * 1000, error_f, color='darkorange')
axes[2].set_xlabel('Time [ms]')
axes[2].set_ylabel('Error [Hz]')
axes[2].set_title('Instantaneous Frequency Estimation Error')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフからは以下のことが読み取れます。
- FM 波形（上段）を見ると、搬送波の周波数が時間とともに変動している — ベースバンド信号が正のとき（瞬時周波数が高いとき）は波の間隔が密になり、負のとき（瞬時周波数が低いとき）は疎になっています。これは FM 変調の直感的な特徴そのものです。
- 推定された瞬時周波数（中段の黒実線）が真値（赤破線）と非常によく一致している — 搬送波周波数 $f_c = 500\,\mathrm{Hz}$ を中心に $\pm k_f = \pm 150\,\mathrm{Hz}$ の範囲で正弦的に変動する瞬時周波数が正確に追跡されています。
- 誤差（下段）は信号の中央部で 1 Hz 未満 — 端点付近を除けば、非常に高精度な推定が実現できています。数値微分と FFT の端点効果が誤差の主因です。
続いて、SSB 変調の位相法を Python で実装してみましょう。
Python 実装 3: SSB 信号生成（位相法）
ヒルベルト変換を用いた位相法で SSB 信号を生成し、そのスペクトルを DSB と比較して片側波帯だけが残ることを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
from scipy.fft import fft, fftfreq
# --- パラメータ設定 ---
fs = 10000 # サンプリング周波数 [Hz]
T = 0.5 # 信号長 [秒]
N = int(fs * T)
t = np.arange(N) / fs
fc = 1000 # 搬送波周波数 [Hz]
fm = 100 # ベースバンド周波数 [Hz]
# --- ベースバンド信号（単一トーン）---
s = np.cos(2 * np.pi * fm * t)
# --- ヒルベルト変換 ---
s_hilbert = np.imag(hilbert(s)) # s(t) のヒルベルト変換 = sin(2πfm*t)
scipy.signal.hilbert で解析信号を計算し、その虚部を取ることで $\cos(2\pi f_m t)$ に対する 90° 位相シフト信号 $\sin(2\pi f_m t)$ を得ています。この信号を搬送波と組み合わせて、DSB・USB・LSB の各信号を生成します。
# --- DSB-SC 信号 ---
x_dsb = s * np.cos(2 * np.pi * fc * t)
# --- SSB 信号（USB: 上側波帯）---
x_usb = s * np.cos(2 * np.pi * fc * t) - s_hilbert * np.sin(2 * np.pi * fc * t)
# --- SSB 信号（LSB: 下側波帯）---
x_lsb = s * np.cos(2 * np.pi * fc * t) + s_hilbert * np.sin(2 * np.pi * fc * t)
USB では $\hat{s}(t)\sin(2\pi f_c t)$ を引くことで下側波帯が打ち消され、LSB では足すことで上側波帯が打ち消されます。次にこれら3信号のスペクトルを計算・比較し、片側波帯だけが残ることを確認します。
# --- スペクトル計算 ---
def compute_spectrum(sig, fs):
"""片側振幅スペクトルを計算"""
N = len(sig)
X = fft(sig)
freqs = fftfreq(N, 1/fs)
# 正の周波数のみ
pos = freqs >= 0
magnitude = 2.0 / N * np.abs(X[pos])
return freqs[pos], magnitude
freq_dsb, mag_dsb = compute_spectrum(x_dsb, fs)
freq_usb, mag_usb = compute_spectrum(x_usb, fs)
freq_lsb, mag_lsb = compute_spectrum(x_lsb, fs)
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
axes[0].plot(freq_dsb, mag_dsb, color='steelblue')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('DSB-SC Spectrum')
axes[0].set_xlim(700, 1300)
axes[0].grid(True, alpha=0.3)
axes[1].plot(freq_usb, mag_usb, color='crimson')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('USB (Upper Sideband) Spectrum')
axes[1].set_xlim(700, 1300)
axes[1].grid(True, alpha=0.3)
axes[2].plot(freq_lsb, mag_lsb, color='forestgreen')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('LSB (Lower Sideband) Spectrum')
axes[2].set_xlim(700, 1300)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このスペクトル比較から、位相法による SSB 変調の原理が明確に確認できます。
- DSB-SC（上段）では $f_c \pm f_m = 1000 \pm 100\,\mathrm{Hz}$ の 2 つのピークが対称的に現れる — これが両側波帯であり、900 Hz（LSB）と 1100 Hz（USB）に等しい振幅のスペクトル線があります。
- USB（中段）では 1100 Hz のピークだけが残り、900 Hz のピークが消えている — ヒルベルト変換による位相法が下側波帯を正確にキャンセルしています。残ったスペクトル線の振幅は DSB の場合と同じであり、エネルギーが保存されていることがわかります。
- LSB（下段）では逆に 900 Hz のピークだけが残っている — 符号を反転するだけで上下の選択が切り替わるのが位相法の利点です。
単一トーンの例で原理を確認しましたが、実際の音声信号のように複数の周波数成分を含む信号に対しても、位相法は全帯域にわたってヒルベルト変換が正確であれば同様に動作します。
ここまでの 3 つの実装で主要な応用を確認できました。最後に、複合的な信号に対してヒルベルト変換を適用する総合的な例を見てみましょう。
Python 実装 4: チャープ信号の解析
実用的な例として、周波数が時間とともに線形に増加するチャープ信号（線形 FM 信号）に対してヒルベルト変換を適用し、包絡線と瞬時周波数の同時推定を行います。チャープ信号はレーダーのパルス圧縮や音響測定で広く使われています。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp
# --- パラメータ設定 ---
fs = 5000 # サンプリング周波数 [Hz]
T = 1.0 # 信号長 [秒]
t = np.arange(0, T, 1/fs)
f0 = 100 # 開始周波数 [Hz]
f1 = 800 # 終了周波数 [Hz]
# --- チャープ信号の生成（振幅にガウス窓を適用）---
gaussian_env = np.exp(-((t - T/2) ** 2) / (2 * 0.15**2)) # ガウス包絡線
x_chirp = gaussian_env * chirp(t, f0=f0, f1=f1, t1=T, method='linear')
ガウス窓を包絡線として適用することで信号の両端が滑らかにゼロへ減衰し、FFT ベースのヒルベルト変換における端点効果を緩和しています。次に解析信号を計算し、包絡線と瞬時周波数を同時に推定します。
# --- ヒルベルト変換 ---
z = hilbert(x_chirp)
envelope = np.abs(z)
inst_phase = np.unwrap(np.angle(z))
f_inst = np.gradient(inst_phase, 1/fs) / (2 * np.pi)
# 真の瞬時周波数（線形チャープ）
f_inst_true = f0 + (f1 - f0) * t / T
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 9))
# チャープ信号と包絡線
axes[0].plot(t, x_chirp, color='steelblue', alpha=0.6, label='Chirp signal')
axes[0].plot(t, envelope, 'r-', linewidth=2, label='Hilbert envelope')
axes[0].plot(t, gaussian_env, 'k--', linewidth=2, label='True envelope')
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Chirp Signal with Gaussian Envelope')
axes[0].legend(loc='upper right')
axes[0].grid(True, alpha=0.3)
# 瞬時周波数
axes[1].plot(t, f_inst_true, 'r--', linewidth=2, label='True $f_i(t)$')
axes[1].plot(t, f_inst, 'k-', linewidth=1.5, label='Estimated $f_i(t)$')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Frequency [Hz]')
axes[1].set_title('Instantaneous Frequency of Chirp Signal')
axes[1].legend(loc='upper left')
axes[1].grid(True, alpha=0.3)
# 時間-周波数表現（瞬時振幅で色付け）
axes[2].scatter(t, f_inst, c=envelope, cmap='inferno', s=1, alpha=0.8)
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Frequency [Hz]')
axes[2].set_title('Time-Frequency Representation (color = amplitude)')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このチャープ信号の解析結果からは、ヒルベルト変換の強力さが一目でわかります。
- 包絡線推定（上段） — ガウス窓で振幅変調されたチャープ信号に対して、ヒルベルト変換による包絡線（赤実線）が真のガウス包絡線（黒破線）にほぼ完全に一致しています。信号の両端（振幅がほぼゼロの領域）で若干の乖離が見られますが、これは SNR が低い領域では避けられない現象です。
- 瞬時周波数推定（中段） — 100 Hz から 800 Hz へ線形に増加する真の瞬時周波数を、推定値が忠実に追跡しています。包絡線の中心部（振幅が大きい領域）では推定精度が特に高く、端部では数値微分の不安定さから若干の振動が見られます。
- 時間-周波数表現（下段） — 瞬時周波数を縦軸、時刻を横軸、色で振幅を表現すると、スペクトログラムのような時間-周波数マップが得られます。チャープ信号の周波数が時間とともに線形に上昇する様子が鮮明に可視化されています。
ヒルベルト変換の重要な性質のまとめ
ここまでの議論を踏まえ、ヒルベルト変換の主要な数学的性質を整理しておきましょう。
線形性
ヒルベルト変換は線形変換です。すなわち、任意の信号 $x_1(t)$, $x_2(t)$ と定数 $a$, $b$ に対して、
$$ \mathcal{H}\{a \, x_1(t) + b \, x_2(t)\} = a \, \mathcal{H}\{x_1(t)\} + b \, \mathcal{H}\{x_2(t)\} $$
これは $H(f) = -j\,\mathrm{sgn}(f)$ との乗算が線形操作であることから明らかです。
二重適用（自己逆性）
ヒルベルト変換を 2 回適用すると、元の信号の符号が反転します。
$$ \mathcal{H}\{\mathcal{H}\{x(t)\}\} = -x(t) $$
周波数領域で確認すると、$H(f) \cdot H(f) = [-j\,\mathrm{sgn}(f)]^2 = j^2 \cdot \mathrm{sgn}^2(f) = -1$（$f \neq 0$）となるので、2 回適用すると $-1$ が掛かり符号が反転します。言い換えれば、ヒルベルト変換を 4 回適用すると元に戻ります（$H^4 = 1$）。これは 90° の位相シフトを 4 回繰り返すと $360°$ に達するという直感と合致します。
エネルギー保存（パーセバルの等式）
ヒルベルト変換は信号のエネルギーを変えません。
$$ \int_{-\infty}^{\infty} |\hat{x}(t)|^2 \, dt = \int_{-\infty}^{\infty} |x(t)|^2 \, dt $$
これは $|H(f)| = 1$（$f \neq 0$）であるため、パワースペクトル密度が保存されることから直ちに導かれます。
直交性
元の信号 $x(t)$ とそのヒルベルト変換 $\hat{x}(t)$ は直交します。
$$ \int_{-\infty}^{\infty} x(t) \, \hat{x}(t) \, dt = 0 $$
これは解析信号の構成において「実部と虚部が独立な情報を担う」ことの数学的表現です。
DC 成分の除去
定数信号（DC 成分）のヒルベルト変換はゼロです。
$$ \mathcal{H}\{c\} = 0 \quad (c = \text{定数}) $$
$H(0) = 0$ なので、DC 成分（$f = 0$）はヒルベルト変換で消滅します。このため、包絡線の推定に先立って信号の DC オフセットを除去しておくのが良い実践です。
これらの性質は、ヒルベルト変換を信号処理の各場面で応用する際の理論的基盤となります。
ヒルベルト変換の適用上の注意点
理論が綺麗に成り立つ条件と、実際の信号処理で遭遇する問題点をまとめておきましょう。
狭帯域信号の仮定
包絡線と瞬時周波数の分離が意味を持つのは、信号が狭帯域（narrowband）である場合に限られます。広帯域信号では包絡線の概念自体が曖昧になり、瞬時周波数が負の値をとることもあります。一般に、信号の帯域幅が搬送波周波数に比べて十分小さい（$B \ll f_c$）とき、ヒルベルト変換による振幅・位相の分離は良好に機能します。
端点効果
FFT ベースの離散ヒルベルト変換では、信号の両端に偽の振動が現れることがあります。これは FFT が信号を周期的とみなすために生じる現象であり、対策としてゼロパディングや窓関数の適用が有効です。また、リアルタイム処理では FIR 近似のヒルベルト変換フィルタを用いることがあり、フィルタ長に応じた遅延が生じます。
数値精度
瞬時周波数の計算では位相の数値微分を行うため、ノイズが増幅されやすいという問題があります。サンプリング周波数に対して信号周波数が十分低い（オーバーサンプリングされている）場合はこの影響は小さいですが、ナイキスト周波数に近い成分を含む信号では精度が悪化します。移動平均などの平滑化や、Savitzky-Golay フィルタによる微分の安定化が実用上有効です。
多成分信号
信号が複数の独立な周波数成分を含む場合（例: 2 つの異なる周波数の正弦波の重ね合わせ）、解析信号から得られる包絡線や瞬時周波数は各成分の合成的な効果を反映し、物理的に解釈が難しくなることがあります。このような場合には、経験的モード分解（EMD）などでまず信号を単一成分に分離してからヒルベルト変換を適用する手法（Hilbert-Huang Transform）が用いられます。
まとめ
本記事では、ヒルベルト変換の理論を基礎から体系的に解説しました。
- ヒルベルト変換は「全帯域 90° 位相シフタ」である — 信号に含まれるすべての周波数成分の位相を一律に $\pi/2$ ずらし、振幅は変えない線形変換です
- 時間領域ではコーシーの主値積分 $\hat{x}(t) = \frac{1}{\pi}\mathrm{p.v.}\int \frac{x(\tau)}{t-\tau}d\tau$、周波数領域では $-j\,\mathrm{sgn}(f)$ との乗算として表現されます
- 解析信号 $z(t) = x(t) + j\hat{x}(t)$ を構成することで、負の周波数成分を除去し、包絡線と瞬時位相を分離できます
- 包絡線 $A(t) = |z(t)|$ は AM 信号の復調に、瞬時周波数 $f_i(t) = \frac{1}{2\pi}\frac{d\theta}{dt}$ は FM 信号の復調に直接利用できます
- SSB 変調の位相法はヒルベルト変換を用いて片側波帯だけを生成する方式であり、帯域幅を半分に節約できます
- 離散ヒルベルト変換は FFT を用いて効率的に実装でき、SciPy の
scipy.signal.hilbertで簡単に利用できます
ヒルベルト変換は信号処理の基盤的なツールであり、通信工学、振動解析、音響処理、レーダー信号処理など、幅広い分野で重要な役割を果たしています。特に解析信号の概念は、I/Q 信号処理やソフトウェア無線の基礎として、現代のディジタル通信システムに不可欠です。
次のステップとして、以下の記事も参考にしてください。