概要
・Pythonで、xlsxを用いてExcelファイルを読み込んだ。・numpyを用いてFFT、pylabで結果を表示した。
・np.log2()の精度？に疑問が生じた。
ExcelファイルからFFTする
データ処理にPythonを使おうとしています。今は、4096個の電圧データをFFT(高速フーリエ変換、Fast Fourier Transform)しようとしています。Pythonには、numpy、scipy、pylabといった強力なライブラリが揃っているので、あまり自分で実装しなくとも、簡単にFFTなど演算が出来ます。楽でいいですね。今回は、"testdata.xlsx"というExcelファイルにデータが収まっています。このファイルのSheet1、1列目の1～4096行目に、小数でデータが入っています。この時、これらを取り出して、numpyでFFTを行い、pylabで表示するだけなら以下のスクリプトで出来ます。
# -*- coding: utf-8 -*-
__author__ = 'fenrir_'
import xlrd
import numpy as np
import math
from pylab import *
def rad_to_deg(rad):
return rad / np.pi * 180
if __name__ == "__main__":
# Excelからデータを読み出す。Sheet1の一番左の行。
filename = r"C:\Users\******\testdata.xlsx"
book = xlrd.open_workbook(filename)
sheet_1 = book.sheet_by_index(0)
# PythonのListから、numpy.ndarrayに変換
col_array_in_numpy = np.array(sheet_1.col_values(0))
# 以後、FFTの処理
start = 0
N = 4096 # データ点数
fs = 2000000
# Hanning窓を掛ける
hanning_window = np.hanning(N)
windowed_data = hanning_window * col_array_in_numpy
fft_result = np.fft.fft((windowed_data[start:start + N]))
freq_list = np.fft.fftfreq(N, 1.0 / fs)
# 振幅・位相のスペクトルを求める。
amplitude_spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in fft_result]
phase_spectrum = [rad_to_deg(np.arctan2(int(c.imag), int(c.real))) for c in fft_result]
# 元データの波形を表示する
subplot(311)
plot(range(start, start + N), col_array_in_numpy[start:start + N])
axis([start, start + N, np.min(col_array_in_numpy), np.max(col_array_in_numpy)])
xlabel("time [sample]")
ylabel("amplitude")
# 振幅スペクトルを表示する
subplot(312)
plot(freq_list, amplitude_spectrum, marker='o', linestyle='-')
axis([0, fs / 2, np.min(amplitude_spectrum), np.max(amplitude_spectrum)])
yscale("log") # 縦軸は対数軸にする
xlabel("frequency [Hz]")
ylabel("amplitude spectrum")
# 位相スペクトルを表示する
subplot(313)
plot(freq_list, phase_spectrum, marker='o', linestyle='-')
axis([0, fs / 2, -180, 180])
xlabel("frequency [Hz]")
ylabel("phase spectrum")
show()
結果、以下のようなグラフが生成されます。(振幅スペクトルは、まだ最大値を基準としたdB表示にはしていません)。  |
| 窓関数なしの場合の、入力データ(上段)と、そのFFT結果 |
 |
