BRL-CAD
ditsplit.c
Go to the documentation of this file.
1 /* D I T S P L I T . 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/ditsplit.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 void
39 irfft(double *x, int n /* length */)
40 {
41  int i, j, k, n1, n2, n4, n8;
42  int i0, i1, i2, i3, i4, i5, i6, i7, i8;
43  int is, id;
44  double t1, t2, t3, t4, t5;
45  double cc1, ss1, cc3, ss3, e, a, a3;
47 
48  irfft_adds = irfft_mults = 0;
49 
50  /* L shaped butterflies */
51  n2 = n << 1;
52  for (k = 2; k < n; k <<= 1) {
53  is = 0;
54  id = n2;
55  n2 = n2 >> 1;
56  n4 = n2 >> 2;
57  n8 = n4 >> 1;
58  e = 2.0*M_PI / n2;
59  l17:
60  for (i = is; i < n; i += id) {
61  i1 = i + 1;
62  i2 = i1 + n4;
63  i3 = i2 + n4;
64  i4 = i3 + n4;
65 
66  t1 = x[i1-1] - x[i3-1];
67  x[i1-1] += x[i3-1];
68  x[i2-1] *= 2.0;
69  x[i3-1] = t1 - 2.0 * x[i4-1];
70  x[i4-1] = t1 + 2.0 * x[i4-1];
71  irfft_adds += 4; irfft_mults += 3;
72 
73  if (n4 == 1)
74  continue;
75  i1 += n8;
76  i2 += n8;
77  i3 += n8;
78  i4 += n8;
79 
80  t1 = (x[i2-1] - x[i1-1]) * M_SQRT1_2;
81  t2 = (x[i4-1] + x[i3-1]) * M_SQRT1_2;
82  x[i1-1] += x[i2-1];
83  x[i2-1] = x[i4-1] - x[i3-1];
84  x[i3-1] = -2.0 * (t2 + t1);
85  x[i4-1] = 2.0 * (-t2 + t1);
86  irfft_adds += 6; irfft_mults += 4;
87 
88  }
89  is = 2 * id - n2;
90  id = 4 * id;
91  if (is < n-1)
92  goto l17;
93 
94  a = e;
95  for (j = 2; j <= n8; j++) {
96  a3 = 3.0 * a;
97  cc1 = cos(a);
98  ss1 = sin(a);
99  cc3 = cos(a3);
100  ss3 = sin(a3);
101  a = j * e;
102  is = 0;
103  id = 2 * n2;
104  l40:
105  for (i = is; i < n; i += id) {
106  i1 = i + j;
107  i2 = i1 + n4;
108  i3 = i2 + n4;
109  i4 = i3 + n4;
110  i5 = i + n4 - j + 2;
111  i6 = i5 + n4;
112  i7 = i6 + n4;
113  i8 = i7 + n4;
114 
115  t1 = x[i1-1] - x[i6-1];
116  x[i1-1] += x[i6-1];
117  t2 = x[i5-1] - x[i2-1];
118  x[i5-1] += x[i2-1];
119  t3 = x[i8-1] + x[i3-1];
120  x[i6-1] = x[i8-1] - x[i3-1];
121  t4 = x[i4-1] + x[i7-1];
122  x[i2-1] = x[i4-1] - x[i7-1];
123  t5 = t1 - t4;
124  t1 += t4;
125  t4 = t2 - t3;
126  t2 += t3;
127  x[i3-1] = t5 * cc1 + t4 * ss1;
128  x[i7-1] = - t4 * cc1 + t5 * ss1;
129  x[i4-1] = t1 * cc3 - t2 * ss3;
130  x[i8-1] = t2 * cc3 + t1 * ss3;
131  irfft_adds += 16; irfft_mults += 8;
132 
133  }
134  is = 2 * id - n2;
135  id = 4 * id;
136  if (is < n-1)
137  goto l40;
138  }
139  }
140 
141  /* Length two butterflies */
142  is = 1;
143  id = 4;
144 l70:
145  for (i0 = is; i0 <= n; i0 += id) {
146  i1 = i0 + 1;
147 
148  t1 = x[i0-1];
149  x[i0-1] = t1 + x[i1-1];
150  x[i1-1] = t1 - x[i1-1];
151  irfft_adds += 2;
152 
153  }
154  is = 2 * id - 1;
155  id = 4 * id;
156  if (is < n)
157  goto l70;
158 
159  /* Digit reverse counter */
160  j = 1;
161  n1 = n - 1;
162  for (i = 1; i <= n1; i++) {
163  if (i < j) {
164  t1 = x[j-1];
165  x[j-1] = x[i-1];
166  x[i-1] = t1;
167  }
168  k = n/2;
169  while (k < j) {
170  j -= k;
171  k /= 2;
172  }
173  j += k;
174  }
175 
176  /* scale result */
177  for (i = 0; i < n; i++)
178  x[i] /= (double)n;
179 }
180 
181 /*
182  * Local Variables:
183  * mode: C
184  * tab-width: 8
185  * indent-tabs-mode: t
186  * c-file-style: "stroustrup"
187  * End:
188  * ex: shiftwidth=4 tabstop=8
189  */
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 irfft(double *x, int n)
Definition: ditsplit.c:39
#define M_SQRT1_2
Definition: fft.h:38