高速フーリエ変換からパワースペクトルを求める (ID:52778)

名前

ホームページ(ブログ、Twitterなど)のURL (省略可)

本文 /***** FFT & IFFT *****/ #include <stdio.h> #include <stdlib.h> #include <math.h> /* #define PI 3.14159265358979323846 */ /* sin table */ static void make_sintbl(n, sintbl) int n; float *sintbl; { int i, n2, n4, n8; double c, s, dc, ds, t; n2 = n / 2; n4 = n / 4; n8 = n / 8; t = sin(PI / n); dc = 2 * t * t; ds = sqrt(dc * (2 - dc)); t = 2 * dc; c = sintbl[n4] = 1; s = sintbl[0] = 0; for (i = 1; i < n8; i++) { c -= dc; dc += t * c; s += ds; ds -= t * s; sintbl[i] = s; sintbl[n4 - i] = c; } if (n8 != 0) sintbl[n8] = sqrt(0.5); for (i = 0; i < n4; i++) sintbl[n2 - i] = sintbl[i]; for (i = 0; i < (n2 + n4); i++) sintbl[i + n2] = - sintbl[i]; } /* Bit Reversal */ /* FFT main routine */ int fft(n, x, y) int n; float *x; float *y; { static int last_n = 0; static int *bitrev = NULL; static float *sintbl = NULL; int i, j, k, ik, h, d, k2, n4, inverse; float t, s, c, dx, dy; /* preparation */ if (n < 0) { n = -n; inverse = 1; } else inverse = 0; n4 = n / 4; if (n != last_n || n == 0) { last_n = n; if (sintbl != NULL) free(sintbl); if (bitrev != NULL) free(bitrev); if (n == 0) return 0; sintbl = malloc((n + n4) * sizeof(float)); bitrev = malloc(n * sizeof(int)); if (sintbl == NULL || bitrev == NULL) { fprintf(stderr, "Insufficient Memory\n"); return 1; } make_sintbl(n, sintbl); make_bitrev(n, bitrev); } for (i = 0; i < n; i++) { j = bitrev[i]; if (i < j) { t = x[i]; x[i] = x[j]; x[j] = t; t = y[i]; y[i] = y[j]; y[j] = t; } } for (k = 1; k < n; k = k2) { h = 0; k2 = k + k; d = n / k2; for (j = 0; j < k; j++) { c = sintbl[h + n4]; if (inverse) s = - sintbl[h]; else s = sintbl[h]; for (i = j; i < n; i += k2) { ik = i + k; dx = s * y[ik] + c * x[ik]; dy = c * y[ik] - s * x[ik]; x[ik] = x[i] - dx; x[i] += dx; y[ik] = y[i] - dy; y[i] += dy; } h += d; } } if (! inverse) for (i = 0; i < n; i++) { x[i] /= n; y[i] /= n; } return 0; } /*** Main Program ***/ #define N 256 int main(argc, argv) int argc; char *argv[]; { FILE *fp; int i; static float x1[N], y1[N], x2[N], y2[N], x3[N], y3[N]; if ((fp = fopen(argv[1], "r")) ==NULL){ printf ("File Open Error \"%s \" file.\n", argv[1]); exit(1); } for (i = 0; i < N; i++) { fscanf(fp, "%f", &x1[i]); x2[i] = x1[i] ; y1[i] = y2[i] = 0; } if (fft(N, x2, y2)) return EXIT_FAILURE; for (i = 0; i < N; i++) { x3[i] = x2[i]; y3[i] = y2[i]; } if (fft(-N, x3, y3)) return EXIT_FAILURE; printf(" Original data Fourier Transformed Inverse Transformed\n"); for (i = 0; i < N; i++) printf("%4d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n", i, x1[i], y1[i], x2[i], y2[i],log10(x2[i]*x2[i]+y2[i]y2[i])); return EXIT_SUCCESS; } 下から２行目のprintf()にlogの式を入れればパワースペクトルはもとめられますか？？？？？