2.1 Waveファイルを読んで離散Fourier変換して周波数成分を調べる

授業で、Mathematica を用いて、WAVE ファイル ``guitar-5-3.wav'' を読み込み、適当な時点からの1秒分のデータ (項数 $N$ をサンプリング周波数と同じにとった) を離散Fourier変換して、 $\vert C_0\vert$ , $\vert C_1\vert$ , $\dots$ , $\vert C_{N-1}\vert$ をプロットする、という実験を行った ([1] の第4章)。

Python を使ってそれと同じことをしてみよう。

次のプログラムでは、 グラフ描画のために matplotlib というモデュールを用いている ¹。

http://nalab.mind.meiji.ac.jp/~mk/fourier/testwavemodule3.py

   1 # testwavemodule3.py
   2 # encoding: utf-8
   3 # http://nalab.mind.meiji.ac.jp/~mk/lecture/fourier-2017/testwavemodule3.py
   4 # これは Python 3.x 用のプログラム
   5 
   6 import wave
   7 import numpy as np
   8 import matplotlib.pyplot as plt
   9 
  10 fname = 'guitar-5-3.wav'
  11 wfile = wave.open(fname, 'r')
  12 
  13 numch = wfile.getnchannels()
  14 samplewidth = wfile.getsampwidth()
  15 samplerate = wfile.getframerate()
  16 numsamples = wfile.getnframes()
  17 
  18 print("チャンネル数 = ", numch)
  19 print("サンプル幅 (バイト数) = ", samplewidth)
  20 print("サンプリングレート(Hz) =", samplerate)
  21 print("サンプル数 =", numsamples)
  22 print("録音時間 =", numsamples / samplerate)
  23 
  24 # すべてのフレームを読み込む (bytesオブジェクトになる)
  25 buf = wfile.readframes(numsamples)
  26 wfile.close()
  27 
  28 # numpy の ndarray に変換する
  29 if samplewidth == 2:
  30     data = np.frombuffer(buf, dtype='int16')
  31     data = data / 32768.0
  32 elif samplewidth == 4:
  33     data = np.frombuffer(buf, dtype='int32')
  34     data = data / 2147483648.0
  35 
  36 # ステレオの場合は左チャンネルだけを取り出す
  37 # (0 から最後まで2つおきに、つまり 0,2,4,.. 番目を取り出す)
  38 if numch == 2:
  39     #l_channel = data[0::2]
  40     #r_channel = data[1::2]
  41     data = data[0::2]
  42 
  43 # 62000番目から1秒分 (samplerate 個) 取り出し、離散フーリエ変換する
  44 start = 62000
  45 N = samplerate
  46 c = np.fft.fft(data[start:start+N])
  47 c = abs(c)
  48 
  49 plt.subplot(2,1,1)
  50 plt.title('data')
  51 plt.plot(range(start, start+N), data[start:start+N])
  52 
  53 plt.subplot(2,1,2)
  54 plt.title('frequency spectrum')
  55 # 対数目盛りにするには次の3行をアンコメントする(注釈を外す)
  56 #ax = plt.gca()
  57 #ax.set_yscale('log')
  58 #ax.set_xscale('log')
  59 
  60 plt.plot(c, linestyle='-')
  61 #横軸を周波数(Hz)にするには、次の2行をアンコメントする
  62 #freqList = np.fft.fftfreq(N, d=1.0/samplerate)
  63 #plt.plot(freqList, c, linestyle='-')
  64 
  65 plt.tight_layout()  # タイトルの被りを防ぐ
  66 plt.show()

57行目で単に

plt.plot(c, linestyle='-')

としているので、授業でやった Mathematica を用いた実験と同様に、 $\left\vert C_0\right\vert$ , $\left\vert C_1\right\vert$ , $\dots$ , $\left\vert C_{N-1}\right\vert$ が順に プロットされるが、横軸の値の説明が少し面倒である ($N$ をサンプリング周波数に等しく取っているので、 $1$ 番目は「直流成分」、 $n$ 番目 ( $2\le n\le N/2$ ) は周波数 $n-1$ (Hz) の成分ということになる)。

図 1: 下のグラフは $\vert C_0\vert$ , $\vert C_1\vert$ , $\dots$ , $\vert C_{N-1}\vert$ をプロット したもの

Numpy で用意されている numpy.fft.fftfreq() を使って

freqList = np.fft.fftfreq(N, d=1.0/samplerate)
plt.plot(freqList, c, linestyle='-')

とすると、横軸の絶対値が周波数(Hz)を表すようになる (詳しいことは省略するので、必要があれば Numpy のドキュメントを読むこと)。

図 2: numpy.fft.fftfreq() を利用して横軸を周波数 (ただし符 号付き) に直す

同様に、1つめのグラフ (PCMデータのプロット) の横軸も、 単位が時間になるように直すべきかもしれない (横軸、縦軸の説明は必須で、即答するためには、 説明がしやすいものにしておくのが望ましい)。