STFTの完全再構成とCOLA制約の理論と実装

音声のノイズ除去アプリを思い浮かべてください。録音された音を細かい時間フレームに切り分け、それぞれをフーリエ変換して周波数成分を眺め、不要な成分を削り、最後に時間波形へ戻して再生します。この「切り分けて → 周波数で加工して → つなぎ戻す」という流れの中核にあるのが短時間フーリエ変換(STFT)です。ところが、ここには見落とされがちな落とし穴があります。何も加工せずにSTFTして逆STFTしただけなのに、戻ってきた音が「ヴォン、ヴォン」と周期的に音量が揺れてしまうことがあるのです。

この音量のうねりは、信号処理の世界で振幅変調(amplitude modulation)と呼ばれる現象です。原因は加工処理ではなく、フレームを切り出すときに掛ける窓関数の重なり方にあります。窓を一定の間隔でずらしながら足し合わせたとき、その総和が場所によって濃くなったり薄くなったりすると、それがそのまま音量のムラとして残ります。逆に、窓の総和がどこでも一定になるように設計すれば、加工しない限り元の信号がそっくりそのまま戻ってきます。この「窓のシフト和が一定」という条件こそが、本記事の主役である COLA(Constant-OverLap-Add)制約 です。

COLAを理解することは、単なる理論的興味にとどまりません。具体的な応用先を2つ挙げると、第一に音声強調・ノイズ抑圧があります。スペクトルサブトラクションやウィーナーフィルタは周波数領域でゲインを掛けますが、COLAが成立していないと加工以前の段階で波形が歪み、加工結果の評価ができません。第二にボコーダーやタイムストレッチ(音の高さを変えずに再生速度を変える処理)があります。位相ボコーダーは解析窓と合成窓を分けて使いますが、両者の積がCOLAを満たすことが完全再構成の前提になります。本記事では、このCOLA条件をオーバーラップ加算の式から省略なく導出し、条件が破れたときに何が起こるかを数式とPythonの両面から確かめていきます。

STFTの基本操作: フレーム分割と窓掛け

上の図は、STFT処理の3段階を示しています。元信号(上段)を一定間隔で切り出してフレームに分割し(中段)、各フレームにハニング窓を掛けて両端をゼロへ落とすと(下段)、各フレームが滑らかな断片になります。この「フレーム分割 → 窓掛け → DFT」というサイクルがSTFTの基本操作であり、逆方向の「IDFT → フレームを元の位置に戻す → 重ねて足す」がオーバーラップ加算です。COLAは、この最後の「重ねて足す」段階で振幅ムラが生じないための条件を規定します。

本記事の内容

  • STFTと逆STFT(オーバーラップ加算)の復習と直感的理解
  • COLA制約 $\sum_m w(n – mH) = c$ の導出(省略なし)
  • COLAが破れると振幅変調が残ることの数式的説明
  • ハニング窓の50%/75%オーバーラップ、NOLA条件、合成窓の選び方
  • Pythonで各窓・各ホップ長のCOLA和を計算し、再構成SNRを評価

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

STFTと逆変換をおさらいする

COLA制約はSTFTの逆変換の構造から自然に現れます。そこでまず、STFTがどんな操作で、どうやって元に戻すのかを直感的に押さえておきましょう。

長い音声信号 $x[n]$ をそのまま一括でフーリエ変換すると、「いつ・どの周波数成分が含まれていたか」という時間情報が消えてしまいます。会話の冒頭で高い声、後半で低い声があったとしても、全体のスペクトルには両方が混ざって現れるだけです。そこで、信号を短い時間フレームに区切り、各フレームごとにフーリエ変換します。これがSTFTの発想です。

フレームを区切るとき、ただブツ切りにすると切れ目で不連続が生じ、スペクトルに偽の高周波(スペクトル漏れ)が出ます。これを抑えるために、フレームの両端をなめらかに0へ落とす窓関数 $w[n]$ を掛けます。さらに、窓で両端を削った分の情報損失を補うため、フレームを少しずつ重ねながら(オーバーラップさせながら)スライドさせます。

代表的な窓関数の形状比較

上の図は、信号処理でよく使われる4種類の窓関数の形状を比較しています。矩形窓は両端まで値が1で急峻に打ち切られ、スペクトル漏れが大きいという弱点を持ちます。ハニング窓・ハミング窓・ブラックマン窓は、中央付近が高く両端に向かって滑らかにゼロへ落ちる釣鐘型で、漏れを大幅に低減できます。COLAの観点では、この「両端がゼロに落ちる性質」がオーバーラップしたときの窓の補い合いを可能にする鍵になります。

STFTの定義

長さ $N$ の窓関数 $w[n]$ ($n = 0, 1, \dots, N-1$)、ホップ長(フレームをずらす間隔) $H$ サンプルとしたとき、第 $m$ フレームのSTFTは次のように定義されます。

$$ X_m[k] = \sum_{n=0}^{N-1} w[n] \, x[n + mH] \, e^{-j 2\pi k n / N} $$

ここで $k = 0, 1, \dots, N-1$ は周波数ビンの番号、$m = 0, 1, 2, \dots$ はフレーム番号です。中身を分解すると、$x[n + mH]$ は信号を $mH$ サンプルだけ進めて読み出す部分、$w[n]$ はそれに窓を掛ける部分、残りの指数項は通常の離散フーリエ変換(DFT)です。つまりSTFTは「窓を掛けて切り出した断片をDFTする」操作を、フレームを $H$ ずつずらしながら繰り返したものに他なりません。

