はじめに

出来る限り公式ドキュメントに従った実装を心掛けていますが、不慣れなところは見逃していただけると幸いです。

概観

短時間フーリエ変換

PythonによるFFTを用いたパワースペクトル推定では、時間領域信号をすべて周波数領域信号へ変換してその信号の性質を明らかにしようとする場合に用いられるパワースペクトルについて解説しました。しかし、音声・音楽など時間とともに変化する信号を解析する際には時間情報も重要な情報です。そのように時間情報、周波数情報の二つの情報をある一つの信号から抽出する際に用いられる解析を短時間フーリエ変換と呼びます。

パワースペクトルを求めることをスペクトル解析と呼ばれることが多いです。¹

窓関数

短時間フーリエ変換(以下、STFT)については別記事にて解説しますが、STFTでは窓関数という関数が重要です。今回はそんな窓関数に焦点を当てた記事を書きます。一般的な信号処理の教科書には窓関数についてページが割り当てられていますが、それらの窓関数はせいぜいblackman, hann,hammingウィンドウぐらいしか解説されていません、僕の学校の授業でもそんな感じでした。他にもたくさんあるはずですし、何なら特徴を数値的に表現して使い分けができるようになりたいです。この記事ではそんなことができるようになれたらなと思いながら書きました。

窓関数の主な利用方法は簡単で、ある信号に対して要素積をとるだけです。例えばある信号$x(t)$に対して窓関数$w(t)$をかけた(適用した)後の関数$y(t)$それらの関係式は以下のようになります。

y(t) = w(t)x(t)

ここで、ある信号$x(t)$をフーリエ変換したことを$X(\omega)=\mathcal{F}[x(t)]$と記述することにします。また、ある信号$X(\omega)$を逆フーリエ変換したことを$x(t)=\mathcal{F^{-1}}[X(\omega)]$と記述することにします。今から、

Y(\omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} W(\omega-\omega') X(\omega') d\omega' = \frac{1}{2\pi} \int_{-\infty}^{\infty} W(\omega') X(\omega-\omega') d\omega'

が成立した時それぞれ$x(t)$、$w(t)$、$y(t)$にはある信号$x(t)$に対して窓関数$w(t)$をかけた(適用した)後の関数$y(t)$であるという関係つまり、$y(t) = w(t)x(t)$が成立していることを示します。

\begin{align}
y(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} Y(\omega) e^{j\omega t}d\omega &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{1}{2\pi} \int_{-\infty}^{\infty} W(\omega-\omega') X(\omega') d\omega' e^{j\omega t}d\omega\\
&= \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{1}{2\pi} \int_{-\infty}^{\infty} W(\omega-\omega') X(\omega') d\omega' e^{j(\omega-\omega') t} e^{j\omega' t} d\omega\\
&= \frac{1}{4\pi^2} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} W(\omega-\omega') e^{j(\omega-\omega') t} X(\omega') e^{j\omega' t} d\omega' d\omega\\
&= \frac{1}{4\pi^2} \int_{-\infty}^{\infty} W(\omega-\omega') e^{j(\omega-\omega') t} d\omega \int_{-\infty}^{\infty}  X(\omega') e^{j\omega' t} d\omega'\\
&= \frac{1}{4\pi^2} \int_{-\infty}^{\infty} W(\beta) e^{j\beta t} d\beta \int_{-\infty}^{\infty}  X(\omega') e^{j\omega' t} d\omega'\\
&= \frac{1}{2\pi} \int_{-\infty}^{\infty} W(\beta) e^{j\beta t} d\beta \frac{1}{2\pi} \int_{-\infty}^{\infty}  X(\omega') e^{j\omega' t} d\omega'\\
&= w(t)x(t)
\end{align}

ちなみに以上の式変形から逆も成立します。上記の関係式から以下のことが言えます。

時間領域で窓関数をかけるということは、周波数領域で畳み込み積分を行うということである。

また窓関数を意識しない（$w(t) = 1$の場合、一般的に矩形窓と表現されます。）も周波数領域では$W(\omega)$の畳み込みになっています。なので$W(\omega)$の性質をしっかりと把握したうえで$w(t) = 1$という窓を用いるのか、もしくはその他の窓を用いるのかを考えるべきです。

窓の特徴について

Wikipediaに非常に良い画像があったので拝借しましたが、窓関数の性質は大きく以下の三つの要素によって変わります。

メインローブの大きさ(正確には違います。)
メインローブの太さ
サイドローブの大きさ(正確には違います。)

