BRL-CAD
splitditc.c
Go to the documentation of this file.
1 /* S P L I T D I T C . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file libfft/splitditc.c
21  *
22  * Real valued, split-radix, decimation in time FFT code generator.
23  *
24  */
25 
26 #include "common.h"
27 
28 #include <stdio.h>
29 
30 #include "fft.h"
31 
32 
33 /* used by fftc.c */
35 
36 
37 void
38 splitdit(int N, int M)
39 {
40  int i0, i1, i2, i3;
41  int a0, a1, a2, a3, b0, b1, b2, b3;
42  int s, d;
43  double a, aa3, e;
44  double cc1, ss1, cc3, ss3;
45  int i, j, k, ni;
46  int n2, n4;
47  rfft_adds = rfft_mults = 0;
48 
49  printf("/*\n"
50  " * BRL-CAD\n"
51  " *\n"
52  " * This file is a generated source file.\n"
53  " * See %s for license and distribution details.\n"
54  " */\n", __FILE__);
55 
56  printf("/*\n"
57  " * Machine-generated Real Split Radix Decimation in Time FFT\n"
58  " */\n\n");
59  printf("#include \"fft.h\"\n\n");
60  printf("void\n");
61  printf("rfft%d(register double X[])\n", N);
62  printf("{\n");
63  printf(" register double t0, t1, c3, d3, c2, d2;\n");
64  printf(" register int i;\n");
65 
66  /* bit reverse counter */
67  printf("\n /* bit reverse */\n");
68  j = 1;
69  ni = N - 1;
70  for (i = 1; i <= ni; i++) {
71  if (i < j) {
72  printf(" t0 = X[%d];\n", j-1);
73  printf(" X[%d] = X[%d];\n", j-1, i-1);
74  printf(" X[%d] = t0;\n", i-1);
75  }
76  k = N/2;
77  while (k < j) {
78  j -= k;
79  k /= 2;
80  }
81  j += k;
82  }
83 
84  /* length two transforms */
85  printf("\n /* length two xforms */\n");
86  for (s = 1, d = 4; s < N; s = 2*d-1, d *= 4) {
87  for (i0 = s; i0 <= N; i0 += d) {
88  i1 = i0 + 1;
89  printf(" t0 = X[%d];\n", i0-1);
90  printf(" X[%d] += X[%d];\n", i0-1, i1-1);
91  printf(" X[%d] = t0 - X[%d];\n", i1-1, i1-1);
92  rfft_adds += 2;
93  }
94  }
95 
96  /* other butterflies */
97  printf("\n /* other butterflies */\n");
98  n2 = 2;
99  for (k = 2; k <= M; k++) {
100  n2 *= 2;
101  n4 = n2/4;
102 
103  /* without mult */
104  for (s = 1, d = 2*n2; s < N; s = 2*d-n2+1, d *= 4) {
105  for (i0 = s; i0 < N; i0 += d) {
106  i1 = i0 + n4;
107  i2 = i1 + n4;
108  i3 = i2 + n4;
109  printf(" t0 = X[%d] + X[%d];\n", i3-1, i2-1);
110  printf(" X[%d] = X[%d] - X[%d];\n", i3-1, i2-1, i3-1);
111  printf(" X[%d] = X[%d] - t0;\n", i2-1, i0-1);
112  printf(" X[%d] += t0;\n", i0-1);
113  rfft_adds += 4;
114  }
115  }
116  if (n4 < 2) continue;
117  /* with 2 real mult */
118  for (s = n4/2+1, d = 2*n2; s < N; s = 2*d-n2+n4/2+1, d *= 4) {
119  for (i0 = s; i0 < N; i0 += d) {
120  i1 = i0 + n4;
121  i2 = i1 + n4;
122  i3 = i2 + n4;
123  printf(" t0 = (X[%d]-X[%d])*M_SQRT1_2;\n", i2-1, i3-1);
124  printf(" t1 = (X[%d]+X[%d])*M_SQRT1_2;\n", i2-1, i3-1);
125  printf(" X[%d] = t1 - X[%d];\n", i2-1, i1-1);
126  printf(" X[%d] = t1 + X[%d];\n", i3-1, i1-1);
127  printf(" X[%d] = X[%d] - t0;\n", i1-1, i0-1);
128  printf(" X[%d] += t0;\n", i0-1);
129  rfft_mults += 2; rfft_adds += 6;
130  }
131  }
132  e = 2.0*M_PI/n2;
133  a = e;
134  if (n4 < 4) continue;
135  for (j = 2; j <= n4/2; j++) {
136  aa3 = 3*a;
137  cc1 = cos(a);
138  ss1 = sin(a);
139  cc3 = cos(aa3);
140  ss3 = sin(aa3);
141  a = j * e;
142  /* with 6 real mult */
143  for (s = j, d = 2*n2; s < N; s = 2*d-n2+j, d *= 4) {
144  for (a0 = s; a0 < N; a0 += d) {
145  b1 = a0 + n4;
146  a1 = b1-j-j+2;
147  b0 = a1 + n4;
148  a2 = b1 + n4;
149  a3 = a2 + n4;
150  b2 = b0 + n4;
151  b3 = b2 + n4;
152  printf(" c2 = X[%d]*%.24g - X[%d]*%.24g;\n", a2-1, cc1, b2-1, ss1);
153  printf(" d2 = -(X[%d]*%.24g + X[%d]*%.24g);\n", a2-1, ss1, b2-1, cc1);
154  printf(" c3 = X[%d]*%.24g - X[%d]*%.24g;\n", a3-1, cc3, b3-1, ss3);
155  printf(" d3 = -(X[%d]*%.24g + X[%d]*%.24g);\n", a3-1, ss3, b3-1, cc3);
156  rfft_mults += 8; rfft_adds += 4;
157  printf(" t0 = c2 + c3;\n");
158  printf(" c3 = c2 - c3;\n");
159  printf(" t1 = d2 - d3;\n");
160  printf(" d3 += d2;\n");
161  printf(" X[%d] = -X[%d] - d3;\n", a2-1, b0-1);
162  printf(" X[%d] = -X[%d] + c3;\n", b2-1, b1-1);
163  printf(" X[%d] = X[%d] + c3;\n", a3-1, b1-1);
164  printf(" X[%d] = X[%d] - d3;\n", b3-1, b0-1);
165  printf(" X[%d] = X[%d] + t1;\n", b1-1, a1-1);
166  printf(" X[%d] = X[%d] - t0;\n", b0-1, a0-1);
167  printf(" X[%d] += t0;\n", a0-1);
168  printf(" X[%d] -= t1;\n", a1-1);
169  rfft_adds += 12;
170  }
171  }
172  }
173  }
174 
175  /* For some reason the Imag part is coming out with the wrong
176  * sign, so we reverse it here! We need to figure this out!
177  */
178  printf("\n /* reverse Imag part! */\n");
179  printf(" for (i = %d/2+1; i < %d; i++)\n", N, N);
180  printf("\tX[i] = -X[i];\n");
181  printf("}\n");
182 }
183 
184 /*
185  * Local Variables:
186  * mode: C
187  * tab-width: 8
188  * indent-tabs-mode: t
189  * c-file-style: "stroustrup"
190  * End:
191  * ex: shiftwidth=4 tabstop=8
192  */
if lu s
Definition: nmg_mod.c:3860
#define N
Definition: randmt.c:39
#define M_PI
Definition: fft.h:35
Header file for the BRL-CAD common definitions.
int rfft_mults
Definition: splitditc.c:34
int rfft_adds
Definition: splitditc.c:34
void splitdit(int N, int M)
Definition: splitditc.c:38
#define M
Definition: msr.c:52