PythonでISTFT

今回は、前回のSTFTの逆の処理を行うISTFTについて実験します。

ISTFTについて
1. ISTFTとは
2. PythonでISTFT
試してみる
バンドパスフィルタを作る
まとめ

ISTFTについて

ISTFTとは

ISTFTとは、Inverse of Short-Time Fourier Transform＝逆短時間フーリエ変換の略で、要するにSTFTの逆の処理を行うものです。つまり、時刻ごと・周波数ごとの振幅・位相の情報から信号を復元します。難しい理論はともかく、プログラムを作るだけならその程度の理解で十分です。

三角関数の重ね合わせ、と思っておけばいいのかな。

Pythonでは、scipy.signal モジュールの istft() メソッドで簡単にISTFTを行うことができます。
（stft()メソッドと同様、istft() メソッドは scipyのver.0.19.0以降で使うことができます。これより古いものがインストールされている場合はバージョンアップして下さい）

istft() メソッドの基本的な使い方は以下の通りです。

t, data = scipy.signal.istft(S, fs, window, nperseg)

引数

S
時間ごと・周波数ごとの成分の振幅・位相の情報を与えます。2次元の配列（numpy.ndarray）で、各要素は振幅と位相を表す複素数の値です。scipy.signal.stft() メソッドの返却値と同じ、

[ [S(t1,f1), S(t1,f2), …, S(t1, fn) ],
[S(t2,f1), S(t2,f2), …, S(t2, fn) ],
　…
[S(tn,f1), S(tn,f2), …, S(tn, fn) ] ]

という型式です。

fs
サンプリング周波数です。浮動小数点数で指定します。省略可能で、デフォルト値は1.0です。

window
窓関数の形状です。通常は省略で良いでしょう。

nsperseg
セグメントの大きさを表します。省略された場合はサンプリング周波数とSの周波数軸のデータ数から算出されます。

返却値

t
周波数のリスト（numpy.ndarray）。

data
変換結果の波形データのリスト（numpy.array）。

試してみる

それでは、実際に scipy.signal.istft() を使って音声信号を分析してみましょう。

サンプルプログラム

実験に使うプログラムは以下の通りです。

内容はごく簡単で、stft() した結果を istft() に渡せば、元の信号に戻るはず…というものです。まぁどうしても変換で情報量が落ちるので完全な復元にはならないのでしょうが…

import numpy
import soundfile
from scipy.signal import stft,istft

src, samplingfreq = soundfile.read("input.wav", dtype="int16")

f, t, spectrum = stft(src, fs=samplingfreq)
t, dst = istft(spectrum, fs=samplingfreq)

soundfile.write("output.wav", numpy.array(dst, dtype="int16"), samplingfreq, format="WAV", subtype="PCM_16")

コード解説

STFTの時と同じ部分については省略。

フーリエ変換と逆変換

f, t, spectrum = stft(src, fs=samplingfreq)
t, dst = istft(spectrum, fs=samplingfreq)

stft() メソッドでフーリエ変換した結果をすぐ istft() メソッドに代入しています。これで（誤差は含むものの）元の信号が復元されます。

stft() の返却値のうち f と t は、samplingfreq と spectrum の要素数から一意に計算できるので（STFTの記事参照）、istft() の引数として与える必要がありません。

ファイル出力

ファイル出力には、PySoundFileモジュールの write() メソッドを使用しています。

使い方は以下の通りです。

soundfile.write(path, data, samplingfreq, format=format, subtype=subtype)

引数

path
保存するファイルのパスを指定します。

data
音声信号のサンプリングデータ（モノラルなら1次元、ステレオなら2次元のリスト系データ）を指定します。個々の要素の値はint16・int32・float32・float64のいずれかでなければなりません。

samplingfreq
サンプリング周波数を指定します。

format
ファイルの形式を指定します。wavファイルの場合は”WAV”を指定します。他に”AIFF”、”AU”、”FLAC”、”OGG”などが指定できるようです。

subtype
各データの型を指定します。WAVファイルの場合はリニアPCMの”PCM_16″や”PCM_U8″、差分PCMの”MS_ADPCM”などが指定できます。

soundfile.write("output.wav", numpy.array(dst, dtype="int16"), samplingfreq, format="WAV", subtype="PCM_16")

