FFT
#define _CRT_SECURE_NO_WARNINGS #define _USE_MATH_DEFINES #include <math.h> #include <stdio.h> #define N 4096 #define FFT 1 #define IFFT -1 void data_in(int, float *, float*); void fft_sub(int, int, float *, float *); float f[N]; void main(void) { int i, n, flag; float fr[N], fi[N]; time_t t1, t2; printf("データ個数を入力 n="); scanf("%d", &n); printf("\n"); data_in(n, fr, fi); time(&t1); flag = FFT; fft_sub(n, flag, fr, fi); time(&t2); printf("time = %f\n", difftime(t2, t1)); printf("データ個数 n=%d \n", n); printf("k f fr fi\n"); for (i = 0; i < n; i++) { printf("%3d %8.8f %8.8f %8.8f\n", i, f[i], fr[i], fi[i]); } printf("%d\n", calCount); } void data_in(int n, float fr[], float fi[]) { //方形波 int i; /* for (i = 0; i < n; i++) { fr[i] = 0.0; fi[i] = 0.0; if (i < n / 4) fr[i] = 1.0; if (i == n / 4) fr[i] = 0.5; if (i == 3 * n / 4) fr[i] = 0.5; if (i > 3 * n / 4) fr[i] = 1.0; if (i > n / 4 && i < 3 * n / 4) fr[i] = 0.0; f[i] = fr[i]; } */ //4_2の波形 } void fft_sub(int n, int flag, float *fr, float *fi) { int i, m, iw, j, k, k0, k2, l, lp, lp2, n2, kk; double arg, darg, c, s, wr, wi, xa, ya; //Bit reversal j = 0; for (i = 0; i < n - 1; i++) { if (i < j) { xa = fr[i]; fr[i] = fr[j]; fr[j] = xa; ya = fi[i]; fi[i] = fi[j]; fi[j] = ya; } n2 = n / 2; while (j >= n2) { j = j - n2; n2 = n2 / 2; } j += n2; } //fft m = 0; n2 = n; while (n2 != 1) { m++; n2 = n2 / 2; } for (l = 1; l <= m; l++)//段 { lp = pow(2.0, (double)l); lp2 = lp / 2; darg = -(double)flag*M_PI / (double)lp2; arg = 0; for (j = 0; j < lp2; j++) { c = cos(arg); s = sin(arg); arg += darg; for (i = j; i < n; i += lp) { iw = i + lp2; /*バタフライ演算部分を自作する*/ } } } if (flag == FFT) { for (i = 0; i < n; i++) { fr[i] /= (double)n; fi[i] /= (double)n; } } }
ソースコード
課題
パソコン実践5.1を演習4.2の問題の波形に対して行った結果を報告せよ.DFTとFFTの値を比較すること.提出レポートには,ソースコード,結果を示すこと.
提出方法
締切日時:6月19日授業時
提出物:課題の問題とプログラムコード,結果を含めたレポートを印刷したもの.