窓関数一覧

今回はscipy.signal.windowsに実装されている窓関数の中から、いくつかの例をあげてその窓関数のパワースペクトル(ただし、図において縦軸が$dB SPL$でないことに注意してください)を見てみることにします。

正確な話をすると、これ以降の図はすべてメインローブの大きさを基準とした(最大値を基準とした)デシベル値なのでサイドローブレベルが前述のメインローブの大きさに対応していると考えてください。

cosine : window with a simple cosine shape

blackman : Blackman window

hann : Hann window

hamming : Hamming window

blackmanharris : minimum 4-term Blackman-Harris window

nuttall : minimum 4-term Blackman-Harris window according to Nuttall

flattop : flat top window

考察

いかがでしょうか、見た感じ窓には二種類あるような気がします。

パターン1. サイドローブレベルは小さいが、中心周波数($0 kHz$)から離れるほどサイドローブが急激に小さくなる。
パターン2. サイドローブレベルは大きいが、中心周波数($0 kHz$)から離れてもサイドローブの大きさは変化しない。

cosine, blackman, hannウィンドウはパターン1に該当するような気がしませんか？そして、hamming, blackmanharris, nuttall, flattopウィンドウはパターン2に該当するような気がしませんか？この種類の分け方は僕が勝手に見た目で分けただけなので学術的には何かしら用語があるかもしれません(ないかもしれません笑)。とりあえず、この記事では窓の特徴を数値的に判断し、現行scipyに実装されている窓関数の使い分けができるようにすることを目標にしましょう。

ここで、どんな値を計算したら窓の特徴を分析したことになるでしょうか？この記事では以下の数値を計算することにします。

メインローブの幅(以下、メイン幅) [単位 : $Hz$]
メインローブのサイドローブの大きさの差(以下、サイドローブレベル) [単位 : $dB$]
サイドローブの平均値 [単位 : $dB$]

前提

図について

Googleで検索していただければわかりますが、基本的には$dB$に変換して表示してある図がほとんどです。たまに、フーリエ変換しただけの図を掲載しているサイトもあります。正直言って、どっちが良いのかはわかりません多分趣味だと思います。ここで、フーリエ変換しているだけの図を掲載する意味としては、窓関数は「周波数領域で畳み込み積分を行う」ということなので、窓関数はフィルタであるという解釈も出来ます。なので、フィルタとして考えるときは、そのまんまの複素振幅スペクトルの方が分かりやすいかもしれません。

しかし、前述の通り窓関数は短時間フーリエ変換で用いられることが多いです。その場合は大きさを$dB$に変換する場合がほとんどですし、周波数の分解能を考慮して窓を選択するはずです。ここで重要なのがメインローブ、サイドローブの大きさと幅です。それらがわかりやすいのはやはり本記事で採用しているようなメインローブを基準として$dB$変換を行った図($dB SPL$は最小可聴音を基準として$dB$変換を行ってます)だと考えます。本記事の図が$dB$に変換して表示しているのはそのためです。

図を生成する際の注意点

窓関数の信号点数は$2^7 = 128$点で生成しています。
サンプリング周波数は$44.1kHz$として計算しています。
上記の条件で計算すると、約$0.003$秒間の信号として得ることができます。

また、パワースペクトルの分解能(周波数ビンと呼ばれる場合もありますが、本記事ではパワースペクトルの分解能と表現します)は約$344Hz$です。信号点数を $2^{20} \simeq 1050000$ 点にした場合、約24秒の信号として得ることができ、周波数分解能は約$0.04Hz$です。

上記のパワースペクトルの分解能が低下する問題を解決するために、scipy.fft.fftにも記されている通りに、zeros padding²を行います。今回はパワースペクトルの分解能が$1Hz$になるように要素数を指定しました。

If it is larger, the input is padded with zeros. If n is not given, the length of the input along the axis specified by axis is used.

メインローブの大きさは、最も大きい値なのでnp.max(np.abs(c))によって算出しています。

ローブの大きさ

それぞれメインローブとサイドローブ(以下、ローブ)の大きさを算出する際はどのようにしたらよいでしょうか？図を見ればわかる通り、ローブの頂点は極大値であるはずです。極値は一般的に一階微分の正負が反転する場所に現れます(いつか数学的なテクニックをまとめてみたいですね)。窓関数の変化は周波数によって変化しています。なので、窓関数を$W(\omega)$つまり、周波数によって変化する一変数関数と考えて極値の分布を考えます。

