FFT その7 - Maximaでこうぞうりきがく

みんな大好きFFTのお話です
以前MaximaのFFT関数を紹介しましたが，そもそもMaximaは数式処理システムなので大量の数値計算は苦手です（´・ω・｀)
なので今回はC言語で組んだFFTプログラムを使ってみたいと思います

FFT本体のサブルーチンは奥村晴彦著「C言語による最新アルゴリズム事典」に記載のものを使わせて頂きました
サポートサイトでサブルーチンのソースコードを配布されています（make_sintbl, make_bitrev, fftの3つです）

FFTにはCooley-Tukey型FFT アルゴリズムを用いています
入力データのサンプル数は2のべき乗として下さい
計算結果をパイプ処理でgnuplotに出力します
gnuplotフォルダはc:\直下に置かれていることを前提とします(それ以外はpathを通して下さい)
gnuplot-win32のver4.4.3で動作確認しました

以下にソースコードを示します（make_sintbl, make_bitrev, fftについては省略します）
【fft.c】

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
#define NMAX 65536
#define SMAX 1024
#define GNUPLOT_PATH "c:/gnuplot/binary/pgnuplot.exe"
FILE *fp1, *fp2, *fp3, *gp;

int   N;
float X1[NMAX], Y1[NMAX];
float X2[NMAX], Y2[NMAX];
float DT, DF;

void output_fft(void);
void gnuplot(void);

int main(void)
{
	if(!input_fft()) return EXIT_FAILURE;

	else if(fft(N, X2, Y2)) return EXIT_FAILURE;

	else{
		output_fft();
		gnuplot();
		return EXIT_SUCCESS;
	}
}

int input_fft(void)
{
	int	 i = 0, j;
	char buff[SMAX], *tok;

	/*  fft.dat  */
	if((fp1 = fopen("fft.dat", "r")) == NULL){
		printf("can't open input data file\n");
		getch();
		exit(1);
	}
	fgets(buff, SMAX, fp1);
	tok = strtok(buff, "\n");
	DT = atof(tok);

	while(!feof(fp1)){
		fgets(buff, SMAX, fp1);
		tok = strtok(buff, "\t");
		if(tok != NULL) X1[i] = atof(tok);
		else break;
		tok = strtok(NULL, "\n");
		if(tok != NULL) Y1[i] = atof(tok);
		else break;
		i++;
	}
	fclose(fp1);

	/*  error check  */
	if((i & (i-1))){
		printf(" error : %d is not powers of 2\n", i);
		getch();
		return 0;
	}
	N = i;
	for(j = 0; j < i; j++) {
		X2[j] = X1[j];
		Y2[j] = Y1[j];
	}
	DF = 1/(i*DT);
	return 1;
}

void output_fft(void)
{
	int i;	char *s[7];

	s[0] = "\n FFT : fast Fourier transform with Cooley-Tukey algorithm\n\n";
	s[1] = " number of samples\t: %d\n";
	s[2] = " sampling interval\t: %g [s]\n";
	s[3] = " sampling time\t\t: %g [s]\n";
	s[4] = " sampling rate\t\t: %g [Hz]\n";
	s[5] = " nyquist frequency\t: %g [Hz]\n";
	s[6] = " resolution of freq\t: %d\n";
	s[7] = " minimum frequency\t: %g [Hz]\n\n";

	/*  stdout  */
	printf(s[0]);
	printf(s[1], N);
	printf(s[2], DT);
	printf(s[3], N*DT);
	printf(s[4], 1/DT);
	printf(s[5], 0.5/DT);
	printf(s[6], N/2);
	printf(s[7], DF);

	/*  fft.out  */
	fp2 = fopen("fft.out", "w");

	fprintf(fp2, s[0]);
	fprintf(fp2, s[1], N);
	fprintf(fp2, s[2], DT);
	fprintf(fp2, s[3], N*DT);
	fprintf(fp2, s[4], 1/DT);
	fprintf(fp2, s[5], 0.5/DT);
	fprintf(fp2, s[6], N/2);
	fprintf(fp2, s[7], DF);

	fprintf(fp2, " TIME\tINPUT\t\t\t\t\tFFT\n");
	fprintf(fp2, "\t\t real\t\timage\t\t real\t\timage\n");
	for (i = 0; i < N; i++)
		fprintf(fp2, "%5d | %10.3e %10.3e | %10.3e %10.3e\n",
		i, X1[i], Y1[i], X2[i], Y2[i]);

	fclose(fp2);
}

void gnuplot(void)
{
	int i, n = N/2;	float x2, y2;

	/*  gnu.out  */
	fp3 = fopen("gnu.out", "w");
	fprintf(fp3, "%10.3e\t%10.3e\n", 0., 0.);
	for (i = 1; i < n; i++){
		x2 = pow(X2[i], 2);
		y2 = pow(Y2[i], 2);
		fprintf(fp3, "%10.3e\t%10.3e\n", i*DF, 2*sqrt(x2 + y2));
	}
	fclose(fp3);

	/*  gnuplot  */
	if((gp = _popen(GNUPLOT_PATH, "w")) == NULL){
		fprintf(stderr,"can't open gnuplot\n");
		exit(1);
	}
	fprintf(gp, "set logscale xy\n");
	fprintf(gp, "set xlabel \"Hz\"\n");
	fprintf(gp, "set ylabel \"amp\"\n");
	fprintf(gp, "plot \"%s\" with lines linetype 1\n","gnu.out");
	fflush(gp);

	printf("\n press any key to exit...\n");
	getch();

	_pclose(gp);
}