BRL-CAD
ditsplitc.c
Go to the documentation of this file.
1 /* D I T S P L 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/ditsplitc.c
21  *
22  * Split Radix, Decimation in Frequency, Inverse Real-valued FFT.
23  *
24  * Input order:
25  * [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
26  *
27  * Transactions on Acoustics, Speech, and Signal Processing, June 1987.
28  *
29  */
30 
31 #include "common.h"
32 
33 #include <stdio.h>
34 
35 #include "fft.h"
36 
37 
38 /* used by ifftc.c */
40 
41 
42 void
43 ditsplit(int n /* length */, int m /* n = 2^m */)
44 {
45  int i, j, k, n1, n2, n4, n8;
46  int i0, i1, i2, i3, i4, i5, i6, i7, i8;
47  int is, id;
48  double cc1, ss1, cc3, ss3, e, a, a3;
49  irfft_adds = irfft_mults = 0;
50 
51  printf("/*\n"
52  " * BRL-CAD\n"
53  " *\n"
54  " * This file is a generated source file.\n"
55  " * See %s for license and distribution details.\n"
56  " */\n", __FILE__);
57 
58  printf("/*\n"
59  " * Machine-generated Real Split Radix Decimation in Freq Inverse FFT\n"
60  " */\n\n");
61 
62  printf("#include \"fft.h\"\n\n");
63  printf("void\n");
64  printf("irfft%d(register double x[])\n", n);
65  printf("{\n");
66  printf(" register double t1, t2, t3, t4, t5;\n");
67  printf(" register int i;\n");
68 
69  /* L shaped butterflies */
70  printf("\n /* L shaped butterflies */\n");
71  n2 = n << 1;
72  for (k = 1; k < m; k++) {
73  is = 0;
74  id = n2;
75  n2 = n2 >> 1;
76  n4 = n2 >> 2;
77  n8 = n4 >> 1;
78  e = 2.0*M_PI / n2;
79  l17:
80  for (i = is; i < n; i += id) {
81  i1 = i + 1;
82  i2 = i1 + n4;
83  i3 = i2 + n4;
84  i4 = i3 + n4;
85 
86  printf(" t1 = x[%d] - x[%d];\n", i1-1, i3-1);
87  printf(" x[%d] += x[%d];\n", i1-1, i3-1);
88  printf(" x[%d] *= 2.0;\n", i2-1);
89  printf(" x[%d] = t1 - 2.0 * x[%d];\n", i3-1, i4-1);
90  printf(" x[%d] = t1 + 2.0 * x[%d];\n", i4-1, i4-1);
91  irfft_adds += 4; irfft_mults += 3;
92 
93  if (n4 == 1)
94  continue;
95  i1 += n8;
96  i2 += n8;
97  i3 += n8;
98  i4 += n8;
99 
100  printf(" t1 = (x[%d] - x[%d]) * M_SQRT1_2;\n", i2-1, i1-1);
101  printf(" t2 = (x[%d] + x[%d]) * M_SQRT1_2;\n", i4-1, i3-1);
102  printf(" x[%d] += x[%d];\n", i1-1, i2-1);
103  printf(" x[%d] = x[%d] - x[%d];\n", i2-1, i4-1, i3-1);
104  printf(" x[%d] = -2.0 * (t2 + t1);\n", i3-1);
105  printf(" x[%d] = 2.0 * (-t2 + t1);\n", i4-1);
106  irfft_adds += 6; irfft_mults += 4;
107 
108  }
109  is = 2 * id - n2;
110  id = 4 * id;
111  if (is < n-1)
112  goto l17;
113 
114  a = e;
115  for (j = 2; j <= n8; j++) {
116  a3 = 3.0 * a;
117  cc1 = cos(a);
118  ss1 = sin(a);
119  cc3 = cos(a3);
120  ss3 = sin(a3);
121  a = j * e;
122  is = 0;
123  id = 2 * n2;
124  l40:
125  for (i = is; i < n; i += id) {
126  i1 = i + j;
127  i2 = i1 + n4;
128  i3 = i2 + n4;
129  i4 = i3 + n4;
130  i5 = i + n4 - j + 2;
131  i6 = i5 + n4;
132  i7 = i6 + n4;
133  i8 = i7 + n4;
134 
135  printf(" t1 = x[%d] - x[%d];\n", i1-1, i6-1);
136  printf(" x[%d] += x[%d];\n", i1-1, i6-1);
137  printf(" t2 = x[%d] - x[%d];\n", i5-1, i2-1);
138  printf(" x[%d] += x[%d];\n", i5-1, i2-1);
139  printf(" t3 = x[%d] + x[%d];\n", i8-1, i3-1);
140  printf(" x[%d] = x[%d] - x[%d];\n", i6-1, i8-1, i3-1);
141  printf(" t4 = x[%d] + x[%d];\n", i4-1, i7-1);
142  printf(" x[%d] = x[%d] - x[%d];\n", i2-1, i4-1, i7-1);
143  printf(" t5 = t1 - t4;\n");
144  printf(" t1 += t4;\n");
145  printf(" t4 = t2 - t3;\n");
146  printf(" t2 += t3;\n");
147  printf(" x[%d] = t5 * %.24f + t4 * %.24f;\n", i3-1, cc1, ss1);
148  printf(" x[%d] = - t4 * %.24f + t5 * %.24f;\n", i7-1, cc1, ss1);
149  printf(" x[%d] = t1 * %.24f - t2 * %.24f;\n", i4-1, cc3, ss3);
150  printf(" x[%d] = t2 * %.24f + t1 * %.24f;\n", i8-1, cc3, ss3);
151  irfft_adds += 16; irfft_mults += 8;
152 
153  }
154  is = 2 * id - n2;
155  id = 4 * id;
156  if (is < n-1)
157  goto l40;
158  }
159  }
160 
161  /* Length two butterflies */
162  printf("\n /* Length two butterflies */\n");
163  is = 1;
164  id = 4;
165 l70:
166  for (i0 = is; i0 <= n; i0 += id) {
167  i1 = i0 + 1;
168 
169  printf(" t1 = x[%d];\n", i0-1);
170  printf(" x[%d] = t1 + x[%d];\n", i0-1, i1-1);
171  printf(" x[%d] = t1 - x[%d];\n", i1-1, i1-1);
172  irfft_adds += 2;
173 
174  }
175  is = 2 * id - 1;
176  id = 4 * id;
177  if (is < n)
178  goto l70;
179 
180  /* Digit reverse counter */
181  printf("\n /* bit reverse */\n");
182  j = 1;
183  n1 = n - 1;
184  for (i = 1; i <= n1; i++) {
185  if (i < j) {
186  printf(" t1 = x[%d];\n", j-1);
187  printf(" x[%d] = x[%d];\n", j-1, i-1);
188  printf(" x[%d] = t1;\n", i-1);
189  }
190  k = n/2;
191  while (k < j) {
192  j -= k;
193  k /= 2;
194  }
195  j += k;
196  }
197 
198  printf("\n /* scale result */\n");
199  printf(" for (i = 0; i < %d; i++)\n", n);
200  printf("\tx[i] /= %f;\n", (double)n);
201  printf("}\n");
202 }
203 
204 /*
205  * Local Variables:
206  * mode: C
207  * tab-width: 8
208  * indent-tabs-mode: t
209  * c-file-style: "stroustrup"
210  * End:
211  * ex: shiftwidth=4 tabstop=8
212  */
int irfft_adds
Definition: ditsplitc.c:39
#define M_PI
Definition: fft.h:35
Header file for the BRL-CAD common definitions.
int irfft_mults
Definition: ditsplitc.c:39
void ditsplit(int n, int m)
Definition: ditsplitc.c:43