|
|
FFTWを使う
高速なFFTライブラリーFFTWを補助ルーチン(fftw.dll)を用いて呼び出します。
下記よりコンパイル済みFFTWライブラリーFFTW 3.3.4(32bit版)をダウンロードしてください
実行には、libfftw3-3.dllを使用します。
libfftw3-3.dllをbasic.exeと同じフォルダに入れてください。
(libfftw3f-3.dll 及び libfftw3l-3.dll は使用していません)
OPTION ARITHMETIC NATIVE
OPTION CHARACTER BYTE
OPTION BASE 0
!'LET N=2^20
LET N=1000000 !'2のペキ乗(2^n)でなくてもいい。
!'DIM RR(N),II(N)
LET T=TIME
LET RE$=REPEAT$(CHR$(0),N*8)
LET IM$=REPEAT$(CHR$(0),N*8)
FOR I=0 TO N-1
!' LET RR(I)=I
!' LET II(I)=0
!' LET RE$(8*I+1:8*I+8)=PACKDBL$(RR(I)) !'実部
!' LET IM$(8*I+1:8*I+8)=PACKDBL$(II(I)) !'虚部
LET RE$(8*I+1:8*I+8)=PACKDBL$(I)
NEXT I
PRINT TIME-T;"秒"
LET T=TIME
CALL FFTW(N,RE$,IM$,0) !'変換
PRINT TIME-T;"秒"
LET T=TIME
CALL FFTW(N,RE$,IM$,-1) !'逆変換
PRINT TIME-T;"秒"
!'FOR I=0 TO N-1
!' LET RR(I)=INT(UNPACKDBL(RE$(8*I+1:8*I+8))/N+.5)
!' LET II(I)=INT(UNPACKDBL(IM$(8*I+1:8*I+8))/N+.5)
!' PRINT RR(I);II(I)
!'NEXT I
END
EXTERNAL SUB FFTW(N,RE$,IM$,SW)
OPTION ARITHMETIC NATIVE
ASSIGN "fftw.dll","fftwsub"
END SUB
--------------------------------------------------------------------------------------
fftw.cpp
#include <cstdlib>
#include <fftw3.h>
using namespace std;
#pragma comment(lib, "libfftw3-3.lib")
// #pragma comment(lib, "libfftw3f-3.lib")
// #pragma comment(lib, "libfftw3l-3.lib")
extern "C" __declspec(dllexport) void fftwsub(int n,double *re,double *im,int sw)
{
int i;
fftw_complex *a, *b;
fftw_plan plan;
a= (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * n);
b= (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * n);
if (a==NULL || b==NULL) exit(1);
if (sw==0)
plan = fftw_plan_dft_1d(n, a, b, FFTW_FORWARD , FFTW_ESTIMATE);
else
plan = fftw_plan_dft_1d(n, a, b, FFTW_BACKWARD, FFTW_ESTIMATE);
for(i = 0;i<n;i++){
a[i][0] = re[i];
a[i][1] = im[i];
}
fftw_execute(plan);
for(i = 0;i<n;i++){
re[i]=b[i][0];
im[i]=b[i][1];
}
fftw_destroy_plan(plan);
fftw_free(a);
fftw_free(b);
}
|
|