逆変換はオーバーラップ加算で行う

では、各フレームのスペクトル $X_m[k]$ から元の信号 $x[n]$ をどう復元するのでしょうか。まず各フレームを逆DFT(IDFT)して時間波形に戻します。

$$ y_m[n] = \frac{1}{N} \sum_{k=0}^{N-1} X_m[k] \, e^{+j 2\pi k n / N}, \quad n = 0, 1, \dots, N-1 $$

このとき、加工していなければ $y_m[n] = w[n] \, x[n + mH]$ が成り立ちます。逆DFTはDFTを完全に元に戻すので、フレームの中身である「窓を掛けた信号断片」がそのまま戻ってくるからです。

各フレーム $y_m[n]$ は元の時間軸上では位置 $mH$ から始まる長さ $N$ の波形です。これらは互いに重なっています。重なった部分を単純に足し合わせる操作を オーバーラップ加算(Overlap-Add, OLA) と呼びます。元の時間インデックス $n$ における再構成信号 $\hat{x}[n]$ は

$$ \hat{x}[n] = \sum_{m} y_m[n – mH] $$

です。ここで $y_m[n – mH]$ は、第 $m$ フレームを正しい時間位置 $mH$ に戻して評価することを意味します(フレーム内のローカル座標 $\ell = n – mH$ が $0 \le \ell \le N-1$ のときだけ寄与します)。

この「足し合わせる」操作のとき、窓の重なり方がそのまま再構成のクオリティを左右します。次のセクションで、加工なしの理想的な場合に $\hat{x}[n]$ が $x[n]$ にどこまで一致するかを式の上で追いかけ、COLA条件を導きます。

オーバーラップ加算の仕組みとCOLAエンベロープ

上の図は、ハニング窓(50%重複)の各フレームがどのように重なり合い、最終的な窓の総和(OLAエンベロープ)を作るかを示しています。5枚の窓(上段)が50%ずつ重なって足し合わさると、下段の赤い曲線(OLAエンベロープ)は破線の定数1(c=1)に正確に重なります。この「総和が均一」という状態こそがCOLA成立であり、端の過渡区間を除けばどの位置でも同じ重みで元信号を回収できることを意味します。

COLA制約を導出する

ここからが本記事の核心です。「窓のシフト和が一定」という条件がなぜ完全再構成と等価なのかを、式を一行ずつ追って導出します。ゴールは、再構成信号 $\hat{x}[n]$ が元信号 $x[n]$ の定数倍になるための窓とホップ長の条件を求めることです。

加工なしのときの再構成信号

加工を一切しない場合、前節で見たように $y_m[n] = w[n] \, x[n + mH]$ です。これをオーバーラップ加算の式に代入します。再構成信号は

$$ \hat{x}[n] = \sum_{m} y_m[n – mH] $$

でした。フレーム $m$ のローカル座標で $y_m[\ell] = w[\ell] \, x[\ell + mH]$ なので、$\ell = n – mH$ を代入すると

$$ y_m[n – mH] = w[n – mH] \, x\big[(n – mH) + mH\big] = w[n – mH] \, x[n] $$

ここで重要なのは、$x$ の引数が $(n – mH) + mH = n$ となって $m$ に依存しなくなる点です。つまり、どのフレームから来た寄与も、時間位置 $n$ では同じ値 $x[n]$ を運んでいます。違うのは掛かっている窓の値 $w[n – mH]$ だけです。これをオーバーラップ加算の総和に戻すと

$$ \hat{x}[n] = \sum_{m} w[n – mH] \, x[n] = x[n] \sum_{m} w[n – mH] $$

となります。$x[n]$ は $m$ に依存しないので和の外へ括り出せました。

COLA条件の登場

上の結果を見つめましょう。再構成信号は

$$ \hat{x}[n] = x[n] \cdot \underbrace{\sum_{m} w[n – mH]}_{\displaystyle W(n)} $$

という形をしています。$W(n) \equiv \sum_{m} w[n – mH]$ を 窓のシフト和 あるいは OLAエンベロープ と呼びます。これは「位置 $n$ に重なっている全フレームの窓の値を足したもの」です。

ここで決定的な観察ができます。もし $W(n)$ が $n$によらず一定値 $c$ になるならば

$$ \hat{x}[n] = c \cdot x[n] $$

となり、再構成信号は元信号のちょうど $c$ 倍です。$c$ で割り戻せば完全に一致します。逆に、$W(n)$ が $n$ によって変動すれば、その変動パターンが $x[n]$ に掛け算され、波形が時間とともに伸び縮みします。これがまさに振幅変調です。以上から次の条件が得られます。

$$ \begin{equation} \sum_{m=-\infty}^{\infty} w[n – mH] = c \quad (\text{定数}, \ \forall n) \end{equation} $$

この条件を COLA(Constant-OverLap-Add)制約 と呼びます。窓 $w$ とホップ長 $H$ がこの式を満たすとき、ペア $(w, H)$ は「COLAを満たす」と言います。完全再構成は加工しない限りこの一本の条件に集約されます。

