Branch data Line data Source code
1 : : /*
2 : : * fourier.cpp - fourier transformation class implementation
3 : : *
4 : : * Copyright (C) 2005, 2006, 2009 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <stdlib.h>
31 : : #include <string.h>
32 : : #include <cmath>
33 : :
34 : : #include "consts.h"
35 : : #include "object.h"
36 : : #include "complex.h"
37 : : #include "vector.h"
38 : : #include "fourier.h"
39 : :
40 : : // Helper macro.
41 : : #define Swap(a,b) { nr_double_t t; t = a; a = b; b = t; }
42 : :
43 : : namespace qucs {
44 : :
45 : : using namespace fourier;
46 : :
47 : : /* The function performs a 1-dimensional fast fourier transformation.
48 : : Each data item is meant to be defined in equidistant steps. The
49 : : number of data items needs to be of binary size, e.g. 64, 128. */
50 : 74 : void fourier::_fft_1d (nr_double_t * data, int len, int isign) {
51 : :
52 : : // bit reversal method
53 : : int i, j, m, n;
54 : 74 : n = 2 * len;
55 : 74 : j = 0;
56 [ + + ]: 1258 : for (i = 0; i < n; i += 2) {
57 [ + + ]: 1184 : if (j > i) { // was index already swapped ?
58 : 444 : Swap (data[j], data[i]); // swap real part
59 : 444 : Swap (data[j+1], data[i+1]); // swap imaginary part
60 : : }
61 : 1184 : m = len;
62 [ + + ][ + + ]: 2294 : while (m >= 2 && j >= m) { // calculate next swap index
[ + + ]
63 : 1110 : j -= m;
64 : 1110 : m >>= 1;
65 : : }
66 : 1184 : j += m;
67 : : }
68 : :
69 : : // Danielson-Lanzcos algorithm
70 : : int mmax, istep;
71 : : nr_double_t wt, th, wr, wi, wpr, wpi, ti, tr;
72 : 74 : mmax = 2;
73 [ + + ]: 370 : while (n > mmax) {
74 : 296 : istep = mmax << 1;
75 : 296 : th = isign * (2 * M_PI / mmax);
76 : 296 : wt = sin (0.5 * th);
77 : 296 : wpr = -2.0 * wt * wt;
78 : 296 : wpi = sin (th);
79 : 296 : wr = 1.0;
80 : 296 : wi = 0.0;
81 [ + + ]: 1406 : for (m = 1; m < mmax; m += 2) {
82 [ + + ]: 3478 : for (i = m; i <= n; i += istep) {
83 : 2368 : j = i + mmax;
84 : 2368 : tr = wr * data[j] - wi * data[j-1];
85 : 2368 : ti = wr * data[j-1] + wi * data[j];
86 : 2368 : data[j] = data[i] - tr;
87 : 2368 : data[j-1] = data[i-1] - ti;
88 : 2368 : data[i] += tr;
89 : 2368 : data[i-1] += ti;
90 : : }
91 : 1110 : wr = (wt = wr) * wpr - wi * wpi + wr;
92 : 1110 : wi = wi * wpr + wt * wpi + wi;
93 : : }
94 : 296 : mmax = istep;
95 : : }
96 : 74 : }
97 : :
98 : : /* The function transforms two real vectors using a single fast
99 : : fourier transformation step. The routine works in place. */
100 : 0 : void fourier::_fft_1d_2r (nr_double_t * r1, nr_double_t * r2, int len) {
101 : : int n3, n2, j;
102 : : nr_double_t rep, rem, aip, aim;
103 : 0 : n3 = 1 + (n2 = len + len);
104 : :
105 : : // put the two real vectors into one complex vector
106 [ # # ]: 0 : for (j = 1; j <= n2; j += 2) {
107 : 0 : r1[j] = r2[j-1];
108 : : }
109 : :
110 : : // transform the complex vector
111 : 0 : _fft_1d (r1, len, 1);
112 : :
113 : : // separate the two transforms
114 : 0 : r2[0] = r1[1];
115 : 0 : r1[1] = r2[1] = 0.0;
116 [ # # ]: 0 : for (j = 2; j <= len; j += 2) {
117 : : // use symmetries to separate the two transforms
118 : 0 : rep = 0.5 * (r1[j] + r1[n2-j]);
119 : 0 : rem = 0.5 * (r1[j] - r1[n2-j]);
120 : 0 : aip = 0.5 * (r1[j+1] + r1[n3-j]);
121 : 0 : aim = 0.5 * (r1[j+1] - r1[n3-j]);
122 : : // ship them out in two complex vectors
123 : 0 : r1[j+1] = aim;
124 : 0 : r2[j+1] = -rem;
125 : 0 : r1[j] = r1[n2-j] = rep;
126 : 0 : r2[j] = r2[n2-j] = aip;
127 : 0 : r1[n3-j] = -aim;
128 : 0 : r2[n3-j] = rem;
129 : : }
130 : 0 : }
131 : :
132 : : /* The following function transforms two vectors yielding real valued
133 : : vectors using a single inverse fast fourier transformation step.
134 : : The routine works in place as well. */
135 : 0 : void fourier::_ifft_1d_2r (nr_double_t * r1, nr_double_t * r2, int len) {
136 : : nr_double_t r, i;
137 : 0 : int j, jj, nn = len + len;
138 : :
139 : : // put the two complex vectors into one complex vector
140 [ # # ]: 0 : for (j = 0, jj = 0; j < nn; j += 2) {
141 : 0 : r = r1[j] - r2[j+1];
142 : 0 : i = r1[j+1] + r2[j];
143 : 0 : r1[jj++] = r;
144 : 0 : r1[jj++] = i;
145 : : }
146 : :
147 : : // transform the complex vector
148 : 0 : _fft_1d (r1, len, -1);
149 : :
150 : : // split the transforms into two real vectors
151 [ # # ]: 0 : for (j = 0; j < nn; j += 2) {
152 : 0 : r2[j] = r1[j+1];
153 : 0 : r1[j+1] = r2[j+1] = 0.0;
154 : : }
155 : 0 : }
156 : :
157 : : /* This function performs a 1-dimensional fast fourier transformation
158 : : on the given vector 'var'. If 'sign' is -1 the inverse fft is
159 : : computed, if +1 the fft itself is computed. It returns a vector of
160 : : binary size (as necessary for a fft). */
161 : 0 : vector fourier::fft_1d (vector var, int isign) {
162 : 0 : int i, n, len = var.getSize ();
163 : :
164 : : // compute necessary binary data array size
165 : 0 : int size = 2;
166 [ # # ]: 0 : while (size < len) size <<= 1;
167 : :
168 : : // put data items (temporary array) in place
169 : : nr_double_t * data;
170 : 0 : data = (nr_double_t *) calloc (2 * size * sizeof (nr_double_t), 1);
171 [ # # ]: 0 : for (n = i = 0; i < len; i++, n += 2) {
172 : 0 : data[n] = real (var (i)); data[n+1] = imag (var (i));
173 : : }
174 : :
175 : : // run 1-dimensional fft
176 : 0 : _fft_1d (data, size, isign);
177 : :
178 : : // store transformed data items in result vector
179 : 0 : vector res = vector (size);
180 [ # # ]: 0 : for (n = i = 0; i < size; i++, n += 2) {
181 : 0 : res (i) = nr_complex_t (data[n], data[n+1]);
182 [ # # ]: 0 : if (isign < 0) res (i) /= size;
183 : : }
184 : :
185 : : // free temporary data array
186 : 0 : free (data);
187 : 0 : return res;
188 : : }
189 : :
190 : : /* The function performs a 1-dimensional discrete fourier
191 : : transformation. Each data item is meant to be defined in
192 : : equidistant steps. */
193 : 0 : void fourier::_dft_1d (nr_double_t * data, int len, int isign) {
194 : 0 : int k, n, size = 2 * len * sizeof (nr_double_t);
195 : 0 : nr_double_t * res = (nr_double_t *) calloc (size, 1);
196 : : nr_double_t th, c, s;
197 [ # # ]: 0 : for (n = 0; n < 2 * len; n += 2) {
198 : 0 : th = n * M_PI / 2 / len;
199 [ # # ]: 0 : for (k = 0; k < 2 * len; k += 2) {
200 : 0 : c = cos (k * th);
201 : 0 : s = isign * sin (k * th);
202 : 0 : res[n] += data[k] * c + data[k+1] * s;
203 : 0 : res[n+1] += data[k+1] * c - data[k] * s;
204 : : }
205 : : }
206 : 0 : memcpy (data, res, size);
207 : 0 : free (res);
208 : 0 : }
209 : :
210 : : /* The function performs a 1-dimensional discrete fourier
211 : : transformation on the given vector 'var'. If 'sign' is -1 the
212 : : inverse dft is computed, if +1 the dft itself is computed. */
213 : 3 : vector fourier::dft_1d (vector var, int isign) {
214 : 3 : int k, n, len = var.getSize ();
215 : 3 : vector res = vector (len);
216 [ + + ]: 861 : for (n = 0; n < len; n++) {
217 : 858 : nr_double_t th = - isign * 2 * M_PI * n / len;
218 : 858 : nr_complex_t val = 0;
219 [ + + ]: 267596 : for (k = 0; k < len; k++)
220 [ + - ]: 266738 : val += var (k) * std::polar (1.0, th * k);
221 [ - + ][ - + ]: 858 : res (n) = isign < 0 ? val / (nr_double_t) len : val;
222 : : }
223 : 3 : return res;
224 : : }
225 : :
226 : : // Helper functions.
227 : 0 : vector fourier::ifft_1d (vector var) {
228 [ # # ]: 0 : return fft_1d (var, -1);
229 : : }
230 : 0 : vector fourier::idft_1d (vector var) {
231 [ # # ]: 0 : return dft_1d (var, -1);
232 : : }
233 : 0 : void fourier::_ifft_1d (nr_double_t * data, int len) {
234 : 0 : _fft_1d (data, len, -1);
235 : 0 : }
236 : 0 : void fourier::_idft_1d (nr_double_t * data, int len) {
237 : 0 : _dft_1d (data, len, -1);
238 : 0 : }
239 : :
240 : : /* The function performs a n-dimensional fast fourier transformation.
241 : : Each data item is meant to be defined in equidistant steps. The
242 : : number of data items needs to be of binary size, e.g. 64, 128 for
243 : : each dimension. */
244 : 0 : void fourier::_fft_nd (nr_double_t * data, int len[], int nd, int isign) {
245 : :
246 : : int i, i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
247 : : int ibit, k1, k2, n, np, nr, nt;
248 : :
249 : : // compute total number of complex values
250 [ # # ]: 0 : for (nt = 1, i = 0; i < nd; i++) nt *= len[i];
251 : :
252 : : // main loop over the dimensions
253 [ # # ]: 0 : for (np = 1, i = nd - 1; i >= 0; i--) {
254 : 0 : n = len[i];
255 : 0 : nr = nt / (n * np);
256 : 0 : ip1 = np << 1;
257 : 0 : ip2 = ip1 * n;
258 : 0 : ip3 = ip2 * nr;
259 : :
260 : : // bit-reversal method
261 [ # # ]: 0 : for (i2rev = 1, i2 = 1; i2 <= ip2; i2 += ip1) {
262 [ # # ]: 0 : if (i2 < i2rev) {
263 [ # # ]: 0 : for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
264 [ # # ]: 0 : for (i3 = i1; i3 <= ip3; i3 += ip2) {
265 : 0 : i3rev = i2rev + i3 - i2;
266 : 0 : Swap (data[i3-1], data[i3rev-1]);
267 : 0 : Swap (data[i3], data[i3rev]);
268 : : }
269 : : }
270 : : }
271 : 0 : ibit = ip2 >> 1;
272 [ # # ][ # # ]: 0 : while (ibit >= ip1 && i2rev > ibit) {
[ # # ]
273 : 0 : i2rev -= ibit;
274 : 0 : ibit >>= 1;
275 : : }
276 : 0 : i2rev += ibit;
277 : : }
278 : :
279 : : // Danielson-Lanzcos algorithm
280 : : nr_double_t ti, tr, wt, th, wr, wi, wpi, wpr;
281 : 0 : ifp1 = ip1;
282 [ # # ]: 0 : while (ifp1 < ip2) {
283 : 0 : ifp2 = ifp1 << 1;
284 : 0 : th = isign * 2 * M_PI / (ifp2 / ip1);
285 : 0 : wt = sin (0.5 * th);
286 : 0 : wpr = -2.0 * wt * wt;
287 : 0 : wpi = sin (th);
288 : 0 : wr = 1.0;
289 : 0 : wi = 0.0;
290 [ # # ]: 0 : for (i3 = 1; i3 <= ifp1; i3 += ip1) {
291 [ # # ]: 0 : for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
292 [ # # ]: 0 : for (i2 = i1; i2 <= ip3; i2 += ifp2) {
293 : 0 : k1 = i2;
294 : 0 : k2 = k1 + ifp1;
295 : 0 : tr = wr * data[k2-1] - wi * data[k2];
296 : 0 : ti = wr * data[k2] + wi * data[k2-1];
297 : 0 : data[k2-1] = data[k1-1] - tr;
298 : 0 : data[k2] = data[k1] - ti;
299 : 0 : data[k1-1] += tr;
300 : 0 : data[k1] += ti;
301 : : }
302 : : }
303 : 0 : wr = (wt = wr) * wpr - wi * wpi + wr;
304 : 0 : wi = wi * wpr + wt * wpi + wi;
305 : : }
306 : 0 : ifp1 = ifp2;
307 : : }
308 : 0 : np *= n;
309 : : }
310 : 0 : }
311 : :
312 : : // Helper functions.
313 : 0 : void fourier::_ifft_nd (nr_double_t * data, int len[], int nd) {
314 : 0 : _fft_nd (data, len, nd, -1);
315 : 0 : }
316 : :
317 : : // Shuffles values of FFT around.
318 : 0 : vector fourier::fftshift (vector var) {
319 : 0 : int i, n, len = var.getSize ();
320 : 0 : vector res = vector (len);
321 : 0 : n = len / 2;
322 [ # # ]: 0 : for (i = 0; i < len / 2; i++) {
323 : 0 : res (i) = var (n + i);
324 : 0 : res (i + n) = var (i);
325 : : }
326 : 0 : return res;
327 : : }
328 : :
329 : : } // namespace qucs
|