以上を踏まえて、極値のインデックスを求めるための関数を以下ののように定義しました。また、cosineウィンドウの極値をプロットした場合を例にして表示てみます。

私は数学が得意ではないので、詳しい説明は数学が得意な人に聞いてみてください。

def get_extremum_idx(y:np.array) -> np.array :
    diffY = np.diff(y)
    diffY_sign = np.sign(diffY)
    diffY_sign_diff = np.diff(diffY_sign)
    diffY_sign_diff_abs = np.abs(diffY_sign_diff)
    return np.where( diffY_sign_diff_abs==np.max(diffY_sign_diff_abs) )[0] + 1

cosine : window with a simple cosine shape

パワースペクトルの対称性

$y(t)\in \mathbb{R}$をフーリエ変換した際の、結果を$\mathcal{F}[y(t)]=Y(\omega)$とします。ここで、$\mathcal{F}[y(t)]=Y(\omega)=\mathrm{Re}(\omega)-i\mathrm{Im}(\omega)$と表現できるとします。$y(t)\in \mathbb{R}$なので、$\omega>0$について

\begin{align}
\mathrm{Re}(\omega) &= \int_{-\infty}^{\infty} y(t) \cos(\omega t) dt \\
\mathrm{Im}(\omega) &= \int_{-\infty}^{\infty} y(t) \sin(\omega t) dt
\end{align}

となります。ここで、

\begin{align}
\mathrm{Re}(-\omega) &= \int_{-\infty}^{\infty} y(t) \cos(-\omega t) dt = \int_{-\infty}^{\infty} y(t) \cos(\omega t) dt = \mathrm{Re}(\omega) \\
\mathrm{Im}(-\omega) &= \int_{-\infty}^{\infty} y(t) \sin(-\omega t) dt = -\int_{-\infty}^{\infty} y(t) \sin(\omega t) dt = - \mathrm{Im}(\omega)
\end{align}

よって、

Y(-\omega)=\mathrm{Re}(-\omega)-i\mathrm{Im}(-\omega)=\mathrm{Re}(\omega)+i\mathrm{Im}(\omega)=\overline{Y(\omega)}

また、

\mathrm{ABSOLUTE}(Y(-\omega)) = \sqrt{ \mathrm{Re^2}(\omega) + \mathrm{Im^2}(\omega) } = \mathrm{ABSOLUTE}(Y(\omega))

ただし、

Y(0) = \int_{-\infty}^{\infty} y(t) e^{j0 t}dt = \int_{-\infty}^{\infty} y(t)dt

なので、$y(t)\in \mathbb{R}$をフーリエ変換した際、パワースペクトルは$0Hz$を除いて正負で対称であることがわかります。

パワースペクトルは0Hzを除いて正負で対称である。

特徴の計算

ローブの極値を求めることが出来たので、特徴的な値を求めます。それぞれの値の定義は下の図を参照してください。また、コーディングする際、対称性を考えて今後の値を求めます。

つまり、サイドローブの平均値を求める際は、下の図に示した黒い点の値を求めることになりますね。

cosine : window with a simple cosine shape

結果

窓関数	メイン幅[kHz]	サイドローブレベル[dB]	サイドローブの平均値[dB]
cosine	1.04	23.00	-72.43
blackman	2.09	96.22	-106.29
hann	1.39	31.47	-97.64
hamming	1.40	45.37	-54.20
blackmanharris	2.78	99.86	-113.30
nuttall	2.81	105.97	-100.10
flattop	3.52	108.27	-94.83

さいごに

表の結果から、blackmanharris、nuttall 、flattopが良い結果を出せそうな気がします。そこで、それぞれの特性を棒グラフで比較してみましょう。

メイン幅の比較～小さければ小さいほど良い～

サイドローブレベルの比較～大きければ大きいほど良い～

サイドローブの平均値の比較～小さければ小さいほど良い～

これらの結果から、基本的にはblackmanharrisを使えば問題はなさそうですが、サイドローブレベルにこだわりを持つ場合はnuttall 、flattopを利用することも検討するべきかもしれません。

番外編～紹介しきれなかった窓たち～

barthann ：a modified Bartlett-Hann window

bohman ：Bohman window

boxcar ：boxcar or rectangular window

parzen ：Parzen window

taylor ：Taylor window

triang ：triangular window

tukey ：Tukey window, also known as a tapered cosine window