aboutsummaryrefslogtreecommitdiffstatshomepage
path: root/linden/indra/newview/llfft.cpp
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--linden/indra/newview/llfft.cpp577
1 files changed, 0 insertions, 577 deletions
diff --git a/linden/indra/newview/llfft.cpp b/linden/indra/newview/llfft.cpp
deleted file mode 100644
index 2b80ad9..0000000
--- a/linden/indra/newview/llfft.cpp
+++ /dev/null
@@ -1,577 +0,0 @@
1/**
2 * @file llfft.cpp
3 * @brief FFT function implementations
4 *
5 * Copyright (c) 2003-2007, Linden Research, Inc.
6 *
7 * Second Life Viewer Source Code
8 * The source code in this file ("Source Code") is provided by Linden Lab
9 * to you under the terms of the GNU General Public License, version 2.0
10 * ("GPL"), unless you have obtained a separate licensing agreement
11 * ("Other License"), formally executed by you and Linden Lab. Terms of
12 * the GPL can be found in doc/GPL-license.txt in this distribution, or
13 * online at http://secondlife.com/developers/opensource/gplv2
14 *
15 * There are special exceptions to the terms and conditions of the GPL as
16 * it is applied to this Source Code. View the full text of the exception
17 * in the file doc/FLOSS-exception.txt in this software distribution, or
18 * online at http://secondlife.com/developers/opensource/flossexception
19 *
20 * By copying, modifying or distributing this software, you acknowledge
21 * that you have read and understood your obligations described above,
22 * and agree to abide by those obligations.
23 *
24 * ALL LINDEN LAB SOURCE CODE IS PROVIDED "AS IS." LINDEN LAB MAKES NO
25 * WARRANTIES, EXPRESS, IMPLIED OR OTHERWISE, REGARDING ITS ACCURACY,
26 * COMPLETENESS OR PERFORMANCE.
27 */
28
29/*
30 * Fast Fourier Transform
31 *
32 */
33#include "llviewerprecompiledheaders.h"
34
35#include "llfft.h"
36#include "llerror.h"
37
38/*
39**
40*********************************************************************
41** Forward and inverse discrete Fourier transforms on complex data **
42*********************************************************************
43**
44**
45** forward_fft(array, rows, cols)
46** COMPLEX *array;
47** S32 rows, cols;
48**
49** inverse_fft(array, rows, cols)
50** COMPLEX *array;
51** S32 rows, cols;
52**
53** These entry points compute the forward and inverse DFT's, respectively,
54** of a single-precision COMPLEX array.
55**
56** The result is a COMPLEX array of the same size, returned in
57** the same space as the input array. That is, the original array is
58** overwritten and destroyed.
59**
60** Rows and columns must each be an integral power of 2.
61**
62** These routines return integer value -1 if an error was detected,
63** 0 otherwise
64**
65** This implementation of the DFT uses the transform pair defined as follows.
66**
67** Let there be two COMPLEX arrays each with n rows and m columns
68** Index them as
69** f(x,y): 0 <= x <= m - 1, 0 <= y <= n - 1
70** F(u,v): -m/2 <= u <= m/2 - 1, -n/2 <= v <= n/2 - 1
71**
72** Then the forward and inverse transforms are related as
73**
74** Forward:
75**
76** F(u,v) = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1}
77** f(x,y) \exp{-2\pi i (ux/m + vy/n)}
78**
79**
80** Inverse:
81**
82** f(x,y) = 1/(mn) \sum_{u=-m/2}^{m/2-1} \sum_{v=-n/2}^{n/2-1}
83** F(u,v) \exp{2\pi i (ux/m + vy/n)}
84**
85** Therefore, the transforms have these properties:
86** 1. \sum_x \sum_y f(x,y) = F(0,0)
87** 2. m n \sum_x \sum_y |f(x,y)|^2 = \sum_u \sum_v |F(u,v)|^2
88**
89*/
90
91
92//DPCOMPLEX *stageBuff; /* buffer to hold a row or column at a time */
93//COMPLEX *bigBuff; /* a pointer to the input array */
94
95
96/*
97 * These macros move complex data between bigBuff and
98 * stageBuff
99 */
100
101inline void LoadRow(DPCOMPLEX* stageBuff, COMPLEX* bigBuff, U32 row, U32 cols)
102{
103 for (U32 j = row*cols, k = 0 ; k < cols ; j++, k++)
104 {
105 stageBuff[k].re = bigBuff[j].re;
106 stageBuff[k].im = bigBuff[j].im;
107 }
108}
109
110inline void StoreRow(DPCOMPLEX* stageBuff, COMPLEX* bigBuff, U32 row, U32 cols)
111{
112 for (U32 j = row*cols, k = 0 ; k < cols ; j++, k++)
113 {
114 bigBuff[j].re = (F32)stageBuff[k].re;
115 bigBuff[j].im = (F32)stageBuff[k].im;
116 }
117}
118
119inline void LoadCol(DPCOMPLEX* stageBuff, COMPLEX* bigBuff, U32 col, U32 rows, U32 cols)
120{
121 for (U32 j = col,k = 0 ; k < rows ; j+=cols, k++)
122 {
123 stageBuff[k].re = bigBuff[j].re;
124 stageBuff[k].im = bigBuff[j].im;
125 }
126}
127
128inline void StoreCol(DPCOMPLEX* stageBuff, COMPLEX* bigBuff, U32 col, U32 rows, U32 cols)
129{
130 for (U32 j = col,k = 0 ; k < rows ; j+=cols, k++)
131 {
132 bigBuff[j].re = (F32)stageBuff[k].re;
133 bigBuff[j].im = (F32)stageBuff[k].im;
134 }
135}
136
137
138/* do something with an error message */
139inline void handle_error(S8* msg)
140{
141 llerrs << msg << llendl;
142}
143
144
145/*
146** compute DFT: forward if direction==0, inverse if direction==1
147** array must be COMPLEX
148*/
149BOOL fft(const LLFFTPlan& plan, COMPLEX *array, S32 rows, S32 cols, S32 direction)
150{
151 S32 i;
152
153 if (!plan.valid() || plan.rows() != rows || plan.cols() != cols)
154 return FALSE;
155
156 /* compute transform row by row */
157
158 if(cols>1)
159 {
160 for(i=0;i<rows;i++)
161 {
162 LoadRow(plan.buffer(), array, i, cols);
163 FFT842(direction, cols, plan.buffer());
164 StoreRow(plan.buffer(), array, i, cols);
165 }
166 }
167
168 if(rows<2) /* done */
169 {
170 //freeBuffer();
171 return TRUE;
172 }
173
174
175 /* compute transform column by column */
176
177 for(i=0;i<cols;i++)
178 {
179 LoadCol(plan.buffer(), array, i, rows, cols);
180 FFT842(direction, rows, plan.buffer());
181 StoreCol(plan.buffer(), array, i, rows, cols);
182 }
183
184 //freeBuffer();
185
186 return TRUE;
187}
188
189
190
191/*
192** FFT842
193** This routine replaces the input DPCOMPLEX vector by its
194** finite discrete complex fourier transform if in==0.
195** It replaces the input DPCOMPLEX vector by its
196** finite discrete complex inverse fourier transform if in==1.
197**
198** in - FORWARD or INVERSE
199** n - length of vector
200** b - input vector
201**
202** It performs as many base 8 iterations as possible and
203** then finishes with a base 4 iteration or a base 2
204** iteration if needed.
205**
206** Ported from the FORTRAN code in Programming for Digital Signal Processing,
207** IEEE Press 1979, Section 1, by G. D. Bergland and M. T. Dolan
208**
209*/
210void FFT842(S32 in, S32 n, DPCOMPLEX *b)
211{
212 F64 fn, r, fi;
213
214 S32 L[16],L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15;
215 S32 j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14;
216 S32 i, j, ij, ji, ij1, ji1;
217 S32 n2pow, n8pow, nthpo, ipass, nxtlt, lengt;
218
219 n2pow = fastlog2(n);
220 nthpo = n;
221 fn = 1.0 / (F64)nthpo;
222
223
224 if(in==FORWARD)
225 {
226 /* take conjugate */
227 for(i=0;i<n;i++)
228 {
229 b[i].im *= -1.0;
230 }
231 }
232
233 if(in==INVERSE)
234 {
235 /* scramble inputs */
236 for(i=0,j=n/2;j<n;i++,j++)
237 {
238 r = b[j].re;
239 fi = b[j].im;
240 b[j].re = b[i].re;
241 b[j].im = b[i].im;
242 b[i].re = r;
243 b[i].im = fi;
244 }
245 }
246
247 n8pow = n2pow/3;
248
249 if(n8pow)
250 {
251 /* radix 8 iterations */
252 for(ipass=1;ipass<=n8pow;ipass++)
253 {
254 nxtlt = 0x1 << (n2pow - 3*ipass);
255 lengt = 8*nxtlt;
256 R8TX(nxtlt,nthpo,lengt,
257 b,b+nxtlt,b+2*nxtlt,
258 b+3*nxtlt,b+4*nxtlt,b+5*nxtlt,
259 b+6*nxtlt,b+7*nxtlt);
260 }
261 }
262
263 if(n2pow%3 == 1)
264 {
265 /* radix 2 iteration needed */
266 R2TX(nthpo,b,b+1);
267 }
268
269
270 if(n2pow%3 == 2)
271 {
272 /* radix 4 iteration needed */
273 R4TX(nthpo,b,b+1,b+2,b+3);
274 }
275
276
277
278 for(j=1;j<=15;j++)
279 {
280 L[j] = 1;
281 if(j-n2pow <= 0) L[j] = 0x1 << (n2pow + 1 - j);
282 }
283 L15=L[1];L14=L[2];L13=L[3];L12=L[4];L11=L[5];L10=L[6];L9=L[7];
284 L8=L[8];L7=L[9];L6=L[10];L5=L[11];L4=L[12];L3=L[13];L2=L[14];L1=L[15];
285
286 ij = 1;
287
288 for(j1=1;j1<=L1;j1++)
289 for(j2=j1;j2<=L2;j2+=L1)
290 for(j3=j2;j3<=L3;j3+=L2)
291 for(j4=j3;j4<=L4;j4+=L3)
292 for(j5=j4;j5<=L5;j5+=L4)
293 for(j6=j5;j6<=L6;j6+=L5)
294 for(j7=j6;j7<=L7;j7+=L6)
295 for(j8=j7;j8<=L8;j8+=L7)
296 for(j9=j8;j9<=L9;j9+=L8)
297 for(j10=j9;j10<=L10;j10+=L9)
298 for(j11=j10;j11<=L11;j11+=L10)
299 for(j12=j11;j12<=L12;j12+=L11)
300 for(j13=j12;j13<=L13;j13+=L12)
301 for(j14=j13;j14<=L14;j14+=L13)
302 for(ji=j14;ji<=L15;ji+=L14)
303 {
304 ij1 = ij-1;
305 ji1 = ji-1;
306
307 if(ij-ji<0)
308 {
309 r = b[ij1].re;
310 b[ij1].re = b[ji1].re;
311 b[ji1].re = r;
312 fi = b[ij1].im;
313 b[ij1].im = b[ji1].im;
314 b[ji1].im = fi;
315 }
316 ij++;
317 }
318
319 if(in==FORWARD) // take conjugates & unscramble outputs
320 {
321 for(i=0,j=n/2;j<n;i++,j++)
322 {
323 r = b[j].re;
324 fi = b[j].im;
325 b[j].re = b[i].re;
326 b[j].im = -b[i].im;
327 b[i].re = r;
328 b[i].im = -fi;
329 }
330 }
331
332 if(in==INVERSE) // scale outputs
333 {
334 for(i=0;i<nthpo;i++)
335 {
336 b[i].re *= fn;
337 b[i].im *= fn;
338 }
339 }
340}
341
342
343/*
344** radix 2 iteration subroutine
345*/
346void R2TX(S32 nthpo, DPCOMPLEX *c0, DPCOMPLEX *c1)
347{
348 S32 k,kk;
349 F64 *cr0, *ci0, *cr1, *ci1, r1, fi1;
350
351 cr0 = &(c0[0].re);
352 ci0 = &(c0[0].im);
353 cr1 = &(c1[0].re);
354 ci1 = &(c1[0].im);
355
356 for(k = 0; k < nthpo; k += 2)
357 {
358 kk = k*2;
359
360 r1 = cr0[kk] + cr1[kk];
361 cr1[kk] = cr0[kk] - cr1[kk];
362 cr0[kk] = r1;
363 fi1 = ci0[kk] + ci1[kk];
364 ci1[kk] = ci0[kk] - ci1[kk];
365 ci0[kk] = fi1;
366 }
367}
368
369
370/*
371** radix 4 iteration subroutine
372*/
373void R4TX(S32 nthpo, DPCOMPLEX *c0, DPCOMPLEX *c1, DPCOMPLEX *c2, DPCOMPLEX *c3)
374{
375 S32 k,kk;
376 F64 *cr0, *ci0, *cr1, *ci1, *cr2, *ci2, *cr3, *ci3;
377 F64 r1,r2,r3,r4,i1,i2,i3,i4;
378
379 cr0 = &(c0[0].re);
380 cr1 = &(c1[0].re);
381 cr2 = &(c2[0].re);
382 cr3 = &(c3[0].re);
383 ci0 = &(c0[0].im);
384 ci1 = &(c1[0].im);
385 ci2 = &(c2[0].im);
386 ci3 = &(c3[0].im);
387
388 for(k = 1; k <= nthpo; k += 4)
389 {
390 kk = (k-1)*2; /* real and imag parts alternate */
391
392 r1 = cr0[kk] + cr2[kk];
393 r2 = cr0[kk] - cr2[kk];
394 r3 = cr1[kk] + cr3[kk];
395 r4 = cr1[kk] - cr3[kk];
396 i1 = ci0[kk] + ci2[kk];
397 i2 = ci0[kk] - ci2[kk];
398 i3 = ci1[kk] + ci3[kk];
399 i4 = ci1[kk] - ci3[kk];
400 cr0[kk] = r1 + r3;
401 ci0[kk] = i1 + i3;
402 cr1[kk] = r1 - r3;
403 ci1[kk] = i1 - i3;
404 cr2[kk] = r2 - i4;
405 ci2[kk] = i2 + r4;
406 cr3[kk] = r2 + i4;
407 ci3[kk] = i2 - r4;
408 }
409}
410
411
412
413/*
414** radix 8 iteration subroutine
415*/
416void R8TX(S32 nxtlt, S32 nthpo, S32 lengt, DPCOMPLEX *cc0, DPCOMPLEX *cc1, DPCOMPLEX *cc2,
417 DPCOMPLEX *cc3, DPCOMPLEX *cc4, DPCOMPLEX *cc5, DPCOMPLEX *cc6, DPCOMPLEX *cc7)
418{
419 S32 j,k,kk;
420 F64 scale, arg, tr, ti;
421 F64 c1,c2,c3,c4,c5,c6,c7;
422 F64 s1,s2,s3,s4,s5,s6,s7;
423 F64 ar0,ar1,ar2,ar3,ar4,ar5,ar6,ar7;
424 F64 ai0,ai1,ai2,ai3,ai4,ai5,ai6,ai7;
425 F64 br0,br1,br2,br3,br4,br5,br6,br7;
426 F64 bi0,bi1,bi2,bi3,bi4,bi5,bi6,bi7;
427
428 F64 *cr0,*cr1,*cr2,*cr3,*cr4,*cr5,*cr6,*cr7;
429 F64 *ci0,*ci1,*ci2,*ci3,*ci4,*ci5,*ci6,*ci7;
430
431 cr0 = &(cc0[0].re);
432 cr1 = &(cc1[0].re);
433 cr2 = &(cc2[0].re);
434 cr3 = &(cc3[0].re);
435 cr4 = &(cc4[0].re);
436 cr5 = &(cc5[0].re);
437 cr6 = &(cc6[0].re);
438 cr7 = &(cc7[0].re);
439
440 ci0 = &(cc0[0].im);
441 ci1 = &(cc1[0].im);
442 ci2 = &(cc2[0].im);
443 ci3 = &(cc3[0].im);
444 ci4 = &(cc4[0].im);
445 ci5 = &(cc5[0].im);
446 ci6 = &(cc6[0].im);
447 ci7 = &(cc7[0].im);
448
449
450 scale = F_TWO_PI/lengt;
451
452 for(j = 1; j <= nxtlt; j++)
453 {
454 arg = (j-1)*scale;
455 c1 = cos(arg);
456 s1 = sin(arg);
457 c2 = c1*c1 - s1*s1;
458 s2 = c1*s1 + c1*s1;
459 c3 = c1*c2 - s1*s2;
460 s3 = c2*s1 + s2*c1;
461 c4 = c2*c2 - s2*s2;
462 s4 = c2*s2 + c2*s2;
463 c5 = c2*c3 - s2*s3;
464 s5 = c3*s2 + s3*c2;
465 c6 = c3*c3 - s3*s3;
466 s6 = c3*s3 + c3*s3;
467 c7 = c3*c4 - s3*s4;
468 s7 = c4*s3 + s4*c3;
469
470 for(k = j; k <= nthpo; k += lengt)
471 {
472 kk = (k-1)*2; /* index by twos; re & im alternate */
473
474 ar0 = cr0[kk] + cr4[kk];
475 ar1 = cr1[kk] + cr5[kk];
476 ar2 = cr2[kk] + cr6[kk];
477 ar3 = cr3[kk] + cr7[kk];
478 ar4 = cr0[kk] - cr4[kk];
479 ar5 = cr1[kk] - cr5[kk];
480 ar6 = cr2[kk] - cr6[kk];
481 ar7 = cr3[kk] - cr7[kk];
482 ai0 = ci0[kk] + ci4[kk];
483 ai1 = ci1[kk] + ci5[kk];
484 ai2 = ci2[kk] + ci6[kk];
485 ai3 = ci3[kk] + ci7[kk];
486 ai4 = ci0[kk] - ci4[kk];
487 ai5 = ci1[kk] - ci5[kk];
488 ai6 = ci2[kk] - ci6[kk];
489 ai7 = ci3[kk] - ci7[kk];
490 br0 = ar0 + ar2;
491 br1 = ar1 + ar3;
492 br2 = ar0 - ar2;
493 br3 = ar1 - ar3;
494 br4 = ar4 - ai6;
495 br5 = ar5 - ai7;
496 br6 = ar4 + ai6;
497 br7 = ar5 + ai7;
498 bi0 = ai0 + ai2;
499 bi1 = ai1 + ai3;
500 bi2 = ai0 - ai2;
501 bi3 = ai1 - ai3;
502 bi4 = ai4 + ar6;
503 bi5 = ai5 + ar7;
504 bi6 = ai4 - ar6;
505 bi7 = ai5 - ar7;
506 cr0[kk] = br0 + br1;
507 ci0[kk] = bi0 + bi1;
508 if(j > 1)
509 {
510 cr1[kk] = c4*(br0-br1) - s4*(bi0-bi1);
511 cr2[kk] = c2*(br2-bi3) - s2*(bi2+br3);
512 cr3[kk] = c6*(br2+bi3) - s6*(bi2-br3);
513 ci1[kk] = c4*(bi0-bi1) + s4*(br0-br1);
514 ci2[kk] = c2*(bi2+br3) + s2*(br2-bi3);
515 ci3[kk] = c6*(bi2-br3) + s6*(br2+bi3);
516 tr = OO_SQRT2*(br5-bi5);
517 ti = OO_SQRT2*(br5+bi5);
518 cr4[kk] = c1*(br4+tr) - s1*(bi4+ti);
519 ci4[kk] = c1*(bi4+ti) + s1*(br4+tr);
520 cr5[kk] = c5*(br4-tr) - s5*(bi4-ti);
521 ci5[kk] = c5*(bi4-ti) + s5*(br4-tr);
522 tr = -OO_SQRT2*(br7+bi7);
523 ti = OO_SQRT2*(br7-bi7);
524 cr6[kk] = c3*(br6+tr) - s3*(bi6+ti);
525 ci6[kk] = c3*(bi6+ti) + s3*(br6+tr);
526 cr7[kk] = c7*(br6-tr) - s7*(bi6-ti);
527 ci7[kk] = c7*(bi6-ti) + s7*(br6-tr);
528 }
529 else
530 {
531 cr1[kk] = br0 - br1;
532 cr2[kk] = br2 - bi3;
533 cr3[kk] = br2 + bi3;
534 ci1[kk] = bi0 - bi1;
535 ci2[kk] = bi2 + br3;
536 ci3[kk] = bi2 - br3;
537 tr = OO_SQRT2*(br5-bi5);
538 ti = OO_SQRT2*(br5+bi5);
539 cr4[kk] = br4 + tr;
540 ci4[kk] = bi4 + ti;
541 cr5[kk] = br4 - tr;
542 ci5[kk] = bi4 - ti;
543 tr = -OO_SQRT2*(br7+bi7);
544 ti = OO_SQRT2*(br7-bi7);
545 cr6[kk] = br6 + tr;
546 ci6[kk] = bi6 + ti;
547 cr7[kk] = br6 - tr;
548 ci7[kk] = bi6 - ti;
549 }
550 }
551 }
552}
553
554
555/* see if exactly one bit is set in integer argument */
556S32 power_of_2(S32 n)
557{
558 S32 bits=0;
559 while(n)
560 {
561 bits += n & 1;
562 n >>= 1;
563 }
564 return(bits==1);
565}
566
567/* get binary log of integer argument; exact if n a power of 2 */
568S32 fastlog2(S32 n)
569{
570 S32 log = -1;
571 while(n)
572 {
573 log++;
574 n >>= 1;
575 }
576 return(log);
577}
generated by cgit v1.1 at 2024-11-15 14:47:36 +0000