高速フーリエ変換からパワースペクトルを求める (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の式を入れればパワースペクトルはもとめられますか？？？？？

←解決時は質問者本人がここをチェックしてください。

戻る