ここで一度立ち止まって、この結果のありがたさを噛みしめておきましょう。STFTは信号を細切れにして窓で削り、しかも互いに重ねながら足し戻すという、一見すると情報が失われそうな操作の連続です。それにもかかわらず、窓のシフト和という極めてシンプルな量さえ一定に保てば、削った情報も重ねた冗長性もすべて辻褄が合って元通りになる、というのがCOLAの主張です。窓の形そのもの(釣鐘型か台形か)や、各フレームでどんな周波数が出たかは一切関係ありません。問われているのは「重なりを足したときに濃淡ができないか」という一点だけです。この単純さこそが、STFTが音声・レーダー・通信の現場で安心して使われている理由でもあります。

周期性から見たCOLA

シフト和 $W(n)$ には隠れた構造があります。窓を $H$ ずつずらして足したものなので、$W(n)$ は周期 $H$ の周期関数です。実際

$$ W(n + H) = \sum_m w[(n + H) – mH] = \sum_m w[n – (m-1)H] = \sum_{m’} w[n – m’H] = W(n) $$

と、和のインデックスを $m’ = m – 1$ と置き換えるだけで $W(n + H) = W(n)$ が示せます。したがって $W(n)$ が定数かどうかは、1周期分 $n = 0, 1, \dots, H-1$ だけ調べれば十分です。実装でCOLA和を確認するときは、この $H$ サンプル分の値がすべて等しいかを見ればよい、ということになります。

このCOLAの周期性は、後で周波数領域からCOLAを言い換えるときにも効いてきます。次は、COLAが破れたときに具体的にどんな歪みが出るのかを掘り下げましょう。

COLA成立と非成立のOLAエンベロープ比較

上の図は、4通りの窓・ホップ長の組み合わせでOLAエンベロープを比較したものです。緑で示したハニング窓50%と75%はエンベロープが完全に水平(定数)となり、COLA成立を視覚的に確認できます。一方、赤で示したハニング窓30%(H=0.7N)は同じ窓でもホップ長が中途半端なためエンベロープが波打ち、ハミング窓50%も僅かに変動しています。$H$ を変えるだけでCOLA成否が逆転するという、ホップ長選定の重要性が明瞭に表れています。

COLAが破れると何が起こるか

COLAが「完全再構成の条件」だと分かりました。では条件を満たさないと、どれくらい・どんなふうに信号が壊れるのでしょうか。直感を養うために、破れたときの歪みを式で表してみます。

再構成は窓エンベロープによる振幅変調になる

加工なしの再構成信号は $\hat{x}[n] = W(n) \, x[n]$ でした。$W(n)$ が定数でない場合、これは「元信号に時間変動するゲイン $W(n)$ を掛けた」信号です。$W(n)$ は周期 $H$ の周期関数なので、フーリエ級数に展開できます。正規化角周波数を $\omega_H = 2\pi / H$ とすると

$$ W(n) = c_0 + \sum_{p=1}^{\infty} \big[ a_p \cos(p \omega_H n) + b_p \sin(p \omega_H n) \big] $$

と書けます。$c_0$ は平均値(DC成分)、第2項以降が変動成分です。これを $\hat{x}[n] = W(n) x[n]$ に代入すると

$$ \hat{x}[n] = c_0 \, x[n] + \sum_{p=1}^{\infty} \big[ a_p \cos(p \omega_H n) + b_p \sin(p \omega_H n) \big] x[n] $$

第1項は欲しかった信号(定数倍)です。問題は第2項です。三角関数と $x[n]$ の積は変調そのものですから、$x[n]$ のスペクトルが $\pm p \omega_H$ だけ周波数シフトしたコピー(側帯波, sideband)を生みます。たとえば $x[n]$ が純音 $\cos(\omega_0 n)$ なら、再構成後には $\omega_0 \pm p\omega_H$ という余計な周波数が出現します。

聴感上は、$\omega_H = 2\pi/H$ がフレームレートに対応する低い周波数なので、「フレームの繰り返し周期での音量のうねり」として知覚されます。これがアプリの音が「ヴォン、ヴォン」と揺れる正体です。COLAを満たせば変動成分 $a_p, b_p$ がすべて0になり、側帯波が消え、$\hat{x}[n] = c_0 x[n]$ という純粋な定数倍だけが残ります。

側帯波という言葉を使ったのには理由があります。これはアナログ通信の振幅変調(AM放送)とまったく同じ数学だからです。AMでは搬送波に音声を掛けて側帯波を作りますが、ここでは原信号 $x[n]$ に窓エンベロープの変動 $W(n) – c_0$ を掛けてしまっている、いわば「意図しない自己変調」が起きています。違いは、AMでは変調が目的なのに対し、STFTの再構成では変調は完全に不要な副作用だという点です。そして変調の深さ(うねりの大きさ)は、窓エンベロープのリプル、すなわち $\sum_{p\geq 1}\sqrt{a_p^2 + b_p^2}$ の大きさで決まります。リプルが小さいほどうねりも小さく、ゼロになればCOLA完全再構成、というわけです。後のPython実験では、まさにこのリプルの大きさと再構成誤差が対応することを数値で確かめます。

COLAが破れると振幅変調が残ることの可視化

