Pythonでフーリエ変換(と逆変換) - 音楽プログラミングの超入門（仮）

音響信号のフーリエ変換

ここでは、離散フーリエ変換のみを扱っています。
信号処理においてフーリエ変換というと、ほとんどの場合、離散フーリエ変換を指しています。

式

フーリエ変換は可逆変換なので、逆変換があります。フーリエ変換・逆変換の定義式はこんな感じです。

$\begin{eqnarray} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ H(f) &=& \int_{-\infty}^{\infty} h(t)e^{-i2\pi ft}\;\;dt, \\ h(t) &=& \int_{-\infty}^{\infty} H(f)e^{i2\pi ft}\;\;df \end{eqnarray}$

ここでは、 $h(t)$ が入力音響信号ですね。

つまり何をしているのか

上の式をぱっと見ても、何やこれよく分からん、という人が多いんじゃないでしょうか。
音響信号のフーリエ変換では、その音の中に何ヘルツの正弦波がどのくらいの大きさ(振幅)で入っているかと、その位相が分かります。もっと簡単に言うと、どの音の高さがどのくらい強く鳴っているかが分かります。
それさえ分かっていれば、深い話に突っ込んでいかない限り、大丈夫でしょう、多分。

何ヘルツまで解析できるのか

信号がサンプリング周波数 $f_s$ で標本化されているとき、その半分未満の周波数成分しか完全に再現することはできません。これを標本化定理といい、 $f_s/2$ をナイキスト周波数といいます。120[Hz]の音が含まれる信号を、200[Hz]でサンプリングしてもその情報は失われてしまいます。

逆に、ある音響信号に含まれる最大周波数の２倍より大きい周波数でサンプリングすれば、その音響信号は完全に再現できることになります。市販CDのサンプリング周波数が44.1[kHz]であるのは、人間の可聴域の上限が20[kHz]程度であるからだそうです。

窓関数

実は、フーリエ変換では、入力信号が周期的な信号のちょうど１周期分であるという前提を置いています。（だいぶはしょった・・。）なので、入力信号の始端と終端が綺麗に繋がるような形になっていることが望ましいですが、大抵そんなことは有りえません。

そのような信号をそのままフーリエ変換すると、元々は含まれていない周波数が現れたり、つまり解析結果に雑音が入ることになります。これを緩和するために、入力信号に、両端に近づくに従って0に近づくような関数をかけることで、擬似的な周期信号とします。この関数を窓関数といいます。音響信号処理で、フーリエ変換を用いるときは、必ずかけた方が良いでしょう。

f:id:yukara_13:20131110161014p:plain,left

上図は、入力信号、窓関数、窓関数をかけた入力信号を示しています。

実際にやってみましょう

# -*- coding: utf-8 -*-
from scipy import arange, hamming, sin, pi
from scipy.fftpack import fft, ifft
from matplotlib import pylab as pl

fs = 8000 # Sampling rate
L = 1024 # Signal length

# 440[Hz]のサイン波を作る。
sine_440 = sin(2. * pi * arange(L) * 440. / fs)
# 600[Hz]のサイン波を作る。
sine_600 = 2 * sin(2. * pi * arange(L) * 600. / fs)
# 800[Hz]のサイン波を作る。
sine_800 = 3 * sin(2. * pi * arange(L) * 800. / fs)

# 全部足す
sig = sine_440 + sine_600 + sine_800

# 窓関数
win = hamming(L)

# フーリエ変換
spectrum_nw = fft(sig) # 窓関数なし
spectrum = fft(sig * win) # 窓関数あり
half_spectrum_nw = abs(spectrum_nw[: L / 2 + 1])
half_spectrum = abs(spectrum[: L / 2 + 1])

# フーリエ逆変換
resyn_sig = ifft(spectrum)
resyn_sig /= win

# 図を表示
fig = pl.figure()
fig.add_subplot(411)
pl.plot(sig)
pl.xlim([0, L])
pl.title("1. Input signal", fontsize = 20)
fig.add_subplot(412)
pl.plot(half_spectrum_nw)
pl.xlim([0, len(half_spectrum_nw)])
pl.title("2. Spectrum (no window)", fontsize = 20)
fig.add_subplot(413)
pl.plot(half_spectrum)
pl.xlim([0, len(half_spectrum)])
pl.title("3. Spectrum (with window)", fontsize = 20)
fig.add_subplot(414)
pl.plot(resyn_sig)
pl.xlim([0, L])
pl.title("4. Resynthesized signal", fontsize = 20)

pl.show()