ニューメリカル・レシピのFFTプログラムをExcelのVBAに移植
ExcelのVBAで、サブルーチンに配列をそのまま渡せることを(いまさら)初めて知った。じゃあ、普通に配列渡して高速フーリエ変換(FFT)とか計算できるな、ってことで手元にあったニューメリカル・レシピ(in C)のCプログラムを移植してみた。
プログラムはこれです。↓
使い方は、上のを標準モジュールとして読み込んだ後、配列data()を使って
call FFT(data nn, isign)
と呼ぶんだけど、ちょっと使い方が変わっていて
data(1) = データ1のreal part,
data(2) = データ1のimag part
data(3) = データ2のreal part
data(4) = データ2のimag part
・・・
とか並べます。データnnまで(つまり配列としては2*nnは必要)。当然nnは2^Mの形です。
isignは+1ならフーリエ変換、-1なら逆フーリエ変換ですよ。
まあ、
Option Explicit
Private Sub CommandButton1_Click()
Dim data(16) As Double
Dim i As Integer
data(1) = 1#
data(2) = 0#
data(3) = 1#
data(4) = 0#
data(5) = 1#
data(6) = 0#
For i = 7 To 16
data(i) = 0
Next i
Call FFT(data, 8, 1)
For i = 1 To 16
Worksheets("Sheet1").Cells(i, 5) = data(i)
Next i
End Sub
Sub FFT(data() As Double, nn As Integer, isign As Integer)
Dim n As Integer, mmax As Integer, m As Integer, j As Integer, istep As Integer, i As Integer
Dim wtemp As Double, wr As Double, wpr As Double, wpi As Double, wi As Double, theta As Double
Dim tempr As Double, tempi As Double
Dim dummy As Double
n = 2 * nn
j = 1
For i = 1 To n - 1 Step 2
If j > i Then
dummy = data(i)
data(i) = data(j)
data(j) = dummy
dummy = data(i + 1)
data(i + 1) = data(j + 1)
data(j + 1) = dummy
End If
m = n / 2
While (m >= 2 And j > m)
j = j - m
m = m / 2
Wend
j = j + m
Next i
mmax = 2
While (n > mmax)
istep = mmax * 2#
theta = isign * (2# * 4# * Atn(1) / mmax)
wtemp = Sin(0.5 * theta)
wpr = -2# * wtemp * wtemp
wpi = Sin(theta)
wr = 1#
wi = 0#
For m = 1 To mmax - 1 Step 2
For i = m To n Step istep
j = i + mmax
tempr = wr * data(j) - wi * data(j + 1)
tempi = wr * data(j + 1) + wi * data(j)
data(j) = data(i) - tempr
data(j + 1) = data(i + 1) - tempi
data(i) = data(i) + tempr
data(i + 1) = data(i + 1) + tempi
Next i
wtemp = wr
wr = wtemp * wpr - wi * wpi + wr
wi = wi * wpr + wtemp * wpi + wi
Next m
mmax = istep
Wend
End Sub
こんな感じで使えばいいかな。間違っていても責任はもてまへん。。。
« ほとんど至るところほとんど至るところといわなきゃいけない。 | トップページ | 仁左衛門の湯 »
「パソコン・インターネット」カテゴリの記事
- ソフトバンクが1.7GHz帯を3Gで使えなくしたので、持っているポケットWiFi GL01Pが使えなくなり、代わりに607HWが送られてきた。(2018.02.19)
- Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その7:メルセンヌツイスタ(2018.02.13)
- Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その6:ルンゲクッタ8次 DOP853(2018.02.12)
- Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その5:フィッティング(非線形)(2018.02.09)
- Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その4:フィッティング(線形)(2018.02.08)
「学問・資格」カテゴリの記事
- 信頼性(寿命)評価に使われる”ワイブル分布のパラメータの算出”をカシオの高精度計算サイトに自作式として作った。(2018.04.20)
- ポアソン分布の例で馬に蹴られて死んだ兵士の例がよく出るが、原典見て実際に計算してみた。(2018.04.19)
- #仮面ライダービルド 第31話の話数を表す数式はメルセンヌ素数!M5=2^5-1=31(2018.04.15)
- ラプラス方程式をある差分化で数値計算すると高木関数になる例(F((2i+1)/^(2k+1))=0.5*(F(i/2^k)+F((i+1/2^k)) + Ck)(2018.04.11)
- 解のない偏微分方程式(2013.01.23)
トラックバック
この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/512682/44381523
この記事へのトラックバック一覧です: ニューメリカル・レシピのFFTプログラムをExcelのVBAに移植:
コメント