電話の向こうの声が二重に聞こえる「やまびこ（エコー）」を、信号だけから検出して消すにはどうすればよいでしょうか。あるいは、録音された母音「あー」から、その人が今どれくらいの高さ（ピッチ）で発声しているかを、自動で数値として取り出すにはどうすればよいでしょうか。どちらの問題も、一見すると性質がまったく違うように見えますが、実は 「ある信号ともう一つの信号の畳み込み（あるいは周期的な繰り返し）を、足し算に変えて分離する」 という同じアイデアで解けます。そのための道具が ケプストラム（cepstrum） です。
ケプストラムは「スペクトルのスペクトル」、つまり対数スペクトルをもう一度フーリエ変換して得られる量です。最初は奇妙に聞こえますが、この一見回りくどい操作には明確な理由があります。スペクトル上で「等間隔に並んだ山（ハーモニクス）」や「リップル（さざ波）」として現れる周期的構造は、もう一度フーリエ変換するとはっきりとしたピークになって浮かび上がるのです。
ケプストラムは次のような場面で活躍します。
- 音声処理のピッチ抽出: 声帯振動の基本周期を、母音波形のケプストラムのピーク位置から正確に検出する（ボコーダー、音声認識の前処理、歌声分析）。
- 反響・多重反射の検出と除去: レーダーやソナー、室内音響、地震波探査において、エコー（時間遅延を伴う反射波）の遅延時間を推定し、ホモモルフィック・デコンボリューションで分離する。
本記事の内容
- なぜ「スペクトルのスペクトル」が役に立つのか（畳み込み → 対数 → 加算の直感）
- ケプストラム $c[n] = \mathrm{IDFT}\{\log|X(e^{j\omega})|\}$ の定義とケフレンシ軸
- ピーク位置が基本周期・エコー遅延に対応する仕組みの導出
- 実ケプストラムと複素ケプストラムの違い
- リフタリング（ケフレンシ領域のフィルタリング）による包絡と微細構造の分離
- Pythonによる合成エコー信号・母音波形の解析とピッチ／遅延ピークの検出
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特に、フーリエ変換が「信号を周波数成分に分解する」という基本と、畳み込みが周波数領域で積になるという畳み込み定理を理解していると、ケプストラムの本質がすっと頭に入ります。
ケプストラムとは — 「スペクトルのスペクトル」という発想
まず、ケプストラムという名前自体が、信号処理の遊び心を表しています。cepstrum は spectrum（スペクトル）の最初の4文字 “spec” を逆に綴った “ceps” から作られた造語です。同じように、周波数を表す frequency をひっくり返した quefrency（ケフレンシ）、フィルタリングを表す filtering をひっくり返した liftering（リフタリング） という言葉も使われます。なぜこんな言葉遊びをするのかというと、ケプストラムは「フーリエ変換をもう一回かける」操作によって、周波数領域でも時間領域でもない、第三の領域へ信号を移すからです。その新しい領域での「軸」や「操作」には、もとの名前を流用しつつ少しずらした名前を与えたのです。
直感をつかむために、簡単な場面を考えましょう。山に向かって「ヤッホー」と叫ぶと、少し遅れて「ヤッホー」というやまびこが返ってきます。マイクで録音した音 $y(t)$ は、もとの声 $x(t)$ と、それが遅れて小さくなったコピー $\alpha\, x(t-\tau)$ の和です。
$$ y(t) = x(t) + \alpha\, x(t – \tau) $$
この $y(t)$ をそのままフーリエ変換しても、もとの声のスペクトルにエコーの影響が「掛け算」として混ざり込み、複雑なスペクトルになります。後で見るように、このときスペクトルには周期 $1/\tau$ の細かいさざ波（リップル）が乗ります。人間の目には「ギザギザしたスペクトル」にしか見えませんが、このギザギザ自体が一つの「周期信号」だと考えると、それをもう一度フーリエ変換すれば、ギザギザの周期に対応するピークが現れるはずです。そのピークの位置こそが、エコーの遅延時間 $\tau$ なのです。
これがケプストラムの核心です。すなわち、
- 信号をフーリエ変換してスペクトルを得る。
- その振幅スペクトルの 対数 を取る（掛け算を足し算に変えるため）。
- それをもう一度フーリエ変換（または逆フーリエ変換）する。
この3ステップで、もとの信号に隠れていた「周期的構造の周期」が、新しい軸（ケフレンシ軸）の上のピークとして読み取れるようになります。ケフレンシ軸の単位は、フーリエ変換をもう一度かけた結果なので、周波数の逆、すなわち 時間 の次元を持ちます。エコー遅延もピッチ周期も、まさに時間の量なので、ケプストラムのピーク位置を読めば直接それらが得られるわけです。
ここまでで「なぜスペクトルのスペクトルなのか」というアイデアの輪郭が見えました。次に、このアイデアを支える最も重要な原理 — 畳み込みが対数を経て加算に変わる仕組み — を数式で丁寧にたどります。
畳み込みが対数スペクトルで加算になる原理
ケプストラムが機能する理論的な土台は、ホモモルフィック（準同型）処理 という考え方にあります。難しそうな名前ですが、やっていることは「掛け算という演算を、対数を通して足し算という演算に翻訳する」だけです。$\log(ab) = \log a + \log b$ という、誰もが知っている対数の性質が主役です。
畳み込みモデルから出発する
声道による発声を例に考えましょう。人間の声は、声帯の振動が作る「励振源」 $e[n]$（パルス列のような周期信号）が、声道（口や鼻の空洞）という線形フィルタ $h[n]$ を通って出てきたものとモデル化できます。出力信号 $x[n]$ は両者の畳み込みです。
$$ x[n] = e[n] * h[n] $$
エコーの場合も同じ形に書けます。もとの音 $s[n]$ に対し、直接波と遅延した反射波からなるインパルス応答 $h[n] = \delta[n] + \alpha\,\delta[n-D]$ を畳み込むと、
$$ x[n] = s[n] * \big(\delta[n] + \alpha\,\delta[n-D]\big) = s[n] + \alpha\, s[n-D] $$
となり、まさにエコー信号が得られます。ここで $D$ はサンプル単位の遅延です。いずれの場合も「分離したい2つの成分の畳み込み」という共通構造を持っています。
フーリエ変換で畳み込みを積にする
畳み込み定理により、時間領域の畳み込みは周波数領域の掛け算になります。$x[n] = e[n] * h[n]$ の離散時間フーリエ変換（DTFT）を取ると、
$$ X(e^{j\omega}) = E(e^{j\omega}) \cdot H(e^{j\omega}) $$
ここで $E$ と $H$ はそれぞれ励振源と声道のスペクトルです。この時点で畳み込みは積に変わりましたが、まだ2つの成分は「掛け算」で絡み合ったままです。掛け算のままでは、片方だけを取り出す（フィルタリングする）ことができません。フィルタリングは本質的に「足し算で混ざった成分を、線形な操作で選り分ける」操作だからです。
対数で積を和にする
そこで対数の出番です。複素数の対数を取ると、絶対値の対数と偏角に分かれます。
$$ \log X(e^{j\omega}) = \log\big(E(e^{j\omega}) \cdot H(e^{j\omega})\big) = \log E(e^{j\omega}) + \log H(e^{j\omega}) $$
絶対値だけに注目すれば、
$$ \log |X(e^{j\omega})| = \log |E(e^{j\omega})| + \log |H(e^{j\omega})| $$
となります。ここが決定的な一歩です。掛け算で絡み合っていた2成分が、対数を取ることで「足し算」になりました。 足し算で重なっているものは、それぞれが異なる場所（後で見るケフレンシ軸上の異なる位置）に分布していれば、線形操作で分離できます。これがホモモルフィック処理の心臓部です。掛け算という「非線形に見える混ざり方」を、対数という橋渡しによって「線形に扱える足し算」へ変換したのです。
もう一度フーリエ変換する理由
最後に、なぜこの対数スペクトルをもう一度フーリエ変換するのでしょうか。理由は、励振源と声道では、対数スペクトル上での「変化の速さ」がまったく違うからです。
- 声道スペクトル $\log|H|$ は、周波数に対してゆっくりと変化する滑らかな曲線（包絡、エンベロープ）です。声道の共鳴（フォルマント）がなだらかな山として現れます。
- 励振源スペクトル $\log|E|$ は、基本周波数の整数倍に細かいピークが等間隔に並んだ、激しく振動する成分です（ハーモニクス構造）。
「周波数軸上でゆっくり変化する成分」と「速く振動する成分」が足し算で重なっている — これはまさに、通常の信号処理で「低周波成分と高周波成分が混ざった信号」を扱う状況とそっくりです。そういう信号は、フーリエ変換して周波数ごとに分けてやれば分離できます。同じことを対数スペクトルに対してやるのです。対数スペクトルを $\omega$ の関数とみなしてフーリエ変換すると、
- ゆっくり変化する声道成分は 低ケフレンシ（原点近く）に集まり、
- 速く振動する励振成分は 高ケフレンシ（周期に対応する位置）にピークを作ります。
こうして、もともと畳み込みで分かちがたく結びついていた声道と励振源が、ケフレンシ軸という新しい軸の上では別々の場所に分かれて見えるようになります。この分離こそがケプストラムのご利益です。
原理が見えたところで、次はケプストラムを正式に定義し、ケフレンシ軸の意味を厳密に詰めていきます。
ケプストラムの数学的定義とケフレンシ軸
ここまでの3ステップ（フーリエ変換 → 対数 → フーリエ変換）を、きちんと式で書き下します。離散信号 $x[n]$ に対する 実ケプストラム（real cepstrum） $c[n]$ は次のように定義されます。
$$ \begin{equation} c[n] = \frac{1}{2\pi}\int_{-\pi}^{\pi} \log\big|X(e^{j\omega})\big|\, e^{j\omega n}\, d\omega \end{equation} $$
ここで $X(e^{j\omega}) = \sum_n x[n] e^{-j\omega n}$ は $x[n]$ のDTFTです。式 $(1)$ の右辺は、関数 $\log|X(e^{j\omega})|$ を $\omega$ の信号とみなしたときの 逆フーリエ変換 にほかなりません。つまり実ケプストラムは「対数振幅スペクトルの逆フーリエ変換」です。
実際の計算では、$N$ 点のDFTを使って次のように離散化します。
$$ \begin{equation} c[n] = \frac{1}{N}\sum_{k=0}^{N-1} \log\big|X[k]\big|\, e^{j 2\pi k n / N}, \qquad X[k] = \sum_{m=0}^{N-1} x[m]\, e^{-j 2\pi k m / N} \end{equation} $$
すなわち、$c[n] = \mathrm{IDFT}\big\{\log|\mathrm{DFT}\{x[n]\}|\big\}$ です。$\log|X[k]|$ は実数列なので、その逆DFTである $c[n]$ も（対称性から）実数列になります。これが「実ケプストラム」と呼ばれる理由です。
ケフレンシ軸 — 変数 $n$ は時間である
ここで最も大切なのは、ケプストラム $c[n]$ の独立変数 $n$ が何を意味するかです。$c[n]$ は「$\omega$ の関数を逆フーリエ変換した結果」なので、$n$ は $\omega$（角周波数 [rad/sample]）に共役な変数、すなわち 時間（サンプル）の次元 を持ちます。この軸を ケフレンシ軸（quefrency axis） と呼びます。
混乱しやすいので整理しましょう。
| 領域 | 横軸 | 単位 | 得られた経緯 |
|---|---|---|---|
| 時間領域 | $n$ | サンプル（秒） | もとの信号 $x[n]$ |
| 周波数領域 | $\omega$ または $k$ | rad/sample（Hz） | $x[n]$ をフーリエ変換 |
| ケプストラム領域 | $n$（ケフレンシ） | サンプル（秒） | $\log|X|$ をフーリエ変換 |
ケフレンシは時間の次元を持つにもかかわらず、もとの時間領域とは別物です。ケフレンシ軸上の位置 $\tau$（秒）は、「対数スペクトルが周波数軸上で $1/\tau$ [Hz] の周期で振動する成分」に対応します。たとえば母音の基本周波数が $f_0 = 125$ Hz であれば、ハーモニクスは周波数軸上で 125 Hz 間隔に並びます。この「125 Hz 間隔の振動」は、ケフレンシ軸上では $\tau = 1/f_0 = 8$ ms の位置にピークを作ります。この $8$ ms はちょうど声帯振動の基本周期です。
ケフレンシ軸の値とサンプル番号・物理時間の対応は次の通りです。サンプリング周波数を $f_s$ とすると、ケフレンシ index $n$ に対応する物理時間（クエフレンシ）は
$$ \tau_n = \frac{n}{f_s} \quad [\text{s}] $$
であり、その位置のピークは周波数領域での周期 $f_s/n$ [Hz] に対応します。
ピーク位置が基本周期・エコー遅延に対応する仕組み
なぜピークが「周期」の位置に立つのかを、最も簡単なエコーモデルで導出してみましょう。直接波1つと、遅延 $D$ サンプル・減衰 $\alpha$（$0<\alpha<1$）の反射波1つからなる信号
$$ x[n] = s[n] + \alpha\, s[n-D] $$
を考えます。インパルス応答は $h[n]=\delta[n]+\alpha\,\delta[n-D]$ で、そのDTFTは
$$ H(e^{j\omega}) = 1 + \alpha\, e^{-j\omega D} $$
です。$x = s*h$ なので $X = S\cdot H$、よって対数振幅は
$$ \log|X(e^{j\omega})| = \log|S(e^{j\omega})| + \log|H(e^{j\omega})| $$
と加算で分かれます。エコー由来の項 $\log|H|$ を詳しく見ます。$|H|^2 = H H^*$ を計算すると、
$$ |H(e^{j\omega})|^2 = (1+\alpha e^{-j\omega D})(1+\alpha e^{j\omega D}) = 1 + \alpha^2 + 2\alpha\cos(\omega D) $$
となります。最後の項 $\cos(\omega D)$ がポイントです。$\log|H|^2$ は $\omega$ について 周期 $2\pi/D$ で振動する 関数になっているのです。すなわち、エコーは対数スペクトルに「周期 $2\pi/D$ のさざ波（リップル）」を刻みます。
この振動成分を、もう一度フーリエ変換すれば、その「周期」に対応するケフレンシ位置 $n=D$ にピークが立ちます。これをきちんと示すため、$\alpha$ が小さいとして対数を級数展開します。$\log(1+u) = u – \tfrac{u^2}{2} + \tfrac{u^3}{3} – \cdots$ を $u = \alpha e^{-j\omega D}$ について使うと、
$$ \log H(e^{j\omega}) = \log(1+\alpha e^{-j\omega D}) = \alpha e^{-j\omega D} – \frac{\alpha^2}{2} e^{-j2\omega D} + \frac{\alpha^3}{3} e^{-j3\omega D} – \cdots $$
この式の各項を見てください。$e^{-j\omega D}$ は逆フーリエ変換すると $\delta[n-D]$、$e^{-j2\omega D}$ は $\delta[n-2D]$ になります。つまり複素ケプストラム（後述）で見れば、エコー成分は
$$ \hat{h}[n] = \alpha\,\delta[n-D] – \frac{\alpha^2}{2}\delta[n-2D] + \frac{\alpha^3}{3}\delta[n-3D] – \cdots $$
となり、ケフレンシ $n=D, 2D, 3D, \dots$ に、振幅が急速に減衰しながら並ぶピーク列を作ります。最も大きいピークは $n=D$、すなわち エコー遅延そのもの の位置に現れます。これが「ケプストラムのピーク位置がエコー遅延に等しい」ことの証明です。
母音のピッチ検出もまったく同じ理屈です。声帯の励振源 $e[n]$ は周期 $P$（基本周期）のパルス列なので、そのスペクトル $|E|$ は基本周波数間隔のハーモニクス（くし状のピーク）を持ちます。対数スペクトルでこのくし状構造は「周期 $2\pi/P$ の振動」となり、ケプストラムでは $n=P$ にピークを作ります。$P$ をサンプル数から秒に直し、その逆数を取れば基本周波数 $f_0 = f_s/P$ が得られます。
導出から、ケプストラムが「周期構造の周期」を直接ピーク位置として取り出す装置であることがはっきりしました。次に、ここで顔を出した「複素ケプストラム」と、最初に定義した「実ケプストラム」の違いを整理します。
実ケプストラムと複素ケプストラム
ケプストラムには大きく分けて2種類あります。式 $(1)$ で定義したのが 実ケプストラム で、対数の「振幅だけ」を使いました。もう一つの 複素ケプストラム（complex cepstrum） は、対数の「振幅と位相の両方」を使います。
複素ケプストラムの定義
複素ケプストラム $\hat{x}[n]$ は、複素対数（複素数の対数）を取ってから逆フーリエ変換します。
$$ \begin{equation} \hat{x}[n] = \frac{1}{2\pi}\int_{-\pi}^{\pi} \log X(e^{j\omega})\, e^{j\omega n}\, d\omega, \qquad \log X(e^{j\omega}) = \log|X(e^{j\omega})| + j\,\angle X(e^{j\omega}) \end{equation} $$
複素対数の実部は対数振幅 $\log|X|$、虚部は位相 $\angle X$ です。実ケプストラムが振幅だけを使うのに対し、複素ケプストラムは位相情報も保持します。
両者の関係は明快です。複素対数の実部が対数振幅、その逆フーリエ変換の「偶対称部分」が実ケプストラムに一致するため、
$$ c[n] = \frac{\hat{x}[n] + \hat{x}[-n]}{2} $$
という関係が成り立ちます。逆に言うと、実ケプストラムは複素ケプストラムから位相情報を捨てて偶対称化したものです。
使い分け — 検出か、再構成か
2つのケプストラムは目的によって使い分けます。
- 実ケプストラム: 計算が簡単で安定です（位相のアンラップが不要）。ピッチ検出やエコー遅延の 検出 のように「ピーク位置さえわかればよい」用途には実ケプストラムで十分です。本記事のピッチ検出・遅延検出はすべて実ケプストラムで行います。
- 複素ケプストラム: 位相を保持しているため、ケフレンシ領域で操作した後に 元の時間信号へ逆変換（再構成） できます。エコーを実際に「除去した波形」を取り出したい、声道と励振源を分離して片方だけを合成し直したい、といった ホモモルフィック・デコンボリューション には複素ケプストラムが必須です。
複素ケプストラムの計算で技術的に厄介なのは、位相 $\angle X(e^{j\omega})$ を $\omega$ について連続になるよう アンラップ（unwrap） する必要がある点です。位相は本来 $2\pi$ の不定性を持つため、$-\pi \sim \pi$ に折り返された値をそのまま使うと不連続が生じ、逆変換が壊れます。NumPy の np.unwrap を使えばこの処理は容易です。なぜアンラップがそこまで重要かというと、複素対数の虚部である位相が不連続だと、それを逆フーリエ変換した結果に高ケフレンシの偽の振動が大量に混入し、本来見たい成分のピークが埋もれてしまうからです。さらに、信号全体に純粋な遅延（線形位相成分）が含まれている場合は、その線形成分をあらかじめ差し引いておく必要もあります。線形位相はケフレンシ軸上で原点付近に巨大なインパルス的成分を生み、解析を乱すためです。これらの前処理を丁寧に行えば、複素ケプストラムは元信号を完全に復元できる可逆な表現になります。
実ケプストラムと複素ケプストラムの関係をもう少し具体的に見ておきましょう。実ケプストラムは振幅情報だけから作られるため、もとの信号の「位相」を捨てています。位相を捨てるということは、波形のタイミング情報を失うということです。たとえば、あるスペクトル振幅を持つ信号は無数にあり（位相の選び方は自由）、そのどれもが同じ実ケプストラムを与えます。それでもピッチやエコー遅延の検出が成立するのは、これらの「周期構造」が振幅スペクトルのリップルとして現れ、位相を見なくても周期そのものは振幅情報だけから読み取れるからです。逆に、エコーを実際に消した波形を合成し直したい場合には、どの瞬間にどの成分が来るかというタイミング、すなわち位相が不可欠になります。だからこそ再構成には複素ケプストラムが要るのです。「何を知りたいか」によって、捨ててよい情報と保持すべき情報が変わる — この見極めがケプストラム活用の勘所です。
最小位相と因果性
複素ケプストラムには美しい性質があります。$x[n]$ が 最小位相系（極・零点がすべて単位円内）であれば、複素ケプストラム $\hat{x}[n]$ は 因果的（$n<0$ でゼロ）になります。逆に最大位相成分は $n<0$ 側に現れます。先ほどのエコー級数展開で $\hat{h}[n]$ が $n=D,2D,\dots>0$ にのみ現れたのは、$|\alpha|<1$ で $h[n]$ が最小位相だったからです。この因果性により、ケフレンシ軸の正負を使って最小位相成分と最大位相成分を分離する、といった高度な処理も可能になります。
実ケプストラムと複素ケプストラムの役割分担が見えました。次は、この新しい領域でフィルタリングを行う「リフタリング」を見ていきます。これがケプストラムを「検出器」から「分離器」へと格上げします。
リフタリング — ケフレンシ領域のフィルタリング
通常の信号処理では、周波数領域で不要な帯域をゼロにする「フィルタリング」を行います。ケプストラムの世界では、同じことを ケフレンシ領域 で行います。これを言葉遊びで リフタリング（liftering） と呼びます。ケフレンシ軸上で特定の範囲を残したり消したりする窓関数を リフター（lifter） と言います。
なぜリフタリングが意味を持つのか、もう一度思い出しましょう。対数スペクトル上では、声道（包絡）はゆっくり変化する成分、励振源（ハーモニクス）は速く振動する成分でした。ケプストラムに移すと、前者は 低ケフレンシ（原点近く）に、後者は 高ケフレンシ（基本周期の位置）に分かれて存在します。したがって、
- ローパスリフター（low-time lifter）: 低ケフレンシだけを残す（高ケフレンシをゼロにする）と、滑らかな声道スペクトル包絡 $\log|H|$ だけが残ります。これを逆変換すればスペクトル包絡が得られ、フォルマント（共鳴周波数）の推定に使えます。
- ハイパスリフター（high-time lifter）: 高ケフレンシだけを残すと、ハーモニクス構造すなわち励振源の情報だけが残ります。ピッチ検出に使えます。
数式で書くと、リフター窓 $\ell[n]$ をケプストラムに掛けて、
$$ c_{\text{lift}}[n] = \ell[n]\, c[n] $$
とし、これをフーリエ変換し直せば、リフタリング後の対数スペクトル
$$ \log|X_{\text{lift}}(e^{j\omega})| = \sum_n c_{\text{lift}}[n]\, e^{-j\omega n} $$
が得られます。低ケフレンシ通過リフター $\ell[n]$（例えば $|n|\le L$ で 1、それ以外 0）を使えば、$\log|X_{\text{lift}}|$ は滑らかな包絡だけになります。
ここで重要な注意があります。実ケプストラムは対称（偶関数的）なので、ケフレンシ index は $0$ から $N-1$ まで並び、$n$ と $N-n$ が対称な対を成します。したがってローパスリフターは「原点近傍の両側」、すなわち $n \le L$ と $n \ge N-L$ の両方を残す必要があります。この点を実装で正しく扱わないと、得られる包絡が歪んでしまいます。
リフタリングを使う際の実践的な注意も挙げておきます。リフターのカットオフ位置をどこに置くかは、分離したい2成分のケフレンシ分布が重ならない位置に取るのが原則です。声道包絡は低ケフレンシ側に集中し、ピッチ周期はそれより高いケフレンシ（本記事の母音では基本周期 8 ms）に現れるので、その間にカットオフを置きます。注意したいのは、包絡そのものも単一のケフレンシ点ではなくある幅を持つことです。フォルマントが複数あると、その間隔に対応するさざ波が包絡側のケフレンシをやや高い方まで広げます。カットオフが低すぎる（後の実装で見るように 2 ms 程度だと）と、近接する2フォルマントが1つの山に融合して分離できず、逆に高すぎるとハーモニクスの一部が包絡側に漏れ込みます。本記事の母音（F1=730, F2=1090 Hz、ピッチ周期 8 ms）では 4 ms 前後が、2フォルマントを分離しつつハーモニクスを排除できる良い妥協点でした。また、リフター窓を矩形（急に 0/1 を切り替える）にすると、再変換した対数スペクトルにリンギング（ギブス現象に似たさざ波）が出ることがあります。これを避けたいときは、矩形ではなく緩やかに減衰する窓（指数窓やハミング状のリフター）を使うと包絡がより滑らかになります。MFCC で用いられる「正弦リフター」も、このリンギング抑制と高ケフレンシ係数の強調を兼ねた工夫の一つです。
ケプストラムの理論（定義・ケフレンシ・実/複素・リフタリング）が一通りそろいました。ここからは Python を使って、合成エコー信号と母音波形に対して実際にケプストラムを計算し、ピーク検出とリフタリングによる分離を体験します。
Pythonでの実装1: ケプストラムの基本計算とエコー遅延の検出
まずは最も理解しやすいエコー検出から始めます。ここでは検出対象を広帯域の合成音（軽く色付けした雑音バースト）とし、人工的にエコー（遅延 5 ms、減衰 0.6）を加え、その実ケプストラムを計算してピークが遅延位置に立つことを確認します。あえて強い周期性を持たない信号を使うのは、エコー検出の本質である「インパルス応答 $\delta[n]+\alpha\delta[n-D]$ が刻むリップル」だけを取り出し、もとの音自身の周期構造がケフレンシ軸上で混信するのを避けるためです（声帯振動のような強い周期は、次節のピッチ検出で正面から扱います）。
最初に、ケプストラムを計算する関数を用意します。これは本記事を通して使い回す中核の関数です。
import numpy as np
import matplotlib.pyplot as plt
def real_cepstrum(x, n_fft=None):
"""実ケプストラム c[n] = IDFT{ log|DFT{x}| } を計算する。
Parameters
----------
x : 1次元配列（解析対象の信号）
n_fft : FFT点数（Noneなら信号長を使用）
Returns
-------
c : 実ケプストラム（実数配列）
"""
if n_fft is None:
n_fft = len(x)
# スペクトルの振幅。log(0)を避けるため微小値を加える
spectrum = np.fft.fft(x, n_fft)
log_mag = np.log(np.abs(spectrum) + 1e-12)
# 対数振幅スペクトルを逆フーリエ変換 → 実部を取る
c = np.fft.ifft(log_mag).real
return c
この関数は定義式 $(2)$ をそのまま実装したものです。np.abs で振幅を取り、np.log で対数にし、np.fft.ifft で逆変換しています。1e-12 は対数が発散しないための保険です。出力の虚部はゼロに近いはずなので .real で実部を取り出します。
次に、エコー付きの合成信号を作ってケプストラムを計算します。
import numpy as np
import matplotlib.pyplot as plt
fs = 8000 # サンプリング周波数 [Hz]
dur = 0.064 # 信号長 [s]（512サンプル）
t = np.arange(int(fs * dur)) / fs
# もとの音: 広帯域の合成音（軽く平滑化した雑音バースト）
# 強い周期性を持たせないことで、エコー由来のリップルだけを際立たせる
rng = np.random.default_rng(1)
s = rng.standard_normal(len(t))
s = np.convolve(s, np.ones(5)/5, mode='same') # 軽く平滑化してスペクトルに色付け
# 窓をかけて端の不連続を抑える
s *= np.hanning(len(s))
# エコーを付加: 遅延 5ms, 減衰 0.6
delay_ms = 5.0
D = int(fs * delay_ms / 1000) # 遅延サンプル数
alpha = 0.6
x = s.copy()
x[D:] += alpha * s[:-D] # x[n] = s[n] + alpha*s[n-D]
print(f"遅延 {delay_ms} ms = {D} サンプル")
ここでは平滑化した雑音バーストという広帯域の波形 s を作り、それを 5 ms（= 40 サンプル）遅らせて 0.6 倍したものを足し込み、エコー信号 x を作っています。出力から、5 ms がちょうど 40 サンプルに対応することが確認できます。続いてケプストラムを計算し、可視化します。
import numpy as np
import matplotlib.pyplot as plt
# 上のセルで定義した real_cepstrum, x, fs, D を使う
c = real_cepstrum(x)
quefrency = np.arange(len(c)) / fs * 1000 # ケフレンシ [ms]
half = len(c) // 2 # 実ケプストラムは対称なので前半だけ見る
fig, axes = plt.subplots(3, 1, figsize=(10, 9))
# (1) 時間波形
axes[0].plot(np.arange(len(x))/fs*1000, x, lw=0.8)
axes[0].set_xlabel('Time [ms]'); axes[0].set_ylabel('Amplitude')
axes[0].set_title('Echo signal x[n] = s[n] + 0.6 s[n-D]')
axes[0].grid(alpha=0.3)
# (2) 対数振幅スペクトル
freq = np.fft.rfftfreq(len(x), 1/fs)
logmag = np.log(np.abs(np.fft.rfft(x)) + 1e-12)
axes[1].plot(freq, logmag, lw=0.9)
axes[1].set_xlabel('Frequency [Hz]'); axes[1].set_ylabel('log|X|')
axes[1].set_title('Log-magnitude spectrum (note the ripple from echo)')
axes[1].grid(alpha=0.3)
# (3) ケプストラム
axes[2].plot(quefrency[:half], c[:half], lw=0.9)
axes[2].axvline(D/fs*1000, color='r', ls='--', label=f'true delay = {D/fs*1000:.1f} ms')
axes[2].set_xlabel('Quefrency [ms]'); axes[2].set_ylabel('c[n]')
axes[2].set_title('Real cepstrum (peak at echo delay)')
axes[2].legend(); axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('cepstrum_echo.png', dpi=150, bbox_inches='tight')
plt.show()
3枚のグラフから、ケプストラムの威力が読み取れます。まず中段の対数振幅スペクトルには、広帯域の全体形状の上に 細かいさざ波（リップル） が等間隔に乗っているのが見えます。これがエコー由来の $2\alpha\cos(\omega D)$ 項です。人間の目にはこのリップルの周期を正確に読むのは困難ですが、下段のケプストラムを見ると、ケフレンシ 5 ms の位置（赤破線）に 鋭く強いピーク（このデモでは $c[D]\approx 0.22$）が立っています。さらに 10 ms（$2D$）付近には 符号が反転した副成分（$c[2D]\approx -0.16$ の谷）が現れ、これは導出で予言した $-\tfrac{\alpha^2}{2}\delta[n-2D]$ の項（負号）に対応します。もとの音は強い周期性を持たないため、ケフレンシ軸にはエコー由来の構造だけが浮かび上がり、スペクトル上で読み取り困難だったエコー遅延が、ケプストラムでは一目瞭然のピークとして現れたわけです。
ピークを自動検出するコードも示します。実用では低ケフレンシ側（包絡成分）を除外してからピークを探すのが定石です。
import numpy as np
# 低ケフレンシ（包絡）を避けて 1ms 以降でピークを探す
c_half = c[:len(c)//2]
qf_ms = np.arange(len(c_half)) / fs * 1000
search = qf_ms > 1.0 # 1ms より大きいケフレンシのみ対象
peak_idx = np.argmax(c_half * search)
print(f"検出されたエコー遅延: {qf_ms[peak_idx]:.2f} ms（真値 {D/fs*1000:.2f} ms）")
このコードは 1 ms 以下の低ケフレンシ（スペクトル包絡が集中する領域）をマスクしてから最大値を探します。出力では検出遅延が真値とぴたり一致する 5.00 ms となり、ケプストラムによる遅延推定が正確に機能していることが確認できます。
エコー検出の手応えがつかめたところで、次は応用の花形であるピッチ検出に進みます。
Pythonでの実装2: 母音波形のピッチ検出
人間の母音は、声帯のパルス的な振動（励振源）が声道のフィルタを通って生まれます。ここでは、より現実に近い母音波形を合成し、そのケプストラムからピッチ（基本周波数）を抽出します。声道フィルタにはフォルマントを模した共鳴を入れ、励振源にはパルス列を使います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
fs = 16000 # サンプリング周波数 [Hz]
dur = 0.04 # フレーム長 [s]
t = np.arange(int(fs*dur)) / fs
f0 = 125.0 # 基本周波数 [Hz]（声帯振動）→ 周期 8ms
P = int(fs / f0) # 基本周期 [サンプル]
# 励振源: 基本周期 P ごとのインパルス列（声帯パルスの単純化）
excit = np.zeros(len(t))
excit[::P] = 1.0
# 声道フィルタ: 2つのフォルマント（共鳴）を持つIIRフィルタ
def formant_filter(x, formants, bws, fs):
"""フォルマント周波数とバンド幅から共鳴フィルタを構成して適用"""
y = x.copy()
for fc, bw in zip(formants, bws):
r = np.exp(-np.pi * bw / fs) # 極の半径
theta = 2*np.pi*fc/fs # 極の角度
a = [1, -2*r*np.cos(theta), r**2] # 2次共鳴の分母
y = sig.lfilter([1], a, y)
return y
# 母音「あ」風: F1=730Hz, F2=1090Hz
vowel = formant_filter(excit, [730, 1090], [80, 90], fs)
vowel *= np.hanning(len(vowel))
print(f"基本周波数 {f0} Hz, 基本周期 {P} サンプル = {1000/f0:.2f} ms")
このコードは、基本周期 P ごとのインパルス列（励振源）を、2つのフォルマント（F1=730 Hz、F2=1090 Hz）を持つ共鳴フィルタに通して母音「あ」風の波形を合成しています。formant_filter は各フォルマントを 2 次の共鳴 IIR フィルタとして実装しています。出力から基本周期が 128 サンプル（8 ms）であることが分かります。次にこの母音のケプストラムを計算します。
import numpy as np
import matplotlib.pyplot as plt
# real_cepstrum と vowel, fs, P, f0 を使う
c = real_cepstrum(vowel)
qf_ms = np.arange(len(c)) / fs * 1000
half = len(c)//2
fig, axes = plt.subplots(3, 1, figsize=(10, 9))
# (1) 母音波形
axes[0].plot(np.arange(len(vowel))/fs*1000, vowel, lw=0.8)
axes[0].set_xlabel('Time [ms]'); axes[0].set_ylabel('Amplitude')
axes[0].set_title('Synthetic vowel waveform')
axes[0].grid(alpha=0.3)
# (2) 対数振幅スペクトル（ハーモニクス構造）
freq = np.fft.rfftfreq(len(vowel), 1/fs)
logmag = np.log(np.abs(np.fft.rfft(vowel)) + 1e-12)
axes[1].plot(freq, logmag, lw=0.8)
axes[1].set_xlim(0, 4000)
axes[1].set_xlabel('Frequency [Hz]'); axes[1].set_ylabel('log|X|')
axes[1].set_title('Log spectrum: slow envelope (formants) + fast ripple (harmonics)')
axes[1].grid(alpha=0.3)
# (3) ケプストラム
axes[2].plot(qf_ms[:half], c[:half], lw=0.9)
axes[2].axvline(1000/f0, color='r', ls='--', label=f'pitch period = {1000/f0:.2f} ms')
axes[2].set_xlim(0, 15)
axes[2].set_xlabel('Quefrency [ms]'); axes[2].set_ylabel('c[n]')
axes[2].set_title('Cepstrum: peak at pitch period')
axes[2].legend(); axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('cepstrum_pitch.png', dpi=150, bbox_inches='tight')
plt.show()
中段の対数スペクトルには、フォルマントによる ゆるやかな山（包絡） の上に、基本周波数間隔で並ぶ くし状の細かいピーク（ハーモニクス） が乗っているのがはっきり見えます。これが「ゆっくり変化する声道成分 + 速く振動する励振成分」の足し算です。下段のケプストラムでは、低ケフレンシ（0〜2 ms 付近）に包絡由来の大きな成分が集中し、その先のケフレンシ 8 ms（赤破線、基本周期）に 明瞭なピッチピーク が立っています。包絡と励振がケフレンシ軸上で見事に分離されているのが読み取れます。
ピッチを自動推定する関数を作り、人間の声の範囲（おおむね 60〜400 Hz）に対応するケフレンシ帯でピークを探します。
import numpy as np
def cepstral_pitch(x, fs, f_min=60, f_max=400):
"""ケプストラムのピーク位置から基本周波数 f0 を推定する。"""
c = real_cepstrum(x)
# f_max → 小さいケフレンシ、f_min → 大きいケフレンシ
q_min = int(fs / f_max) # 探索する最小ケフレンシ index
q_max = int(fs / f_min) # 探索する最大ケフレンシ index
segment = c[q_min:q_max]
peak = np.argmax(segment) + q_min # ピークのケフレンシ index
f0_est = fs / peak
return f0_est, peak
f0_est, peak = cepstral_pitch(vowel, fs)
print(f"推定 f0 = {f0_est:.2f} Hz（真値 {f0:.2f} Hz）, ピーク {peak} サンプル")
この関数は、探したい基本周波数の範囲を逆数を通してケフレンシ index の範囲に変換し（高い周波数ほど小さいケフレンシに対応）、その区間で最大ピークを探します。出力では推定 $f_0$ が真値 125 Hz にきわめて近い値となり、ケプストラム法が頑健なピッチ検出器であることが確認できます。ノイズや倍音の振幅変動に対して、ピーク位置はほとんど動かないのがこの手法の強みです。
ピッチ検出ができたので、次はリフタリングを使って、母音から声道スペクトル包絡（フォルマント）を分離して取り出してみます。
Pythonでの実装3: リフタリングによる包絡と微細構造の分離
最後に、ケプストラムの真骨頂であるリフタリングを実装します。低ケフレンシだけを残すローパスリフターでスペクトル包絡（声道・フォルマント）を、高ケフレンシだけを残すハイパスリフターで励振源（ハーモニクス）を取り出します。
import numpy as np
def liftering(x, fs, lifter_cutoff_ms=4.0, mode='lowpass'):
"""ケプストラムのリフタリングで対数スペクトルを分離する。
mode='lowpass' : 低ケフレンシのみ残す → スペクトル包絡
mode='highpass' : 高ケフレンシのみ残す → 励振（微細構造）
"""
n_fft = len(x)
spectrum = np.fft.fft(x, n_fft)
log_mag = np.log(np.abs(spectrum) + 1e-12)
c = np.fft.ifft(log_mag).real
# リフター窓（実ケプストラムは対称なので両端を残す）
cutoff = int(fs * lifter_cutoff_ms / 1000)
lifter = np.zeros(n_fft)
if mode == 'lowpass':
lifter[:cutoff] = 1
lifter[-cutoff:] = 1 # 対称成分も残す
else: # highpass
lifter[:] = 1
lifter[:cutoff] = 0
lifter[-cutoff:] = 0
c_lift = c * lifter
# 再びフーリエ変換して対数スペクトルへ戻す
log_mag_lift = np.fft.fft(c_lift).real
return log_mag_lift
この関数は、ケプストラムにリフター窓を掛けてから再度フーリエ変換し、リフタリング後の対数スペクトルを返します。mode='lowpass' ではカットオフ未満のケフレンシ（と対称な高index側）だけを残し、mode='highpass' ではその逆を残します。実ケプストラムの対称性を考慮し、両端を扱っている点に注意してください。これを母音に適用して可視化します。
import numpy as np
import matplotlib.pyplot as plt
# vowel, fs を使う
freq = np.fft.rfftfreq(len(vowel), 1/fs)
n_half = len(freq)
log_full = np.log(np.abs(np.fft.rfft(vowel)) + 1e-12)
env = liftering(vowel, fs, lifter_cutoff_ms=4.0, mode='lowpass')[:n_half]
fine = liftering(vowel, fs, lifter_cutoff_ms=4.0, mode='highpass')[:n_half]
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# (1) 元の対数スペクトル と 包絡
axes[0].plot(freq, log_full, lw=0.8, alpha=0.6, label='log|X| (full)')
axes[0].plot(freq, env, 'r-', lw=2, label='envelope (low-time lifter)')
axes[0].set_xlim(0, 4000)
axes[0].set_xlabel('Frequency [Hz]'); axes[0].set_ylabel('log magnitude')
axes[0].set_title('Spectral envelope by low-time liftering (formants)')
axes[0].legend(); axes[0].grid(alpha=0.3)
# (2) 微細構造（励振）
axes[1].plot(freq, fine, 'g-', lw=0.9, label='fine structure (high-time lifter)')
axes[1].set_xlim(0, 4000)
axes[1].set_xlabel('Frequency [Hz]'); axes[1].set_ylabel('log magnitude')
axes[1].set_title('Excitation fine structure (harmonics)')
axes[1].legend(); axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('cepstrum_liftering.png', dpi=150, bbox_inches='tight')
plt.show()
上段では、元のギザギザした対数スペクトル（薄い線）の上に、ローパスリフターで抽出した 滑らかな包絡（赤線） が重なっています。この包絡は、もとの設定（F1=730 Hz, F2=1090 Hz）に対応して、おおよそ 750 Hz と 1075 Hz に2つの山（フォルマント F1、F2）を持っており、声道の共鳴特性を正しく取り出せていることが分かります。ここでカットオフを lifter_cutoff_ms=4.0 に取っているのが効いています。これより小さく（たとえば 2 ms に）すると2つのフォルマントが約 800 Hz の単一の山に融合してしまい、F1 と F2 が分離できません。フォルマント間隔に対応するケフレンシまで包絡側に含めるよう、カットオフを十分に確保する必要があるのです。下段のハイパスリフター出力は、包絡を取り除いた ハーモニクス構造（微細構造） だけを残しており、基本周波数間隔で振動する成分が見えます。一つの母音波形から、声道（包絡）と励振源（微細構造）という、物理的に別々の起源を持つ2成分を完全に分離できたわけです。これがホモモルフィック・デコンボリューションの基本であり、ボコーダーや音声符号化、音声合成の土台になっています。
最後に、雑音に対する頑健性を確認するため、ノイズを加えた母音でピッチ検出の安定性を見ておきます。
import numpy as np
# vowel, fs, f0 を使う。様々なSNRでピッチ検出を試す
rng = np.random.default_rng(0)
print("SNR [dB] | 推定 f0 [Hz]")
for snr_db in [40, 20, 10, 5, 0]:
sig_power = np.mean(vowel**2)
noise_power = sig_power / (10**(snr_db/10))
noisy = vowel + np.sqrt(noise_power) * rng.standard_normal(len(vowel))
f0_est, _ = cepstral_pitch(noisy, fs)
print(f" {snr_db:5d} | {f0_est:6.2f}")
出力を見ると、SNR が 40 dB から 0 dB まで下がっても推定 $f_0$ は真値 125 Hz の近傍（おおむね 125〜127 Hz）に踏みとどまる、という頑健な挙動が観察できます。これは、ケプストラムが「周期構造の周期」というロバストな特徴を捉えており、各ハーモニクスの振幅が雑音で乱れても、それらが作るピーク位置自体は容易には動かないためです。さらに SNR を下げていけば、いずれは雑音由来の偽ピークがピッチピークを上回って推定が破綻しますが、その耐性はかなり高いことが分かります。時間領域の自己相関法と並んで、ケプストラム法が実用的なピッチ検出器として長く使われてきた理由がここにあります。
まとめ
本記事では、ケプストラム分析について、その発想・数学的定義・導出・実装までを一通り解説しました。
- 発想: ケプストラムは「対数スペクトルのフーリエ変換」、すなわちスペクトルのスペクトルです。畳み込み $x=e*h$ をフーリエ変換で積 $X=EH$ に、対数で和 $\log|X|=\log|E|+\log|H|$ に変え、もう一度フーリエ変換することで、もとは分かちがたかった2成分をケフレンシ軸上の別々の場所へ分離します（ホモモルフィック処理）。
- 定義とケフレンシ: 実ケプストラム $c[n]=\mathrm{IDFT}\{\log|X[k]|\}$ の独立変数 $n$ は時間の次元を持つケフレンシで、ケフレンシ位置 $\tau$ は対数スペクトルが $1/\tau$ [Hz] 周期で振動する成分に対応します。
- ピーク位置の意味: エコー $x=s*(\delta+\alpha\delta[n-D])$ では $\log|H|$ に周期 $2\pi/D$ のリップルが生じ、ケプストラムは $n=D,2D,\dots$ にピークを作ります。母音では基本周期 $P$ の位置にピッチピークが立ち、$f_0=f_s/P$ で基本周波数が得られます。
- 実 vs 複素: 検出にはピーク位置だけ分かればよい実ケプストラムが手軽で安定。元の波形を再構成したい分離（ホモモルフィック・デコンボリューション）には位相を保持する複素ケプストラムが必須です。
- リフタリング: ケフレンシ領域のフィルタリング。低ケフレンシ通過で声道スペクトル包絡（フォルマント）、高ケフレンシ通過で励振源（ハーモニクス）を分離できます。
- 実装: 合成エコー信号でケプストラムが遅延 5 ms をピークとして検出し、母音波形では 125 Hz のピッチを雑音下でも頑健に推定し、リフタリングで包絡と微細構造を分離できることを確認しました。
ケプストラムは、メル周波数ケプストラム係数（MFCC）として音声認識・話者認識の標準特徴量へと発展し、さらに地震波探査やレーダーのエコー解析など、分野を超えて使われています。次のステップとして、以下の記事も参考にしてください。