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回繰り返すらしいのですが教えてください。

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

戻る