Diff: fhtdit2.cc

Revision:

0:2076b4d80327

--- /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]);
+            }
+        }
+    }
+}
+// -------------------------
+