FFTS is no longer depended on any other math library, and this should help to verify its numerical accuracy.

1 files changed, 3 insertions, 77 deletions

diff --git a/src/ffts_real.c b/src/ffts_real.c
index f1355c7..f6e6127 100644
--- a/src/ffts_real.c
+++ b/src/ffts_real.c
@@ -34,6 +34,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "ffts_real.h"
 #include "ffts_internal.h"
+#include "ffts_trig.h"
 
 #ifdef HAVE_NEON
 #include <arm_neon.h>
@@ -604,7 +605,6 @@ ffts_plan_t*
 ffts_init_1d_real(size_t N, int sign)
 {
     ffts_plan_t *p;
-    size_t i;
 
     p = (ffts_plan_t*) calloc(1, sizeof(*p) + sizeof(*p->plans));
     if (!p) {
@@ -642,85 +642,11 @@ ffts_init_1d_real(size_t N, int sign)
         goto cleanup;
     }
 
-    if (sign < 0) {
-        /* peel off the first */
-        p->A[0] =  0.5f;
-        p->A[1] = -0.5f;
-#ifdef HAVE_SSE3
-        p->B[0] = -0.5f;
-#else
-        p->B[0] =  0.5f;
-#endif
-        p->B[1] =  0.5f;
-
-        for (i = 1; i < N/4; i++) {
-            float t0 = (float) (0.5 * (1.0 - sin(2.0 * M_PI / N * i)));
-            float t1 = (float) (0.5 * (1.0 * cos(2.0 * M_PI / N * i)));
-            float t2 = (float) (0.5 * (1.0 + sin(2.0 * M_PI / N * i)));
-
-            p->A[    2 * i + 0] =  t0;
-            p->A[N - 2 * i + 0] =  t0;
-            p->A[    2 * i + 1] = -t1;
-            p->A[N - 2 * i + 1] =  t1;
-
 #ifdef HAVE_SSE3
-            p->B[    2 * i + 0] = -t2;
-            p->B[N - 2 * i + 0] = -t2;
+	ffts_generate_table_1d_real_32f(p, sign, 1);
 #else
-            p->B[    2 * i + 0] =  t2;
-            p->B[N - 2 * i + 0] =  t2;
+	ffts_generate_table_1d_real_32f(p, sign, 0);
 #endif
-            p->B[    2 * i + 1] =  t1;
-            p->B[N - 2 * i + 1] = -t1;
-        }
-
-        /* and the last */
-        p->A[2 * i + 0] =  0.0f;
-        p->A[2 * i + 1] =  0.0f;
-#ifdef HAVE_SSE3
-        p->B[2 * i + 0] = -1.0f;
-#else
-        p->B[2 * i + 0] =  1.0f;
-#endif
-        p->B[2 * i + 1] =  0.0f;
-    } else {
-        /* peel of the first */
-        p->A[0] =  1.0f;
-#ifdef HAVE_SSE3
-        p->A[1] =  1.0f;
-#else
-        p->A[1] = -1.0f;
-#endif
-        p->B[0] =  1.0f;
-        p->B[1] =  1.0f;
-
-        for (i = 1; i < N/4; i++) {
-            float t0 = (float) (1.0 * (1.0 - sin(2.0 * M_PI / N * i)));
-            float t1 = (float) (1.0 * (1.0 * cos(2.0 * M_PI / N * i)));
-            float t2 = (float) (1.0 * (1.0 + sin(2.0 * M_PI / N * i)));
-
-            p->A[    2 * i + 0] =  t0;
-            p->A[N - 2 * i + 0] =  t0;
-#ifdef HAVE_SSE3
-            p->A[    2 * i + 1] =  t1;
-            p->A[N - 2 * i + 1] = -t1;
-#else
-            p->A[    2 * i + 1] = -t1;
-            p->A[N - 2 * i + 1] =  t1;
-#endif
-
-            p->B[    2 * i + 0] =  t2;
-            p->B[N - 2 * i + 0] =  t2;
-            p->B[    2 * i + 1] =  t1;
-            p->B[N - 2 * i + 1] = -t1;
-        }
-
-        /* and the last */
-        p->A[2 * i + 0] = 0.0f;
-        p->A[2 * i + 1] = 0.0f;
-        p->B[2 * i + 0] = 2.0f;
-        p->B[2 * i + 1] = 0.0f;
-    }
 
     return p;