上の図は、同じ200Hz純音を素朴なOLAで再構成したときの違いを示しています。上段が元信号、中段はハニング窓50%(COLA成立)での再構成で波形の振幅が均一に保たれています。下段のハニング窓30%(COLA非成立)では、明らかに振幅が周期的に揺れる「うねり」が再構成信号に現れており、窓エンベロープの変動がそのまま振幅変調として波形に写り込んでいます。加工を一切しなくても窓とホップ長の選択ミスだけで音質が劣化するという、COLA成立の実用的な重要性を直感的に示しています。

OLAエンベロープのスペクトル: COLA成立ならDC成分のみ

上の図は、OLAエンベロープを1周期分フーリエ解析したスペクトルです。緑(COLA成立・ハニング50%)は直流(0 cycles/sample)の成分だけがあり、それ以外の周波数は文字通りゼロです。これはエンベロープに変動成分が一切ないことを意味します。一方、赤(COLA非成立・ハニング30%)は直流以外の複数の周波数にも成分があり、これらがそれぞれ $\pm p \omega_H$ の側帯波として元信号のスペクトルに畳み込まれます。「COLA成立 = スペクトルにDCのみ = 変調なし」という対応が、周波数領域から見たCOLAの意味です。

NOLA条件 — 加工する場合の最低ライン

COLAは「加工しないなら完全再構成」を保証する条件でした。しかし実際のSTFTアプリは周波数領域で何かしら加工します。加工して逆変換するとき、別の弱い条件が必要になります。それが NOLA(Nonzero OverLap-Add)条件 です。

$$ \begin{equation} \sum_{m=-\infty}^{\infty} w^2[n – mH] \neq 0 \quad (\forall n) \end{equation} $$

直感的には「どの時間位置 $n$ でも、窓の二乗和が0にならない」という条件です。なぜ二乗和かというと、加工後の標準的な逆STFTでは、解析窓と同じ窓を合成側でも掛けて(これを重み付きオーバーラップ加算, Weighted OLAと呼びます)、最後に窓二乗和エンベロープ $\sum_m w^2[n – mH]$ で割り戻すアルゴリズムを使うからです。割り戻すために、この二乗和エンベロープがどこでも0であってはいけません。もしある位置で二乗和が0だと、その時刻の信号情報がすべての窓に消されてしまい、復元不可能になります。

COLAとNOLAの関係を整理すると、COLAは「無加工なら完全再構成できる」強い条件、NOLAは「加工しても逆STFTがそもそも定義できる(0除算が起きない)」最低限の条件です。窓が非負でどこかで0にならなければNOLAは自動的に満たされます。実際のシステム設計では、解析窓と合成窓の積に対してCOLAを満たすように両者を選ぶのが定石です。次のセクションで、その合成窓の選び方を具体的に見ていきます。

解析窓・合成窓とCOLAの実用設計

これまでは「窓を掛けて切り出し → 逆DFT → そのまま足す」という素朴な逆変換(OLA)を考えてきました。実用の逆STFTでは、合成側にもう一枚窓を掛けるWeighted OLAが標準です。なぜ二枚使うのか、どう選ぶのかを整理します。

なぜ合成窓を掛けるのか

加工なしのOLAでは、フレームの中心は重みが大きく端は重みが小さい状態でそのまま足されます。しかし周波数領域で加工(ゲイン操作)を加えると、逆DFTで戻したフレーム $y_m[n]$ はもはやきれいに窓の掛かった断片ではなくなり、フレーム端で不連続が生じます。この端の不連続をなめらかに減衰させるため、合成側で再び窓 $s[n]$ を掛けてから足します。重み付きオーバーラップ加算の再構成式は

$$ \hat{x}[n] = \frac{\displaystyle \sum_m s[n – mH] \, y_m[n – mH]}{\displaystyle \sum_m s[n – mH] \, w[n – mH]} $$

と書けます。分子は合成窓を掛けて足したもの、分母は解析窓と合成窓の積のシフト和です。加工なしのとき $y_m[n-mH] = w[n-mH] x[n]$ を代入すると、分子は $x[n] \sum_m s[n-mH] w[n-mH]$ となり、分母とちょうど打ち消し合って $\hat{x}[n] = x[n]$ になります。つまり完全再構成の条件は 積 $s[n] w[n]$ がCOLAを満たすこと、すなわち

$$ \sum_m s[n – mH] \, w[n – mH] = c \quad (\forall n) $$

に一般化されます。$s[n] = 1$(合成窓なし)とすれば、これは元の $\sum_m w[n-mH] = c$ に戻ります。解析窓も合成窓も同じ平方根窓 $\sqrt{w[n]}$ にすれば積が $w[n]$ になり、$w$ がCOLAを満たすよう選べばよい、という設計がよく使われます。

ハニング窓のCOLA成立条件

最も実用的な例が ハニング窓(Hann window) です。長さ $N$ のハニング窓は

$$ w[n] = \frac{1}{2}\left(1 – \cos\frac{2\pi n}{N}\right) = \sin^2\frac{\pi n}{N}, \quad n = 0, 1, \dots, N-1 $$

で定義されます(両端で0になめらかに落ちる釣鐘型)。ハニング窓は特定のオーバーラップ率でちょうどCOLAを満たすという美しい性質があります。代表的なのが 50%オーバーラップ($H = N/2$)です。

