4 /***********************************************************
5 fft.c -- FFT (¹â®FourierÊÑ´¹)
6 ***********************************************************/
8 ´Ø¿ô{\tt fft()}¤Î²¼ÀÁ¤±¤È¤·¤Æ»°³Ñ´Ø¿ôɽ¤òºî¤ë.
10 static void make_sintbl(int n, float sintbl[])
13 double c, s, dc, ds, t;
15 n2 = n / 2; n4 = n / 4; n8 = n / 8;
18 dc = 2 * t * t; ds = sqrt(dc * (2 - dc));
19 t = 2 * dc; c = sintbl[n4] = 1; s = sintbl[0] = 0;
20 for (i = 1; i < n8; i++) {
23 sintbl[i] = s; sintbl[n4 - i] = c;
25 if (n8 != 0) sintbl[n8] = sqrt(0.5);
29 for (i = 1; i < n4; i++)
30 sintbl[i] = sin( (double)i / n * ( 2 * PI ) );
32 for (i = 0; i < n4; i++)
33 sintbl[n2 - i] = sintbl[i];
34 for (i = 0; i < n2 + n4; i++)
35 sintbl[i + n2] = - sintbl[i];
38 ´Ø¿ô{\tt fft()}¤Î²¼ÀÁ¤±¤È¤·¤Æ¥Ó¥Ã¥Èȿžɽ¤òºî¤ë.
40 static void make_bitrev(int n, int bitrev[])
44 n2 = n / 2; i = j = 0;
49 while (k <= j) { j -= k; k /= 2; }
54 ¹â®FourierÊÑ´¹ (Cooley--Tukey¤Î¥¢¥ë¥´¥ê¥º¥à).
55 ɸËÜÅÀ¤Î¿ô {\tt n} ¤Ï2¤ÎÀ°¿ô¾è¤Ë¸Â¤ë.
56 {\tt x[$k$]} ¤¬¼ÂÉô, {\tt y[$k$]} ¤¬µõÉô ($k = 0$, $1$, $2$,
57 \ldots, $|{\tt n}| - 1$).
58 ·ë²Ì¤Ï {\tt x[]}, {\tt y[]} ¤Ë¾å½ñ¤
\ad¤µ¤ì¤ë.
59 ${\tt n} = 0$ ¤Ê¤éɽ¤Î¥á¥â¥ê¤ò²òÊü¤¹¤ë.
60 ${\tt n} < 0$ ¤Ê¤éµÕÊÑ´¹¤ò¹Ô¤¦.
61 Á°²ó¤È°Û¤Ê¤ë $|{\tt n}|$ ¤ÎÃͤǸƤӽФ¹¤È,
62 »°³Ñ´Ø¿ô¤È¥Ó¥Ã¥Èȿž¤Îɽ¤òºî¤ë¤¿¤á¤Ë¿¾¯Í¾Ê¬¤Ë»þ´Ö¤¬¤«¤«¤ë.
63 ¤³¤Îɽ¤Î¤¿¤á¤Îµ
\ad²±Îΰè³ÍÆÀ¤Ë¼ºÇÔ¤¹¤ë¤È1¤òÊÖ¤¹ (Àµ¾ï½ªÎ»»þ
65 ¤³¤ì¤é¤Îɽ¤Îµ
\ad²±Îΰè¤ò²òÊü¤¹¤ë¤Ë¤Ï ${\tt n} = 0$ ¤È¤·¤Æ
66 ¸Æ¤Ó½Ð¤¹ (¤³¤Î¤È¤
\ad¤Ï {\tt x[]}, {\tt y[]} ¤ÎÃͤÏÊѤï¤é¤Ê¤¤).
68 int fft_hontai(int n, float x[], float y[])
70 static int last_n = 0; /* Á°²ó¸Æ½Ð¤·»þ¤Î {\tt n} */
71 static int *bitrev = NULL; /* ¥Ó¥Ã¥Èȿžɽ */
72 static float *sintbl = NULL; /* »°³Ñ´Ø¿ôɽ */
73 int i, j, k, ik, h, d, k2, n4, inverse;
74 float t, s, c, dx, dy;
78 n = -n; inverse = 1; /* µÕÊÑ´¹ */
81 if (n != last_n || n == 0) {
83 if (sintbl != NULL) free(sintbl);
84 if (bitrev != NULL) free(bitrev);
85 if (n == 0) return 0; /* µ
\ad²±Îΰè¤ò²òÊü¤·¤¿ */
86 sintbl = malloc((n + n4) * sizeof(float));
87 bitrev = malloc(n * sizeof(int));
88 if (sintbl == NULL || bitrev == NULL) {
89 fprintf(stderr, "µ
\ad²±ÎΰèÉÔÂ
\ad\n"); return 1;
91 make_sintbl(n, sintbl);
92 make_bitrev(n, bitrev);
94 for (i = 0; i < n; i++) { /* ¥Ó¥Ã¥Èȿž */
97 t = x[i]; x[i] = x[j]; x[j] = t;
98 t = y[i]; y[i] = y[j]; y[j] = t;
101 for (k = 1; k < n; k = k2) { /* ÊÑ´¹ */
102 h = 0; k2 = k + k; d = n / k2;
103 for (j = 0; j < k; j++) {
105 if (inverse) s = - sintbl[h];
107 for (i = j; i < n; i += k2) {
109 dx = s * y[ik] + c * x[ik];
110 dy = c * y[ik] - s * x[ik];
111 x[ik] = x[i] - dx; x[i] += dx;
112 y[ik] = y[i] - dy; y[i] += dy;
117 if (! inverse) /* µÕÊÑ´¹¤Ç¤Ê¤¤¤Ê¤én¤Ç³ä¤ë */
118 for (i = 0; i < n; i++) { x[i] /= n; y[i] /= n; }
119 return 0; /* Àµ¾ï½ªÎ» */
123 int fft( int n, Fukusosuu *in, int disp )
126 static int last_n = 0;
128 if( last_n != abs( n ) )
137 x = calloc( n, sizeof( float ) );
138 y = calloc( n, sizeof( float ) );
142 x[i] = (float)in[i].R;
143 y[i] = (float)in[i].I;
145 fft_hontai( n, x, y );
148 (in+i)->R = (double)x[i];
149 (in+i)->I = (double)y[i];
155 fp = fopen( "outputOfFft", "w" );
157 fprintf( fp, "%f %f\n", in[i].R, in[i].I );
159 system( "vi outputOfFft < /dev/tty > /dev/tty " );
164 int fft( int n, Fukusosuu *in, int disp )
171 x = calloc( m, sizeof( float ) );
172 y = calloc( m, sizeof( float ) );
176 x[i] = (float)in[i].R;
177 y[i] = (float)in[i].I;
184 fp = fopen( "inputOfFft", "w" );
185 fprintf( fp, "m=%d\n", m );
187 fprintf( fp, "%f %f\n", x[i], y[i] );
189 system( "vi inputOfFft < /dev/tty > /dev/tty " );
192 fft_hontai( n, x, y );
197 fp = fopen( "outputOfFft", "w" );
198 fprintf( fp, "m=%d\n", m );
200 fprintf( fp, "%f %f\n", x[i], y[i] );
202 system( "vi outputOfFft < /dev/tty > /dev/tty " );
207 (in+i)->R = (double)x[i];
208 (in+i)->I = (double)y[i];