|
［ 簡単な説明 ］ 高速フーリエ変換ライブラリの使用例です。 test16.c は、定数 CPX の定義の有無で、Complex 型入力と、double 型実数部・虚数部別入力を選択できるように、コーディングしています。１次元、２次元のルーチンでそれぞれＦＦＴ、逆ＦＦＴを行ない、結果を表示します。 test17.c は、fft1( )ルーチンを用いて、ベンチマークを行っています。 |
|
（先頭に戻る）
/* test16.c */
#include <stdio.h>
#include "sslib.h"
/*#define CPX /* Complex mode */
#ifndef CPX
void print(double *ar, double *ai, int n);
void print2(double *ar, double *ai, int n, int nmax);
#else
void printx(Complex *a, int n);
void print2x(Complex *a, int n, int nmax);
#endif
#ifndef CPX
void print(double *ar, double *ai, int n)
{
int i;
double *p, *q;
printf(" ar ai\n");
for(i = 0, p = ar, q = ai; i < n; i++) printf("%14.7e %14.7e\n", *p++, *q++);
}
#else
void printx(Complex *a, int n)
{
int i;
Complex *p;
printf(" ar ai\n");
for(i = 0, p = a; i < n; i++, p++) printf("%14.7e %14.7e\n", p->r, p->i);
}
#endif
#ifndef CPX
void print2(double *ar, double *ai, int n, int nmax)
{
int i, j, k, l;
double *p;
printf("< ar >\n");
for(i = k = 0, p = ar; i < n; i++)
{
for(j = 0; j < nmax; j++, k++) printf("%14.7e ", *p++);
putchar('\n');
}
printf("< ai >\n");
for(i = k = 0, p = ai; i < n; i++)
{
for(j = 0; j < nmax; j++, k++) printf("%14.7e ", *p++);
putchar('\n');
}
}
#else
void print2x(Complex *a, int n, int nmax)
{
int i, j, k, l;
Complex *p;
printf("< Real part >\n");
for i = k = 0, p = a; i < n; i++)
{
for(j = 0; j < nmax; j++, k++, p++) printf("%14.7e ", p->r);
putchar('\n');
}
printf("< Imaginary part >\n");
for(i = k = 0, p = a; i < n; i++)
{
for(j = 0; j < nmax; j++, k++, p++) printf("%14.7e ", p->i);
putchar('\n');
}
}
#endif
int main(void)
{
static double ar[8] = { 0., 0., 0., 1., 1., 0., 0., 0.};
static double ai[8] = { 0., 0., 0., 0., 0., 0., 0., 0.};
static double ar2[16] = { 0., 0., 0., 0., 0., 1., 1., 0.,
0., 1., 1., 0., 0., 0., 0., 0.};
static double ai2[16] = { 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0.};
Complex a[8], a2[16];
int flag, i, iter, j, k, n, nmax;
iter = 0;
n = 8;
#ifndef CPX
for(i = 0; i < n; i++) a[i] = tocomplex(ar[i], ai[i]);
print(ar, ai, n);
#else
printx(a, n);
#endif
/* forward FFT */
flag = 0;
#ifndef CPX
fft1(ar, ai, n, iter, flag);
print(ar, ai, n);
#else
fft1x(a, n, iter, flag);
printx(a, n);
#endif
/* reverse FFT */
flag = 1;
#ifndef CPX
fft1(ar, ai, n, iter, flag);
print(ar, ai, n);
#else
fft1x(a, n, iter, flag);
printx(a, n);
#endif
n = nmax = 4;
for(i = k = 0; i < n; i++) for(j = 0; j < n; j++, k++) a2[k] = tocomplex(ar2[k], ai2[k]);
#ifndef CPX
print2(ar2, ai2, n, nmax);
#else
print2x(a2, n, nmax);
#endif
flag = 0;
#ifndef CPX
fft2(ar2, ai2, n, nmax, flag);
print2(ar2, ai2, n, nmax);
#else
fft2x(a2, n, nmax, flag);
print2x(a2, n, nmax);
#endif
flag = 1;
#ifndef CPX
fft2(ar2, ai2, n, nmax, flag);
print2(ar2, ai2, n, nmax);
#else
fft2x(a2, n, nmax, flag);
print2x(a2, n, nmax);
#endif
return 1;
}
|
（先頭に戻る）
ar ai
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
1.0000000e+00 0.0000000e+00
1.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
ar ai
2.0000000e+00 0.0000000e+00
-1.7071068e+00 -7.0710678e-01
1.0000000e+00 1.0000000e+00
-2.9289322e-01 -7.0710678e-01
0.0000000e+00 0.0000000e+00
-2.9289322e-01 7.0710678e-01
1.0000000e+00 -1.0000000e+00
-1.7071068e+00 7.0710678e-01
ar ai
0.0000000e+00 0.0000000e+00
1.9630836e-17 -2.6325784e-17
0.0000000e+00 0.0000000e+00
1.0000000e+00 -5.1675786e-17
1.0000000e+00 0.0000000e+00
-1.9630836e-17 2.9185367e-17
0.0000000e+00 0.0000000e+00
0.0000000e+00 4.8816203e-17
< ar >
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
0.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00
0.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
< ai >
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
< ar >
4.0000000e+00 -2.0000000e+00 0.0000000e+00 -2.0000000e+00
-2.0000000e+00 -1.1102230e-16 0.0000000e+00 2.0000000e+00
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
-2.0000000e+00 2.0000000e+00 0.0000000e+00 1.1102230e-16
< ai >
0.0000000e+00 -2.0000000e+00 0.0000000e+00 2.0000000e+00
-2.0000000e+00 2.0000000e+00 0.0000000e+00 2.2204460e-16
0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00
2.0000000e+00 0.0000000e+00 0.0000000e+00 -2.0000000e+00
< ar >
1.3877788e-17 2.7755576e-17 -1.3877788e-17 -2.7755576e-17
0.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00
0.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00
0.0000000e+00 1.3877788e-17 0.0000000e+00 -1.3877788e-17
< ai >
0.0000000e+00 1.3877788e-17 0.0000000e+00 -1.3877788e-17
-1.4297916e-18 -3.2044950e-17 -1.4297916e-18 2.9185367e-17
0.0000000e+00 -3.0615159e-17 0.0000000e+00 3.0615159e-17
-1.2447996e-17 2.9185367e-17 1.5307579e-17 -2.6325784e-17
|
（先頭に戻る）
/* test17.c */
#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <time.h>
#include "sslib.h"
#define NN 1024 /* number of FFT sampling */
#define N 256 /* number of data */
static double f1[N], f2[N]; /* analysis data */
void sines_2(double *x, int n, double t1, double t2, double t3, double d1, double d2, double d3)
{
int i;
double ii, *p, pi2, sc1, sc2, sc3;
pi2 = 2. * M_PI;
sc1 = pi2 / t1;
sc2 = pi2 / t2;
sc3 = pi2 / t3;
d1 *= pi2;
d2 *= pi2;
d3 *= pi2;
for(i = 0, ii = 0., p = x; i < n; i++, ii++)
*p++ = (sin(ii * sc1 + d1) + sin(ii * sc2 + d2) / 2. + sin(ii * sc3 + d3) / 4.) / 1.75;
}
int main(void)
{
static double ar[8] = { 0., 0., 0., 1., 1., 0., 0., 0. };
static double ai[8] = { 0., 0., 0., 0., 0., 0., 0., 0. };
int i, iter, n, l;
double err, maxerr, *x1, *y1;
clock_t t1, t2;
iter = 3;
n = 8;
printf(" ar ai\n");
for(i = 0; i < n; i++) printf("%14.7e %14.7e\n", ar[i], ai[i]);
fft1(ar, ai, n, iter, 0);
printf("<fourier transform>\n ar ai\n");
for(i = 0; i < n; i++) printf("%14.7e %14.7e\n", ar[i], ai[i]);
fft1(ar, ai, n, iter, 1);
printf("<reverse fourier transform>\n ar ai\n");
for(i = 0; i < n; i++) printf("%14.7e %14.7e\n", ar[i], ai[i]);
printf("\n<Time Trial>\n");
n = i = NN;
l = 0;
while(i > 1)
{
i /= 2;
l++;
}
x1 = (double *)malloc(n * sizeof(double));
y1 = (double *)malloc(n * sizeof(double));
if(x1 == NULL || y1 == NULL)
{
fprintf(stderr, "Error : Out of memory!\n");
exit(1);
}
sines_2(x1, n, 100., 25., 10., 0., 0.5, -0.25);
for(i = 0; i < n; i++) y1[i] = 0.;
t1 = clock();
for(i = 0; i < 500; i++)
{
fft1(x1, y1, n, l, 0);
fft1(x1, y1, n, l, 1);
}
t2 = clock();
printf("fft1 (forward & reverse 50 times) : %ld [sec]\n", t2 - t1);
free((char *)x1);
free((char *)y1);
return 1;
}
|
（先頭に戻る）
ar ai
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
1.0000000e+00 0.0000000e+00
1.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
0.0000000e+00 0.0000000e+00
<fourier transform>
ar ai
2.0000000e+00 0.0000000e+00
-1.7071068e+00 -7.0710678e-01
1.0000000e+00 1.0000000e+00
-2.9289322e-01 -7.0710678e-01
0.0000000e+00 0.0000000e+00
-2.9289322e-01 7.0710678e-01
1.0000000e+00 -1.0000000e+00
-1.7071068e+00 7.0710678e-01
<reverse fourier transform>
ar ai
0.0000000e+00 0.0000000e+00
1.9630836e-17 -2.6325784e-17
0.0000000e+00 0.0000000e+00
1.0000000e+00 -5.1675786e-17
1.0000000e+00 0.0000000e+00
-1.9630836e-17 2.9185367e-17
0.0000000e+00 0.0000000e+00
0.0000000e+00 4.8816203e-17
<Time Trial>
fft1 (forward & reverse 50 times) : 163 [sec]
|
戻る（メインメニュー）
戻る（科学技術計算メニュー）
注意事項