ここでは、ファイル形式にWAV、データは符号付き16ビットリニアPCMを指定しています。

データは istft() の結果をそのまま渡しても一応動くのですがなぜかノイズが酷く、いったん明示的にint16型を要素とする numpy.ndarray に変換してから渡すようにしたら改善しました。 istft() の返却値は float32型を要素とする numpy.ndarray なのですが、内部の型変換でなにか変なことをやってるのかな…？

実験

まず、このファイルを入力してみます。

出力はこうなりました。

もっと音が歪んだりノイズが乗ったりするかと思いきや、区別がつきません。

別の例でも試してみましょう。
今度は複雑な高調波を含む人の声です。変換→逆変換によって高調波の比率が破壊されてしまうと言葉が聞き取れなくなるかもしれません。

出力はこうなりました。

何の問題もなく言葉が聞き取れます。PCスピーカで聞く限り音質が劣化したような感じもしません。

バンドパスフィルタを作る

どうやら stft() した結果を istft() すれば、ほぼもとの音声信号が戻ることが判りました、では次に、stft() の結果を使って信号を加工する『フィルタ』を作ってみます。

最初はごく簡単に、特定の周波数成分をカットするフィルタを作ります。これは、stft() の結果の中から該当する周波数成分の振幅を0（複素数なので0+0j）にしてから istft() すればできそうです。

サンプルプログラム

import numpy
import soundfile
from scipy.signal import stft,istft

cutoff_h = 3400
cutoff_l = 300

src, samplingfreq = soundfile.read("input.wav", dtype="int16")

f, t, spectrum = stft(src, fs=samplingfreq, nperseg=256)
spectrum[f<cutoff_l] = 0.0+0.0j
spectrum[f>cutoff_h] = 0.0+0.0j
t, dst = istft(spectrum, fs=samplingfreq)
soundfile.write("output.wav", numpy.array(dst, dtype="int16"), samplingfreq, format="WAV", subtype="PCM_16")

コード解説

フィルタ処理

f, t, spectrum = stft(src, fs=samplingfreq, nperseg=256)
spectrum[f<cutoff_l] = 0.0+0.0j
spectrum[f>cutoff_h] = 0.0+0.0j
t, dst = istft(spectrum, fs=samplingfreq)

ここのキモとなる記述は、

spectrum[f<cutoff_l] = 0.0+0.0j

でしょう。Python独特の簡潔な表現ですが、CやJava風に書くならば、

for i in range(len(f)):
    if f[i]<cutoff_l:
        spectrum[i]=0.0+0.0j

のような意味になります。

もう少し詳しく見ていくと、まず

f < cutoff_l

の結果は

[ True True True False False False … False ]

のような、True と False が並んだ配列になります。
（ f は値の小さい順に要素が並んでいるので、『○より小さい』という条件を与えると最初のいくつかが True、残りが False になります）

そして、

spectrum[ [True True True False False False … False] ]

は、spectrumの要素のうち、添字に与えたTrue／Falseの配列の、Trueと同じ位置にある要素を取り出す、という意味です。

実際にはspectrumは「『周波数ごとに要素の並んだ配列』が時間ごとに並ぶ2次元配列」ですが、それぞれの時間ごとに同じ周波数に対応する要素すべてに 0.0+0.0j が代入されます。

よって、『cutoff_lよりも低い周波数成分に対応するデータがすべて0.0+0.0j＝存在しなくなる』ということです。

同様に

spectrum[f>cutoff_h] = 0.0+0.0j

では、『cutoff_hよりも高い周波数成分に対応するデータが全て0.0+0.0j＝存在しなくなる』ということです。

つまり、この2行の処理を行うと、信号に含まれていた周波数成分のうち cutoff_l～cutoff_h の範囲のものだけが残ることになります。

このように、特定の範囲の周波数成分だけを残す働きを持つフィルタのことを、バンドパスフィルタといいます。

※ちなみに、残す周波数の上限だけが設定され、それより低い周波数成分が残るフィルタをローパスフィルタ、残す周波数の下限だけが設定され、それより高い周波数成分が残るフィルタをハイパスフィルタ、特定の範囲の周波数を消し、それより高い周波数成分と低い周波数成分を残すフィルタをバンドカットフィルタ、といいます。