はじめに
さて、自由研究の素材を「画像:MNIST」から「音:Sound」に変更すると宣言してから幾年月。
ようやく大詰めを迎えました。いよいよ大御所にしてラスボス「フーリエ変換」FFTの登場です。前回のlibrosaでも内部的には大活躍でしたが、直接対決です。
おおまかな内容
音を解析し、再構築するには「フーリエ変換」（もしくは高速フーリエ変換：FFT)を使うことが多いです。
離散フーリエ変換は
いきなりやってみました
元にする波形は準備（その５）でお世話になったプログラムで作成した９倍音まで含むノコギリ波です。
# coding:utf-8
import wave
import numpy as np
import scipy.fftpack
from scipy import signal
from pylab import *
if __name__ == "__main__":
# あらかじめ準備した ノコギリ波 C音（基底 261.626Hz:9倍音まで)を読み込み
wf = wave.open("Music/saw-sample-c.wav", "r")
fs = wf.getframerate() # サンプリング周波数 上のサンプルでは48K
x = wf.readframes(wf.getnframes())
x = frombuffer(x, dtype="int16") / 32768.0 # -1 - +1に正規化
wf.close()
start = 0 # サンプリングする開始位置
N = 2 ** 12 # FFTのサンプル数 2**12は4096です
# サンプリング数にあったハニング窓を準備します。
win = signal.hann(N)
# FFTを実行してスペクトラム情報を獲得します。
spectrum = scipy.fftpack.fft(x[start:start + N] * win)
freqList = scipy.fftpack.fftfreq(N, d=1.0 / fs)
# 振幅と位相をグラフ表示するために抽出します。
amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in spectrum] # 振幅スペクトル
phaseSpectrum = [np.arctan2(int(c.imag), int(c.real)) for c in spectrum] # 位相スペクトル
# スペクトラム情報から元の波形に復元します。
resyn_sig = ifft(spectrum)
# 1.元の波形を描画
figure(figsize=(8, 8))
subplot(511) # 5行1列のグラフの1番目の位置にプロット
subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.2, hspace=0.5)
plot(range(start, start + N), x[start:start + N])
axis([start, start + N, -1.0, 1.0])
xlabel("Fig 1. Original Wave / time [sample]")
ylabel("amp")
# 2.窓関数を通した波形を描画
subplot(512) # 5行1列のグラフの2番目の位置にプロット
plot(range(start, start + N), x[start:start + N] * win)
axis([start, start + N, -1.0, 1.0])
xlabel("Fig 2. through window Wave / time [sample]")
ylabel("amp")
# 3.振幅スペクトルを描画
subplot(513)
plot(freqList, amplitudeSpectrum, linestyle='-')
axis([0, 3000, 0, 600])
xlabel("Fig 3. frequency [Hz]")
ylabel("amp spect.")
# 4.位相スペクトルを描画
subplot(514)
plot(freqList, phaseSpectrum, linestyle='-')
axis([0, 3000, -np.pi, np.pi])
xlabel("Fig 4. frequency [Hz]")
ylabel("phase spect.")
# 5.元の波形を描画
subplot(515) # 5行1列のグラフの5番目の位置にプロット
plot(range(start, start + N), resyn_sig[start:start + N])
axis([start, start + N, -1.0, 1.0])
xlabel("Fig 5. resynth Wave / time [sample]")
ylabel("amp")
show()
ちなにみ、ちょっとだけscipyのFFTのほうが、numpyのFFTよりも早いという噂がありscipyのFFTを使っています。
図１.元にする波形の のこぎり波です
あらかじめ別のpythonプログラムで生成したノコギリ波（に近い）の波形です。
少し低い ド(C) で基底の周波数は261.626Hzです。９倍音までを重ねています。
音も貼ってみます。4096サンプルなのでとても短いです。
図２.窓関数を通した波形
そのまま使うと両端でクリッピングノイズが発生するので、窓関数で両端をゼロにしてやります。(正確にはサンプリングの「窓」（上記プログラムでは N ）の範囲で周期関数であると仮定して計算するため、終端がゼロでないと元の信号にない成分が導出されるのでこれを防ぐためためです）
これも音を貼ってみます。立ち上がり、立下りがクリップノイズなどがなく、スムーズに聞こえると思います。一層短く聞こえます。
図３.スペクトラムの振幅情報
フーリエ変換(FFT)して得たスペクトラムの振幅情報です。グラフで９倍音までちゃんと検出されていることがわかります。
図４.スペクトラムの位相情報
フーリエ変換(FFT)して得たスペクトラムの位相情報です。おもいのほか重要な情報です。
図５.スペクトラムから元に戻した波形
スペクトラム情報からFFT逆変換(scipy.ifft())を使い元の波形にしています。
ほぼ完ぺきなくらいに復元できてます。ここまで復元できるのは技術的には当たり前なのでしょうが、感覚的にはまるで魔法のようです。元の音と違いはわからないですが、いちおう貼っておきます。
ちなみに、窓関数をつかわないとこうなります。
フーリエ変換のまとめ
フーリエ変換は元になる波形情報からある範囲(窓)を計算し、以下を導出します
- 周波数
- 振幅
- 位相
窓関数を入れて計算することが多いです。
窓関数には用途に応じいろいろな種類がありますが、正弦波、または正弦波の組み合わせの解析にはハニング窓関数がよいとのことでした。
以下がハニング窓関数です。
...
みなかったことにしましょう。
ちなみにグラフはこちらです。
ハニング関数はscipy.signal.hanning()を使えます。scipy.signal.hann()も同じです。似た名前のハミング関数scipy.signal.hamming()は似ているけど別のものなので注意しましょう。（ハミングは両端がゼロでない関数です）
ハニング（ハン窓）関数は気象学者のJulius von Hannにより名づけられました。
ハミング関数はR. W. Hammingにより名づけられました。
低域の正確な検出には多くのサンプルが必要です
いろいろな周波数で確認していると問題点がでてきました。
高域ではまったく問題ない解像度なのですが、低域から中域にかけてサンプル数が少ないため分解能がひくいのです。
たとえば、サンプリング周波数が48kでサンプル数が1024の場合はfreqList=[0.0, 46.875, 93.75, ..]と約46.8Hz刻みです。これでは低域の半音検出もできません。
低音の「ド」は、C=32.703, C#=34.648 で、差は1.945Hz です。sin波の音をCかC#を特定するにはこの周波数付近では1Hzの解像度はほしいところです。
とするとサンプル数（＝窓の大きさ）Nは上限（＝サンプリング周波数）の48000までいってしまいました。これで低音の「ド」を解析させると、1Hz刻みで
amplitudeSpectrum[30] = 90.05 (30.0Hz)
amplitudeSpectrum[31] = 474.23 (31.0Hz)
amplitudeSpectrum[32] = 4315.14 (32.0Hz)
amplitudeSpectrum[33] = 5665.88 (33.0Hz)
amplitudeSpectrum[34] = 1734.12 (34.0Hz)
なので、原音C=32.703Hzは12半音識別できます。(ちなみにB=30.868, C#=34.648)
中域はA=440.000Hz, G#=415.305Hz, A#=466.164 で14Hz程度の識別は最低必要です。例えば「窓」サイズN = 4096 (2**12) で11.72Hz刻みになります。
高域A=3520.000Hz, G#=3322.438Hz, A#=3729.310Hz で209Hz程度の識別は最低必要です。例えば「窓」サイズN = 256 (2**8) で187.5Hz刻みになります。
低音の検出の「窓」は思いっきり広くとり、中音では程々に、高音では結構狭くてもかまわない感じです。（ただ、ビブラートとかチョーキングとか、フレットレスベース感を出したいので精度をどれくらいにするか、各種資源や処理速度とのせめぎあいです）どのように実装に反映するかは、実験しながら決めたいです。
特定の若者の能力の発現させるため20kHz以上の音波を使うSF
があったような、なかったような。ボクは、周波数の確認をしてるときに高い周波数（14000Hz以上）がまったく聞こえないことがわかり、ショックを受けました。あ、ボクはかえるなので人間のスペックと競っても意味ないか。。。
自分はどうかな？と興味のある方は
などで確認できます。音量にはくれぐれも注意してください。（聞こえないからといって決して音量をあげないでください。）
蛇足の中の余談ですが、かえるの耳はちゃんとした内耳構造をもっていますが、哺乳類などが持つ「かたつむり管（蝸牛）」はありません。
かえるの内耳には
- 低音を主に受容する sacculus（10-250 Hzの音を受容）
- 中音域を主に受容する amphibian papilla（100-1000Hzの周波数を受容）
- 高音を主に受容する basilar papilla（1-4 KHzの周波数を受容）
という三種類の音の受容器があります。音の振動により、リンパ液に浸った有毛細胞と蓋膜が、局所的にこすれあうことで音を検出します。
やはり、音の検出についてはボクの耳のように３種類くらいフーリエ変換を使い分けて行うのがよさそうです。
最初のサンプルと同じ曲で、グラフでも最初の７秒すこしはドラムがはいっていないことがわかります。曲と一緒に聴くと、今の音は”
