BRL-CAD
splitdit.c
Go to the documentation of this file.
1 /* S P L I T D 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/splitdit.c
21  *
22  * Real valued, split-radix, decimation in time FFT.
23  *
24  * Data comes out as:
25  *
26  * [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
27  *
28  */
29 
30 #include "fft.h"
31 
32 
34 
35 
36 void
37 rfft(double *X, int N)
38 {
39  int i0, i1, i2, i3;
40  int a0, a1, a2, a3, b0, b1, b2, b3;
41  int s, d;
42  double t0, t1, t2, r1, a, aa3, e, c2, c3, d2, d3;
43  double cc1, ss1, cc3, ss3, xt;
44  int i, j, k, ni;
45  int n2, n4;
46  rfft_adds = rfft_mults = 0;
47 
48  /* bit reverse counter */
49  j = 1;
50  ni = N - 1;
51  for ( i = 1; i <= ni; i++ ) {
52  if ( i < j ) {
53  xt = X[j-1];
54  X[j-1] = X[i-1];
55  X[i-1] = xt;
56  }
57  k = N/2;
58  while ( k < j ) {
59  j -= k;
60  k /= 2;
61  }
62  j += k;
63  }
64 
65  /* length two transforms */
66  for ( s = 1, d = 4; s < N; s = 2*d-1, d *= 4 ) {
67  for ( i0 = s; i0 <= N; i0 += d ) {
68  i1 = i0 + 1;
69  r1 = X[i0-1];
70  X[i0-1] = r1 + X[i1-1];
71  X[i1-1] = r1 - X[i1-1];
72  rfft_adds += 2;
73  }
74  }
75 
76  /* other butterflies */
77  n2 = 2;
78 /* for ( k = 2; k <= M; k++ ) {*/
79  for ( k = 4; k <= N; k <<= 1 ) {
80  n2 *= 2;
81  n4 = n2/4;
82  /* without mult */
83  for ( s = 1, d = 2*n2; s < N; s = 2*d-n2+1, d *= 4 ) {
84  for ( i0 = s; i0 < N; i0 += d ) {
85  i1 = i0 + n4;
86  i2 = i1 + n4;
87  i3 = i2 + n4;
88  t0 = X[i2-1] + X[i3-1];
89  X[i3-1] = X[i2-1] - X[i3-1];
90  X[i2-1] = X[i0-1] - t0;
91  X[i0-1] += t0;
92  rfft_adds += 4;
93  }
94  }
95  if ( n4 < 2 ) continue;
96  /* with 2 real mult */
97  for ( s = n4/2+1, d = 2*n2; s < N; s = 2*d-n2+n4/2+1, d *= 4 ) {
98  for ( i0 = s; i0 < N; i0 += d ) {
99  i1 = i0 + n4;
100  i2 = i1 + n4;
101  i3 = i2 + n4;
102  t1 = (X[i2-1]-X[i3-1])*M_SQRT1_2;
103  t2 = (X[i2-1]+X[i3-1])*M_SQRT1_2;
104  X[i2-1] = t2 - X[i1-1];
105  X[i3-1] = t2 + X[i1-1];
106  X[i1-1] = X[i0-1] - t1;
107  X[i0-1] += t1;
108  rfft_mults += 2; rfft_adds += 6;
109  }
110  }
111  e = 2.0*M_PI/n2;
112  a = e;
113  if ( n4 < 4 ) continue;
114  for ( j = 2; j <= n4/2; j++ ) {
115  aa3 = 3*a;
116  cc1 = cos(a);
117  ss1 = sin(a);
118  cc3 = cos(aa3);
119  ss3 = sin(aa3);
120  a = j * e;
121  /* with 6 real mult */
122  for ( s = j, d = 2*n2; s < N; s = 2*d-n2+j, d *= 4 ) {
123  for ( a0 = s; a0 < N; a0 += d ) {
124  b1 = a0 + n4;
125  a1 = b1-j-j+2;
126  b0 = a1 + n4;
127  a2 = b1 + n4;
128  a3 = a2 + n4;
129  b2 = b0 + n4;
130  b3 = b2 + n4;
131  c2 = X[a2-1]*cc1 - X[b2-1]*ss1;
132  d2 = -(X[a2-1]*ss1 + X[b2-1]*cc1);
133  c3 = X[a3-1]*cc3 - X[b3-1]*ss3;
134  d3 = -(X[a3-1]*ss3 + X[b3-1]*cc3);
135  rfft_mults += 8; rfft_adds += 4;
136  t1 = c2 + c3;
137  c3 = c2 - c3;
138  t2 = d2 - d3;
139  d3 = d2 + d3;
140  X[a2-1] = -X[b0-1] - d3;
141  X[b2-1] = -X[b1-1] + c3;
142  X[a3-1] = X[b1-1] + c3;
143  X[b3-1] = X[b0-1] - d3;
144  X[b1-1] = X[a1-1] + t2;
145  X[b0-1] = X[a0-1] - t1;
146  X[a0-1] += t1;
147  X[a1-1] -= t2;
148  rfft_adds += 12;
149  }
150  }
151  }
152  }
153 
154  /*
155  * For some reason the Imag part is coming out with the wrong
156  * sign, so we reverse it here! We need to figure this out!
157  */
158  for ( i = N/2+1; i < N; i++ )
159  X[i] = -X[i];
160 }
161 
162 /*
163  * Local Variables:
164  * mode: C
165  * tab-width: 8
166  * indent-tabs-mode: t
167  * c-file-style: "stroustrup"
168  * End:
169  * ex: shiftwidth=4 tabstop=8
170  */
if lu s
Definition: nmg_mod.c:3860
#define N
Definition: randmt.c:39
#define M_PI
Definition: fft.h:35
Definition: color.c:49
int rfft_adds
Definition: splitdit.c:33
void rfft(double *X, int N)
Definition: splitdit.c:37
int rfft_mults
Definition: splitdit.c:33
#define M_SQRT1_2
Definition: fft.h:38