FFT on Groovy

すいませんちょい遅れましたが12/21のG* Advent Calendarです。
FFT(高速フーリエ変換)というかDFT(離散フーリエ変換)やります。

使うライブラリ

JTransforms https://github.com/wendykierp/JTransforms

入力波形

input

[0.0, 0.049067674327418015, 0.0980171403295606, 0.14673047445536175, 0.19509032201612825, 0.24298017990326387, 0.29028467725446233, 0.33688985339222005, 0.3826834323650898, 0.4275550934302821, 0.47139673682599764, 0.5141027441932217, 0.5555702330196022, 0.5956993044924334, 0.6343932841636455, 0.6715589548470183, 0.7071067811865475, 0.7409511253549591, 0.7730104533627369, 0.8032075314806448, 0.8314696123025452, 0.8577286100002721, 0.8819212643483549, 0.9039892931234433, 0.9238795325112867, 0.9415440651830208, 0.9569403357322089, 0.970031253194544, 0.9807852804032304, 0.989176509964781, 0.9951847266721968, 0.9987954562051724, 1.0, 0.9987954562051724, 0.9951847266721969, 0.989176509964781, 0.9807852804032304, 0.970031253194544, 0.9569403357322089, 0.9415440651830208, 0.9238795325112867, 0.9039892931234434, 0.881921264348355, 0.8577286100002721, 0.8314696123025453, 0.8032075314806449, 0.7730104533627371, 0.740951125354959, 0.7071067811865476, 0.6715589548470186, 0.6343932841636455, 0.5956993044924335, 0.5555702330196022, 0.5141027441932218, 0.47139673682599786, 0.42755509343028203, 0.3826834323650899, 0.33688985339222033, 0.2902846772544624, 0.24298017990326407, 0.1950903220161286, 0.1467304744553618, 0.09801714032956083, 0.049067674327417966, 1.2246467991473532E-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

FFTコード

dft_sample.groovy

@Grab("com.github.wendykierp:JTransforms:3.1")

import org.jtransforms.fft.DoubleFFT_1D

def numPoints = 64

// numPoints + 1個のデータ
def x = []
0.step(1 + 1.0 / numPoints, 1.0 / numPoints) { x << it }

// FFTに用いるsin波
def f = x.collect { Math.sin( it * Math.PI) }

// FFT用にサイズを2倍にし、0で埋める
f << [0] * (numPoints - 1)
f = f.flatten()
println f

// FFT用にdoubleに詰めなおす
double[] fDouble = new double[numPoints * 2]
f.eachWithIndex { it, i ->
  fDouble[i] = it
}

def fft = new DoubleFFT_1D(numPoints)
fft.realForward(fDouble)
println fDouble

出力(パワースペクトル)

正規化とか特になにもしれいなのであれですが

output

[40.73548387208331, -0.024548622108927276, -13.589407314865438, 8.881784197001252E-16, -2.7244417193377584, 4.440892098500626E-16, -1.172321899279049, 5.551115123125783E-16, -0.6549698903522008, 5.551115123125783E-17, -0.41983327155349126, 6.38378239159465E-16, -0.29324623849482734, 1.6653345369377348E-16, -0.21732034328985128, 4.718447854656915E-16, -0.1682193329044574, -1.2559397966072083E-15, -0.1346525297928969, 9.020562075079397E-17, -0.11070596112458637, -6.626643678231403E-16, -0.09303715106699177, -3.382710778154774E-16, -0.07964083223305846, -9.922618282587337E-16, -0.06925466172876894, -1.734723475976807E-17, -0.06105183635156024, 2.252972114424878E-16, -0.05447253054678732, -1.4181364416110398E-16, -0.04912684976946835, -0.0, -0.04473660722209264, -1.4181364416110398E-16, -0.04109913630139779, 5.583641188300348E-16, -0.03806409939444906, 1.491862189340054E-16, -0.035518222343737055, -9.922618282587337E-16, -0.03337500672361353, 9.107298248878237E-16, -0.03156765146552513, 5.030698080332741E-16, -0.030044092954215715, -5.620504062164855E-16, -0.028763473809871924, -1.2559397966072083E-15, -0.02769359446339291, -6.245004513516506E-16, -0.02680905360195942, 1.1102230246251565E-16, -0.026089880411642563, -5.551115123125783E-16, -0.025520524589636717, 5.551115123125783E-17, -0.025089112000558356, -7.771561172376096E-16, -0.024786902307144976, -4.440892098500626E-16, -0.024607904706760664, 0.0, 1.2246467991473532E-16, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

グラフにするとこんな感じ。波数が小さいところに全部集まってるので、FFT自体はできていそう。

補足

• スレッド数指定などを省略したので、性能を出したい場合はそのあたりを追加したほうが良さそう。テストコードが参考になるはず https://github.com/wendykierp/JTransforms/blob/master/src/test/java/org/jtransforms/fft/DoubleFFT_1DTest.java
• double[] への変換は、Java8のlambdaにGroovyが対応したらもっと簡単になりそう http://stackoverflow.com/questions/6018267/how-to-cast-from-listdouble-to-double-in-java の回答二つ目とか。

GroovyでFFTしてみた感想

最初はFFTW使おうとしていたので結構苦戦しましたが、使いやすいライブラリが分かってみると簡単！