いぬおさんのおもしろ数学実験室

おいしい紅茶でも飲みながら数学、物理、工学、プログラミング、そして読書を楽しみましょう

0でパディングしてFFTしてもよいのか

 FFT高速フーリエ変換)はデータ数が2の累乗のときに使える方法です。DFTにはそういう制限はありません。最近何回かFFTについて記事を書きました。それによると2の累乗が仮定されているのです。しかし、実際にSciPyのFFTに500データを渡しても無事に(?)FFTを実行してくれます。どういう仕組みになっているのか分かりませんが……。2の累乗でないと実行速度は遅くなる、という記事も見かけます。

 サイトによっては「2の累乗でないときには0でパディングするとよい」などと書いてあるんですが、どうもぼくにはこれはまずそうな気がします。例えば500データしかなければ12データ分、0を付け足してFFTを実行する、というわけです。でも、これではもちろん対象データの周期は一方は500データ、他方は512となってしまいますし(当たり前……)、そもそも異なるデータ(0が余分に混ざったデータ)を分析することになります。状況が許すなら、多目に600データ用意しておき、そのうちの512データを取り出してFFTする方が良さそうに思えます。実際にどちらがよいのかぼくには分かりませんが。

 今回、ひとつ実験をやってみました。500データをFFTしたいのですが、そのままFFTを実行(さきに書いたとおり、SciPyは無事に実行してくれる)するのと、0で12データだけパディングしてFFTをするのとで、得られる分析結果がどう異なるのか見てみます。

 以降、ザッと説明するだけなのでFFTをやったことがない人は分からないかも。ごめんなさい!!

 サンプリング周波数は600Hzとしておきます。データ数が500なら、基本周波数は600/500=1.2Hz、データ数が512なら基本周波数は600/512≒1.172Hzです。FFTではこの基本周波数の整数倍の周波数成分が求まるのでした。

#=====================

import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt

fs = 600 #サンプリング周波数 600Hz
t = np.arange(0, 4, 1/fs) #4秒分の時刻の配列 1/600秒刻み
tt500 = np.arange(0, 500)
tt500 = tt500 * (fs / 500) #周波数分析のグラフの横軸
tt512 = np.arange(0, 512)
tt512 = tt512 * (fs / 512) #周波数分析のグラフの横軸

y = np.sin(2*np.pi*50*t) #50Hzの信号
for i in range(500, 512): #y[500]~y[511]を0にする
    y[i] = 0
yf = sp.fftpack.fft(y[:500]) #A 通常のデータ(500個)をFFT
plt.plot(tt500, np.abs(yf))
plt.show()
yf1 = sp.fftpack.fft(y[:512]) #B 12データ分、0でパディングしてからFFT
plt.plot(tt512, np.abs(yf1))
plt.show()

#=====================

結果は次の通り。

 

#A 500データ、パディングなし

f:id:Inuosann:20220301234731p:plain

 

#B 512データ、パディングあり

f:id:Inuosann:20220301234819p:plain

どちらもサンプリング周波数は600Hzなので、グラフの右端は600Hzになっています。分析対象のデータは50Hzの信号ですが、Aは1.2Hz、Bは1.172Hzの整数倍の周波数成分を用いて元データを表さなければならず、単純にピークがひとつだけ……とはなっていません。しかしそれでもとにかく各周波数の成分はグラフのように分かるので、2つを比較はできます。結果、似てはいるけれど当然異なるグラフです。

 

 違う話ですが、ウィンドウ(窓関数)を使うのも近い部分があって、窓関数を使うと元のデータが変わってしまいます。使った方がいいケースもあるのでしょうが、それを周波数分析しても当然元データを分析したのとは違う結果が出ます。

 今回、実はあまり自信がないんですが敢えて記事にしています。今、FFTについて整理しているところです。分かっていたつもりでも曖昧だったり、そんなことばかり。賢い人は違うのでしょうが、これが勉強というものなのかも知れません……。