50%オーバーラップでハニング窓のシフト和が定数になることを確かめましょう。位置 $n$ には2枚の窓が重なります。第0フレームの $w[n]$ と、半フレーム前から始まる第 $-1$ フレームの $w[n + N/2]$ です(ローカル座標で $n + N/2$)。コサイン形を使うと、足し算は

$$ w[n] + w[n + N/2] = \frac{1}{2}\left(1 – \cos\frac{2\pi n}{N}\right) + \frac{1}{2}\left(1 – \cos\frac{2\pi (n + N/2)}{N}\right) $$

第2項のコサインの中身を展開すると $\frac{2\pi n}{N} + \pi$ なので、$\cos(\theta + \pi) = -\cos\theta$ を使って

$$ = \frac{1}{2} – \frac{1}{2}\cos\frac{2\pi n}{N} + \frac{1}{2} – \frac{1}{2}\cos\left(\frac{2\pi n}{N} + \pi\right) = \frac{1}{2} – \frac{1}{2}\cos\frac{2\pi n}{N} + \frac{1}{2} + \frac{1}{2}\cos\frac{2\pi n}{N} $$

コサイン項が打ち消し合って

$$ w[n] + w[n + N/2] = 1 $$

となります。$n$ に依存せず定数1。見事にCOLAが成立しました。同様に 75%オーバーラップ($H = N/4$)では4枚の窓が重なり、シフト和は定数 $2$ になります(後でPythonで確認します)。一般にハニング窓は $H = N/2, N/4, N/8, \dots$ のように $N$ を割る適切な値でCOLAを満たします。

ハニング窓50%オーバーラップのCOLA成立の代数的証明

上の図は、COLA成立の代数的な仕組みを可視化しています。左パネルでは2枚の窓 $w[n]$(オレンジ)と $w[n+N/2]$(青)がどのように補い合って定数1(点線)に収束するかを示しており、重なりが深い中央付近では一方が高く他方が低いという「補完関係」が見て取れます。右パネルは、この補完の根拠となる数式 $\cos(\theta+\pi) = -\cos\theta$ を視覚化したものです。両フレームのコサイン成分(オレンジと青)が全域で符号が反転しており、足し合わせると完全にゼロとなる(黒破線)ことが分かります。このコサイン項の相殺こそが、ハニング窓でCOLAが厳密に成立する数学的根拠です。

矩形窓・その他の窓

矩形窓($w[n] = 1$)は、$H = N$(オーバーラップなし)なら当然シフト和が1で定数なのでCOLAを満たします。しかし矩形窓はスペクトル漏れが大きいので解析には不向きです。ハミング窓は両端が0でない(オフセット0.08が残る)ため、50%オーバーラップでは厳密には定数になりませんが、適切なホップ長を選べば近似的にCOLAを満たします。ブラックマン窓など他の窓も、ホップ長次第でCOLAを満たすかどうかが変わります。

ここまでで理論は出そろいました。COLAが定数和の条件であること、破れると振幅変調が出ること、合成窓を含めた一般化、ハニング窓の具体的なCOLA成立を導きました。あとは、これらが本当に数値で確かめられるのかをPythonで検証しましょう。

解析窓と合成窓の積がCOLAを満たすWOLA設計の可視化

上の図は、3通りのWOLA設計における解析窓・合成窓・その積の形状を示しています。設計Aの平方根ハニング窓は、解析と合成に同じ $\sqrt{w}$ を使うため積が $w$(ハニング)になり、50%オーバーラップでCOLAが保証されます。設計Bは解析がハニング・合成が矩形(合成窓なし相当)で、積がハニングそのものになり同様にCOLAを満たします。設計CのハミングはCOLAを厳密には満たさないため、用途に応じてどの設計を採用するかが実装の要点となります。

Pythonでの実装と検証

理論で「COLAを満たせば完全再構成、破れれば振幅変調」と分かりました。ここからは実際にコードを書いて、(1)各窓・各ホップ長のCOLA和がどうなるか、(2)再構成誤差がCOLA成立で本当に消えるか、をSNRで定量的に確かめます。

COLA和を可視化する

まず、いろいろな窓とホップ長についてシフト和 $W(n) = \sum_m w[n – mH]$ を計算し、定数になるかどうかを目で見てみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import windows

# ===== 窓のシフト和(COLAエンベロープ)を計算する関数 =====
def cola_envelope(window, hop, num_frames=12):
    """窓をhopずつずらして足したシフト和を返す。中央付近の1区間を抽出。"""
    N = len(window)
    total_len = hop * (num_frames - 1) + N
    env = np.zeros(total_len)
    for m in range(num_frames):
        start = m * hop
        env[start:start + N] += window  # オーバーラップ加算
    # 端のフレームが揃っていない区間を避け、中央部を切り出す
    center = total_len // 2
    return env, center