| Hanning窓を適用させた場合の、入力データ(上段)と、そのFFT結果 |
一部データを抜き出す処理を加える
さて、今回はデータ数Nは4096でしたが、これが異なる点数の時を考えます。この際、考慮すべきことが1つあります。それは、FFTに使うデータの個数です。というのも、FFTではNは2の累乗である(N = 2^n)ことが望ましいです。なので、N≠2^nであった場合、FFTに用いるデータ数は2^nとなるように、一部を抜き出さねばなりません。例えば、データ数N = kである場合、kよりも小さく、一番大きな2の累乗の数だけ、データを抜き出したくなります。実際には、log(2, k)以下の一番大きな整数Nを求めれば、log(2, k) > Nよりk > 2^Nとなります。よく対数の練習問題にあるやつですね。なので、このNを求める処理を加えます。
Pythonで関数として実装すると、以下のようになりました。
def floor_2_pow(value): # 2**N < log2(value) < 2**N+1となるような整数を返す。つまり、valueに一番近い2**N log_2_value = np.log2(value) exponentiation = np.floor(log_2_value) # print "value : %10.24f exponentiation : %10.24f" % (log_2_value, exponentiation) return 2 ** exponentiation
さて、これを書いて実行していたのですが、value = 4096の時、2^11 = 2048が返ってきてしまいました。デバッグして各変数を見てみると、log_2_valueは12.0(float64)ですが、exponentiationが11.0になっていました。floor(n)は、n以下の最大の整数を求める関数ですから、11.0ではなく12.0が入っているはずです。実際に、
exponentiation = np.floor(np.float64(12.0))とすれば、exponentiationは12.0になります。なので、問題はnp.floor()以前、np.log2()にありそうです。
ここで、浮動小数点では丸め誤差などの誤差がよく問題になることを思い出しました。なので、np.log2(4096)の返り値の表示精度を上げてみました(コメントアウトされているprint文の、コメントを外した)。その結果は以下のようになりました。
log_2_value : 11.999999999999998223643161 exponentiation : 11.000000000000000000000000
これを見ると、np.log2(4096)が実は12.0ではなく、11.999999…を返していたようです。デバッガのWatch窓は小数第一位までしか表示しないので、私が12.0だと勝手に勘違いしていたみたいです。問題はここにあったようです。
検証
さて、本当にnp.log2()が悪いんでしょうか。よく分からなかったので、幾つかの値でnp,log2()をチェックしてみました。以下のスクリプト(verification_log2.py)を用いました。# -*- coding: utf-8 -*- __author__ = 'fenrir_' import math import numpy as np if __name__ == "__main__": # np.log2()が返す値を検証する。 for num in range(1, 20): val = 2 ** num print "np.log2(%d) : %10.24f math.log(%d, 2.0) : %10.24f" % (val, np.log2(val), val, math.log(val, 2.0)) # np.log10()が返す値を検証する for num in range(0, 20): print "np.log10(%d) : %10.24f" % (num, np.log10(10 ** num))その結果、以下のような結果が得られました。
np.log2(2) : 1.000000000000000000000000 math.log(2, 2.0) : 1.000000000000000000000000
np.log2(4) : 2.000000000000000000000000 math.log(4, 2.0) : 2.000000000000000000000000
np.log2(8) : 2.999999999999999555910790 math.log(8, 2.0) : 3.000000000000000000000000
np.log2(16) : 4.000000000000000000000000 math.log(16, 2.0) : 4.000000000000000000000000
np.log2(32) : 5.000000000000000000000000 math.log(32, 2.0) : 5.000000000000000000000000
np.log2(64) : 5.999999999999999111821580 math.log(64, 2.0) : 6.000000000000000000000000
np.log2(128) : 6.999999999999999111821580 math.log(128, 2.0) : 7.000000000000000000000000
np.log2(256) : 8.000000000000000000000000 math.log(256, 2.0) : 8.000000000000000000000000
np.log2(512) : 9.000000000000000000000000 math.log(512, 2.0) : 9.000000000000000000000000
np.log2(1024) : 10.000000000000000000000000 math.log(1024, 2.0) : 10.000000000000000000000000
np.log2(2048) : 11.000000000000000000000000 math.log(2048, 2.0) : 11.000000000000000000000000
np.log2(4096) : 11.999999999999998223643161 math.log(4096, 2.0) : 12.000000000000000000000000
np.log2(8192) : 12.999999999999998223643161 math.log(8192, 2.0) : 13.000000000000000000000000
np.log2(16384) : 13.999999999999998223643161 math.log(16384, 2.0) : 14.000000000000000000000000
np.log2(32768) : 15.000000000000000000000000 math.log(32768, 2.0) : 15.000000000000000000000000
np.log2(65536) : 16.000000000000000000000000 math.log(65536, 2.0) : 16.000000000000000000000000
np.log2(131072) : 17.000000000000000000000000 math.log(131072, 2.0) : 17.000000000000000000000000
np.log2(262144) : 18.000000000000000000000000 math.log(262144, 2.0) : 18.000000000000000000000000
np.log2(524288) : 19.000000000000000000000000 math.log(524288, 2.0) : 19.000000000000000000000000
うーん、np.log2()はところどころ期待した結果が得られてないようです。math.log()の方は、ちゃんと出てますね。ちなみに、np.log10()ではこのようなことは起きないようです。
これについては、Stack Overflowでも話題になったみたいです。確かに、浮動小数点の精度は有限ですね。
Python NumPy log2 vs MATLAB - Stack Overflow http://stackoverflow.com/questions/17702065/python-numpy-log2-vs-matlab
対策
この誤差は、np.log2()の求解アルゴリズムに依るものなのでしょうか。丸め誤差でしょうか。今のところ、この原因はよく分からないです(私が無知なので、調べてもない)。とにかく、このままではnp.log2(2^n)≠nとなることがあるようなので、value = 2^nの時はnp.log2()を使わないようにしました。処理としては、「valueが2で割り切れるならば、valueを返す」です。最終的に、コードは以下のようになりました。def is_dividable_by_2(value): while True: if value % 2 != 0: return False elif value == 2: return True else: value /= 2 def floor_2_pow(value): # 2**N < log2(value) < 2**N+1となるような整数を返す。つまり、valueに一番近い2**N if is_dividable_by_2(value): return value else: log_2_value = np.log2(value) exponentiation = np.floor(log_2_value) # print "value : %10.24f exponentiation : %10.24f" % (log_2_value, exponentiation) return 2 ** exponentiation
参考サイト
・高速フーリエ変換（FFT） - 人工知能に関する断創録 http://aidiary.hatenablog.com/entry/20110618/1308367728
numpyでのFFTの仕方が載っています。このサイトは、Pythonでの機械学習やデータ処理など、色々有用な例が豊富にあって、とても勉強になります。
・numpy.floor — NumPy v1.8 Manual http://docs.scipy.org/doc/numpy/reference/generated/numpy.floor.html
・numpy.floor — NumPy v1.8 Manual http://docs.scipy.org/doc/numpy/reference/generated/numpy.floor.html
numpy.floor()のReferenceページです。numpy.floor(2.0)は2.0だって、例に書いてありますね。5回くらい「えっ本当だよね？」って見返してしまいました…。
・ライブラリ：unittest - Life with Python http://www.lifewithpython.com/2014/03/unittest.html
・ライブラリ：unittest - Life with Python http://www.lifewithpython.com/2014/03/unittest.html
Pythonでのテストの書き方が掲載されています。ちなみに、main.pyにある関数をtest.pyで使いたい場合は、from main import *と書けばいいみたいです。
・Creating Tests http://www.jetbrains.com/pycharm/webhelp/creating-tests.html
・Creating Tests http://www.jetbrains.com/pycharm/webhelp/creating-tests.html
私が使っているIDE「PyCharm」のHelpです。テストの作り方が載っています。Ctrl + Shift + Tを押すだけです。楽でいいですね…。
・高速フーリエ変換 - Wikipedia http://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B
結局のところ、FFTのデータ数Nは、2の累乗「でなければいけない」のか、「であれば一番早い(他の値でもよし)」なのか、どっちなのか。こういう厳密なところを全然おさえられていないので、頭の悪さ、無知といったものが露呈してしまいますね…。
・高速フーリエ変換 - Wikipedia http://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B
結局のところ、FFTのデータ数Nは、2の累乗「でなければいけない」のか、「であれば一番早い(他の値でもよし)」なのか、どっちなのか。こういう厳密なところを全然おさえられていないので、頭の悪さ、無知といったものが露呈してしまいますね…。
