ニューメリカル・レシピの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
こんな感じで使えばいいかな。間違っていても責任はもてまへん。。。
« ほとんど至るところほとんど至るところといわなきゃいけない。 | トップページ | 仁左衛門の湯 »
「パソコン・インターネット」カテゴリの記事
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- Gemini CLIが使えるようになっていたので早速VSCodeのターミナルから使って、JavaScriptで連立一次方程式を計算するコードを書いてもらった。普通にガウスの消去法で計算するhtmlを作ってくれた。(2025.06.27)
- Google ColabのJulia言語で1次元のGray-Scottモデル(∂u/∂t=u²v-(F+k)u+Du∂²u/∂x²,∂v/∂t=-u²v+F(1-v)+Dv∂²v/∂x²)を計算してパルスが次々分裂する様子を見る。空間6次の差分、時間8次のルンゲクッタ法で計算。(2025.07.08)
「学問・資格」カテゴリの記事
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
- 高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事(2025.07.02)
- 高周波・RFニュース 2025年6月30日 QualcommがAIを用いた6Rxアンテナ解説、Next G Allianceと日本のXGMFが5G,6Gで協力、5G Americasが25Q1で5G加入者増加と発表、TechInsigtsがHuawei Pura 80 Pro+分解、Qorvoが5-7GHzのWi-Fi 7 FEM発表(2025.06.30)
- Gemini CLIが使えるようになっていたので早速VSCodeのターミナルから使って、JavaScriptで連立一次方程式を計算するコードを書いてもらった。普通にガウスの消去法で計算するhtmlを作ってくれた。(2025.06.27)
コメント