N = 256  # 窓長
# 検証する (窓, ホップ長, ラベル) の組
configs = [
    (windows.hann(N, sym=False), N // 2, "Hann, 50% (H=N/2)"),
    (windows.hann(N, sym=False), N // 4, "Hann, 75% (H=N/4)"),
    (windows.hann(N, sym=False), int(N * 0.7), "Hann, 30% (H=0.7N)"),
    (windows.hamming(N, sym=False), N // 2, "Hamming, 50% (H=N/2)"),
]

このコードは、窓を hop サンプルずつずらしながら配列に足し込むことでシフト和を組み立てています。num_frames を十分多くとり、端の立ち上がり・立ち下がりの影響を受けない中央部分だけを後で見るようにしています。次に、この関数を使って各設定のエンベロープを描きます。

fig, ax = plt.subplots(figsize=(11, 6))
for window, hop, label in configs:
    env, center = cola_envelope(window, hop, num_frames=16)
    # 中央付近の hop*4 サンプルを表示(定常部)
    seg = env[center - 2 * hop: center + 2 * hop]
    ax.plot(seg, linewidth=2, label=label)

ax.set_xlabel("Sample index (centered segment)")
ax.set_ylabel(r"OLA envelope  $W(n)=\sum_m w[n-mH]$")
ax.set_title("Constant-OverLap-Add (COLA) envelope for various windows/hops")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("cola_envelope.png", dpi=150, bbox_inches="tight")
plt.show()

このグラフから3つのことが読み取れます。第一に、ハニング窓の50%オーバーラップ(H=N/2)は中央部で平坦な水平線(値1)になっており、COLAが厳密に成立していることが分かります。第二に、ハニング窓の75%(H=N/4)も平坦ですが値が2になっています。これは4枚重なるためで、定数であることに変わりはなくCOLAを満たします。第三に、ハニング窓でもホップ長が中途半端な30%(H=0.7N)だと、エンベロープが波打っており定数になっていません。この波打ちこそが振幅変調を生む元凶です。ハミング窓50%もわずかに波打っており、厳密にはCOLAを満たさないことが見て取れます。

COLA和のリプル(変動量)を定量化する

目で見るだけでなく、エンベロープがどれだけ定数からずれているかを数値化しましょう。定常部での最大値・最小値・相対リプルを計算します。

import numpy as np
from scipy.signal import windows

def cola_ripple(window, hop, num_frames=24):
    """定常部でのエンベロープの平均・最小・最大・相対リプルを返す。"""
    N = len(window)
    total_len = hop * (num_frames - 1) + N
    env = np.zeros(total_len)
    for m in range(num_frames):
        env[m * hop: m * hop + N] += window
    # 全フレームが重なる定常区間 [N, total_len - N] を解析
    steady = env[N: total_len - N]
    mean = steady.mean()
    ripple = (steady.max() - steady.min()) / mean  # 相対リプル
    return mean, steady.min(), steady.max(), ripple

N = 256
tests = [
    (windows.hann(N, sym=False), N // 2, "Hann 50%"),
    (windows.hann(N, sym=False), N // 4, "Hann 75%"),
    (windows.hann(N, sym=False), int(N * 0.7), "Hann 30%"),
    (windows.hamming(N, sym=False), N // 2, "Hamming 50%"),
    (windows.boxcar(N), N, "Rect no-overlap"),
]
print(f"{'config':<18}{'mean':>8}{'min':>8}{'max':>8}{'ripple':>10}")
for window, hop, label in tests:
    mean, lo, hi, rip = cola_ripple(window, hop)
    print(f"{label:<18}{mean:8.4f}{lo:8.4f}{hi:8.4f}{rip:10.2e}")

このコードの出力では、相対リプル(ripple 列)がCOLA成立の判定材料になります。ハニング窓50%・75%、矩形窓オーバーラップなしはリプルが $10^{-15}$ 程度(浮動小数点の丸め誤差レベル)で、実質ゼロ=COLA成立です。一方、ハニング窓30%は数%〜十数%のリプルが出ており、ハミング窓50%も無視できないリプルが残ります。リプルが大きい設定ほど、後で見る再構成誤差が悪化することが予想できます。リプルの大きさと再構成品質が対応するはずだ、という仮説を次の実験で検証します。

NOLA条件: 窓二乗和エンベロープが全位置で0にならないことの確認

上の図は、NOLA条件の確認として各窓の二乗和エンベロープを可視化したものです。ハニング・ハミング・ブラックマンいずれも窓長の内側では二乗和が確実に正の値を保ち、どこでも0になる箇所がありません。矩形窓オーバーラップなしは、ちょうどフレームの切れ目で値が0になるように見えますが、信号長がちょうど窓長の整数倍であれば問題ありません。NOLA条件が破れる典型は、極端に狭いホップ長や特殊な窓形状で信号のある時刻に全フレームの窓がゼロになってしまう場合であり、実用的な窓選択では事実上起こりません。

STFT/逆STFTを実装して再構成SNRを測る

いよいよ本丸です。解析→合成のSTFTパイプラインを自前で実装し、加工なしで再構成したときの誤差をSNRで評価します。まずSTFTと逆STFTの関数を定義します。

import numpy as np

def stft(x, window, hop):
    """信号xをSTFT。各列が1フレームのDFT。"""
    N = len(window)
    num_frames = 1 + (len(x) - N) // hop
    X = np.empty((N, num_frames), dtype=complex)
    for m in range(num_frames):
        seg = x[m * hop: m * hop + N] * window  # 窓掛け
        X[:, m] = np.fft.fft(seg)               # フレームのDFT
    return X

def istft_ola(X, window, hop, length):
    """逆STFT(重み付きオーバーラップ加算 + COLA正規化)。"""
    N = len(window)
    num_frames = X.shape[1]
    y = np.zeros(length)
    win_sumsq = np.zeros(length)  # 窓二乗和エンベロープ(正規化用)
    for m in range(num_frames):
        frame = np.real(np.fft.ifft(X[:, m]))   # 逆DFT
        start = m * hop
        y[start:start + N] += frame * window     # 合成窓を掛けて加算
        win_sumsq[start:start + N] += window ** 2 # 二乗和を蓄積
    # NOLA: 二乗和が0でない位置だけ割り戻す
    nz = win_sumsq > 1e-12
    y[nz] /= win_sumsq[nz]
    return y

ここでは解析窓と合成窓に同じ window を使い、合成後に窓二乗和エンベロープ $\sum_m w^2[n-mH]$ で割り戻すWeighted OLA方式を実装しています。win_sumsq > 1e-12 の判定がNOLA条件のチェックに対応し、二乗和が0の位置では割り算を回避しています。この正規化により、$w^2$ がCOLAを満たさなくても加工なしなら理論上は完全再構成できる作りになっています。

# ===== テスト信号で再構成SNRを測る =====
np.random.seed(0)
fs = 16000
T = 1.0
t = np.arange(int(fs * T)) / fs
# 複数の正弦波 + 軽いノイズ(音声を模した広帯域信号)
x = (np.sin(2*np.pi*220*t) + 0.6*np.sin(2*np.pi*440*t)
     + 0.4*np.sin(2*np.pi*1300*t) + 0.05*np.random.randn(len(t)))

def reconstruction_snr(x, window, hop):
    X = stft(x, window, hop)
    y = istft_ola(X, window, hop, length=len(x))
    N = len(window)
    # 端のフレーム境界効果を避け、内側だけで誤差評価
    s = slice(N, len(x) - N)
    err = x[s] - y[s]
    snr = 10 * np.log10(np.sum(x[s]**2) / np.sum(err**2))
    return snr

このコードは、複数の正弦波とノイズを混ぜた音声風の信号を作り、STFT→逆STFTを通したあと、元信号との差からSNR(信号対誤差比)を計算します。端のフレームは窓が全部重なっておらずCOLAが崩れるので、内側のサンプルだけで評価しています。SNRが高い(数十dB以上)ほど再構成が正確だと判断できます。

窓・ホップ長を変えてSNRを比較する

最後に、いろいろな窓とホップ長の組み合わせで再構成SNRを一覧にし、COLA成立と再構成品質の対応を確かめます。

import numpy as np
from scipy.signal import windows

N = 1024
cases = [
    (windows.hann(N, sym=False),    N // 2, "Hann 50%"),
    (windows.hann(N, sym=False),    N // 4, "Hann 75%"),
    (windows.hann(N, sym=False),    int(N*0.7), "Hann 30%"),
    (windows.hamming(N, sym=False), N // 2, "Hamming 50%"),
    (windows.hamming(N, sym=False), N // 4, "Hamming 75%"),
    (windows.blackman(N, sym=False),N // 4, "Blackman 75%"),
    (windows.boxcar(N),             N,      "Rect no-ov"),
]
print(f"{'config':<16}{'recon SNR [dB]':>16}")
for window, hop, label in cases:
    snr = reconstruction_snr(x, window, hop)
    print(f"{label:<16}{snr:16.1f}")

この出力が本記事の結論を数値で裏づけます。注目すべきは、窓二乗和エンベロープで正規化するWeighted OLA方式を使っているため、加工しない限りNOLAさえ満たせばどの設定でも高いSNR(おおむね 100 dB を超える丸め誤差レベル)が出るという点です。これは「加工なしなら、合成側の正規化が窓のムラを吸収してくれる」ことを意味します。一方で、もし正規化を外した素朴なOLA(合成窓も正規化もなし)にすると、COLAを満たさない設定だけSNRが極端に悪化します。この対比を最後に可視化しましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import windows

def istft_naive(X, window, hop, length):
    """正規化なしの素朴なOLA(合成窓なし)。COLA依存性を見るため。"""
    N = len(window); y = np.zeros(length)
    for m in range(X.shape[1]):
        frame = np.real(np.fft.ifft(X[:, m]))
        y[m*hop: m*hop + N] += frame  # ただ足すだけ
    return y

def snr_naive(x, window, hop):
    X = stft(x, window, hop)
    y = istft_naive(X, window, hop, len(x))
    N = len(window); s = slice(N, len(x) - N)
    # COLA定数cで割り戻してから比較(理想は c で一定倍)
    c = np.sum([window[(np.arange(N) ) % N] for _ in range(1)])  # placeholder
    return x[s], y[s]

N = 1024
labels, ripples, snrs = [], [], []
for window, hop, label in cases:
    # COLAリプル
    env = np.zeros(hop*30 + N)
    for m in range(30): env[m*hop:m*hop+N] += window
    steady = env[N:hop*30]; rip = (steady.max()-steady.min())/steady.mean()
    # 素朴OLAのSNR(定数倍を最小二乗で補正してから誤差評価)
    X = stft(x, window, hop); y = istft_naive(X, window, hop, len(x))
    s = slice(N, len(x)-N)
    alpha = np.dot(x[s], y[s]) / np.dot(x[s], x[s])  # 最適スケール
    err = y[s] - alpha * x[s]
    snr = 10*np.log10(np.sum((alpha*x[s])**2)/np.sum(err**2))
    labels.append(label); ripples.append(rip); snrs.append(snr)

fig, ax1 = plt.subplots(figsize=(11, 6))
xpos = np.arange(len(labels))
ax1.bar(xpos - 0.2, snrs, width=0.4, color="steelblue", label="Naive-OLA recon SNR [dB]")
ax1.set_ylabel("Reconstruction SNR [dB]", color="steelblue")
ax1.set_xticks(xpos); ax1.set_xticklabels(labels, rotation=30, ha="right")
ax2 = ax1.twinx()
ax2.bar(xpos + 0.2, ripples, width=0.4, color="indianred", label="COLA ripple")
ax2.set_ylabel("COLA ripple (max-min)/mean", color="indianred")
ax1.set_title("COLA ripple vs naive-OLA reconstruction SNR")
plt.tight_layout()
plt.savefig("cola_snr.png", dpi=150, bbox_inches="tight")
plt.show()

このグラフから、COLAリプル(赤)と素朴OLAの再構成SNR(青)がきれいに逆相関することが読み取れます。ハニング窓50%・75%や矩形窓オーバーラップなしはリプルがほぼ0で、定数倍を補正したあとのSNRが非常に高く(丸め誤差レベル)、完全再構成が実現しています。逆にハニング窓30%やハミング・ブラックマン窓の一部の設定はリプルが大きく、たとえ最適な定数倍で補正してもSNRが頭打ちになります。これは残ったムラが単なる定数倍では消せない振幅変調(側帯波)であることを示しています。理論で予測した「COLA成立 ⇔ 完全再構成、COLA破れ ⇔ 振幅変調」が、数値実験でそのまま裏づけられました。

COLAリプルと再構成SNRの対応バーチャート

上のバーチャートから、COLAリプル(赤)が実質ゼロの設定(ハニング50%、ハニング75%、矩形オーバーラップなし)では再構成SNR(青)が300 dB超の丸め誤差レベルに達し、完全再構成が数値的に確認できます。一方、ハニング30%はCOLAリプルが0.815と大きく、SNRが10 dB台にとどまります。ハミング窓各設定のリプルはほぼゼロで高SNRが出ていますが、これはハミング窓50%が実は厳密にはCOLAを満たすためです(オフセット0.08が両フレームで打ち消し合い定数1.08になる)。この結果は、「リプルゼロ = 完全再構成」という理論の予測を定量的に裏づけています。

まとめ

本記事では、STFTの解析→合成で原信号を完全に取り戻すための条件であるCOLA制約を、オーバーラップ加算の構造から導出し、Pythonで検証しました。

  • 逆STFTの正体: 各フレームを逆DFTして時間位置に戻し、重なりを足し合わせるオーバーラップ加算。加工なしなら各フレームの寄与は同じ $x[n]$ を運び、違いは掛かった窓の値だけ
  • COLA制約: 再構成信号は $\hat{x}[n] = x[n] \sum_m w[n-mH]$ となる。窓のシフト和 $W(n) = \sum_m w[n-mH]$ が定数 $c$ なら完全再構成($\hat{x} = c\,x$)。これが $\sum_m w[n-mH] = c$ というCOLA条件
  • COLAの周期性: $W(n)$ は周期 $H$ なので、1周期分($H$ サンプル)が定数かを見れば判定できる
  • 破れたときの歪み: $W(n)$ の変動成分が $x[n]$ に掛かり、フレームレート $2\pi/H$ の側帯波を生む振幅変調が残る
  • NOLA条件: 加工する場合の最低限の条件 $\sum_m w^2[n-mH] \neq 0$。窓二乗和で割り戻す逆STFTが0除算を起こさないための条件
  • 合成窓とハニング窓: 解析窓と合成窓の積 $s[n]w[n]$ がCOLAを満たせばよい。ハニング窓は50%(和=1)・75%(和=2)オーバーラップでCOLAを厳密に満たす
  • 数値検証: COLAリプルと再構成SNRは逆相関し、COLA成立時のみ丸め誤差レベルの完全再構成が得られることを確認

COLAは、音声強調・タイムストレッチ・ボコーダーといったSTFTベースの処理すべての土台です。ここで身につけた「窓のシフト和を一定にする」という視点を持っておくと、実装で謎の音量うねりに出会ったとき、真っ先に窓とホップ長のCOLA成立を疑えるようになります。

次のステップとして、以下の記事も参考にしてください。

画像なし
短時間フーリエ変換(STFT)の理論と実装 — 時間-周波数解析の基礎
STFTの定義・スペクトログラム・不確定性原理まで、時間-周波数解析の核となる概念をPythonで実装しながら解説します。COLA条件を理解する前提知識として必読です。
スペクトル漏れと窓関数の選び方(コヒーレント利得・ENBW)の理論と実装
スペクトル漏れの原因とハニング・ハミング・ブラックマン窓の特性(ENBW・コヒーレント利得)を比較し、目的に合った窓の選択指針を解説します。COLA制約と合わせて読むと窓関数の理解が深まります。
画像なし
オーバーラップ加算法・保留法による高速畳み込みの理論と実装
FFTを使った高速畳み込み(オーバーラップ加算法・保留法)の理論とPython実装を解説します。COLAの基礎となるオーバーラップ加算の信号処理的な応用例として読むと理解が深まります。