ニューメリカル・レシピの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
こんな感じで使えばいいかな。間違っていても責任はもてまへん。。。
« ほとんど至るところほとんど至るところといわなきゃいけない。 | トップページ | 仁左衛門の湯 »
「パソコン・インターネット」カテゴリの記事
- 高周波回路シミュレータQucsStudioを使ってみる(その3)Mixed Mode S parameterを計算(2018.10.12)
- 高周波回路シミュレータQucsStudioを使ってみる(その2)SパラメータのTouchStoneフォーマットで出力するには?(2018.10.11)
- 高周波回路シミュレータQucsStudioを使ってみる(その1)まずは何をさておきμの文字化けだけには注意。(2018.10.10)
- 円の弧長,弦長,矢高,半径のどれか2つを与えて残りを計算(カシオの高精度計算サイト自作式)で180°以上、複数解に対応。(2018.10.09)
- macbook proをmacOS Mojaveにアップデート。せっかくなんでダークモードにしてみる。(2018.09.26)
「学問・資格」カテゴリの記事
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
- フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。(2021.02.23)
- 「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP! インピーダンスの大きさと位相がグラフ化できる。(2021.02.12)
コメント