ese 519
Diff: fhtdit2.cc
- Revision:
- 0:2076b4d80327
diff -r 000000000000 -r 2076b4d80327 fhtdit2.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fhtdit2.cc Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,130 @@ +// This file is part of the FXT library. +// Copyright (C) 2010, 2012 Joerg Arndt +// License: GNU General Public License version 3 or later, +// see the file COPYING.txt in the main directory. + +#include "fxttypes.h" +#include "complextype.h" +#include "sincos.h" +#include "revbinpermute.h" +#include "sumdiff.h" + +#include <cmath> // M_PI + + +void +fht_depth_first_dit2(double *f, ulong ldn) +// Radix-2 decimation in time (DIT) FHT. +// Depth-first version. +// Compared to usual fht this one +// - does more trig computations +// - is (far) better memory local +{ + const ulong n = 1UL<<ldn; + revbin_permute(f, n); + + for (ulong ldm=1; ldm<=ldn; ++ldm) + { + const ulong m = (1UL<<ldm); + const ulong mh = (m>>1); + const ulong m4 = (mh>>1); + const double phi0 = M_PI/(double)mh; + + for (ulong r=0; r<n; r+=m) + { + { // j == 0: + ulong t1 = r; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + + if ( m4 ) + { + ulong t1 = r + m4; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + + for (ulong j=1, k=mh-1; j<k; ++j, --k) + { + double s, c; + SinCos(phi0*(double)j, &s, &c); + + ulong tj = r + mh + j; + ulong tk = r + mh + k; + double fj = f[tj]; + double fk = f[tk]; + f[tj] = fj * c + fk * s; + f[tk] = fj * s - fk * c; + + ulong t1 = r + j; + ulong t2 = tj; // == t1 + mh; + sumdiff(f[t1], f[t2]); + + t1 = r + k; + t2 = tk; // == t1 + mh; + sumdiff(f[t1], f[t2]); + } + } + } +} +// ------------------------- + + +void +fht_dit2(double *f, ulong ldn) +// Radix-2 decimation in time (DIT) FHT. +{ + const ulong n = 1UL<<ldn; + revbin_permute(f, n); + + for (ulong ldm=1; ldm<=ldn; ++ldm) + { + const ulong m = (1UL<<ldm); + const ulong mh = (m>>1); + const ulong m4 = (mh>>1); + const double phi0 = M_PI/(double)mh; + + for (ulong r=0; r<n; r+=m) + { + { // j == 0: + ulong t1 = r; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + + if ( m4 ) + { + ulong t1 = r + m4; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + } + + for (ulong j=1, k=mh-1; j<k; ++j, --k) + { + double s, c; + SinCos(phi0*(double)j, &s, &c); + + for (ulong r=0; r<n; r+=m) + { + ulong tj = r + mh + j; + ulong tk = r + mh + k; + double fj = f[tj]; + double fk = f[tk]; + f[tj] = fj * c + fk * s; + f[tk] = fj * s - fk * c; + + ulong t1 = r + j; + ulong t2 = tj; // == t1 + mh; + sumdiff(f[t1], f[t2]); + + t1 = r + k; + t2 = tk; // == t1 + mh; + sumdiff(f[t1], f[t2]); + } + } + } +} +// ------------------------- +