FFTプログラムでパワースペクトルはできたんですが周波数を入れるには？ (ID:65403)

名前

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

本文 /*--------------------------------------------- *プログラム名 ： パワースペクトル *ファイル名 ： ffttest.c *バージョン ： up *作成日 ： 1994/4 *-------------------------------------------- */ #include <stdio.h> #include <process.h> #include <conio.h> #include <ctype.h> #include <stdlib.h> #include <math.h> #define PI 3.14159265358979323826 #define N 256 /* データの個数 */ int fft( int, double *, double * ); int main( void ) { int j,i, n=256; static double x[N], y[N],P[N]; /* data */ FILE *fp1; for(j=0;j<256;j++){ x[j]=sin(2*PI*5*j*0.01)+sin(2*PI*20*j*0.01)+sin(2*PI*35*j*0.01); } for(j=0;j<256;j++) y[j]=0; fft( n, x, y ); for(i=1;i<=128;i++){ P[i]=0.00390625*(x[i]*x[i]+y[i]*y[i]); printf("%f",y[i]); putchar('\n'); } fp1=fopen("fu-rie","w"); fprintf(fp1,"fu-rie \n"); for(i=1;i<=128;i++){ fprintf(fp1,"%3d%15f \n",i,P[i]); fprintf(fp1,"\n"); } fclose(fp1); return( 0 ); } /*----------------------------------------------------------------------- * 関数名: int fft( int n, double x[], double y[] ) * 機 能: ＦＦＴ計算 * 入 力: n データ数 * : x データ配列（実部） * : y データ配列（虚部） * 出 力: 0 * 作成日: 1991/09/24 *---------------------------------------------------TextProgram--------*/ static void make_sintbl( int n, double 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.0; 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]; } static void make_bitrev( int n, int bitrev[] ) { int i, j, k, n2; n2 = n / 2; i = j = 0; for ( ; ; ) { bitrev[i] = j; if ( ++i >= n ) break; k = n2; while ( k <= j ) { j -= k; k /= 2; } j += k; } } int fft( int n, double x[], double y[] ) { static int last_n = 0; /* 前回呼出時のn */ static int *bitrev = NULL; /* ビット反転表 */ static double *sintbl = NULL; /* 三角関数表 */ int i, j, k, ik, h, d, k2, n4, inverse; double t, s, c, dx, dy; /* 準備 */ 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(double) ); bitrev = malloc( n * sizeof(int) ); if ( sintbl == NULL || bitrev == NULL ) { fprintf( stderr, "記憶領域不足\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]; 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; } } return( 0 ); } の fp1=fopen("fu-rie","w"); fprintf(fp1,"fu-rie \n"); for(i=1;i<=128;i++){ fprintf(fp1,"%3d%15f \n",i,P[i]);の 4行目のiの位置に周波数がくるらしいのですが、 周波数がどうだせばいいかわかりません。for分で128回繰り返すらしいのですが教えてください。