summaryrefslogtreecommitdiffstats
path: root/src/ffts_trig.c
diff options
context:
space:
mode:
authorJukka Ojanen <jukka.ojanen@linkotec.net>2015-07-14 23:51:35 +0300
committerJukka Ojanen <jukka.ojanen@linkotec.net>2015-07-14 23:51:35 +0300
commitf571c435ceccc56e79c024f626a91c57f52d94ff (patch)
tree576044bd788a5ba1d6314aa27ed305e6473e3d92 /src/ffts_trig.c
parent9885d87c6335d5b688bdf6b90e78de1add605d63 (diff)
downloadffts-f571c435ceccc56e79c024f626a91c57f52d94ff.zip
ffts-f571c435ceccc56e79c024f626a91c57f52d94ff.tar.gz
FFTS is no longer depended on any other math library, and this should help to verify its numerical accuracy.
Diffstat (limited to 'src/ffts_trig.c')
-rw-r--r--src/ffts_trig.c130
1 files changed, 127 insertions, 3 deletions
diff --git a/src/ffts_trig.c b/src/ffts_trig.c
index 8af96b9..514a1e5 100644
--- a/src/ffts_trig.c
+++ b/src/ffts_trig.c
@@ -118,7 +118,7 @@ ffts_generate_cosine_sine_32f(ffts_cpx_32f *const table, int table_size)
c[0] = 1.0;
s[0] = 0.0;
- /* generate sine and cosine table with maximum error less than 1 ULP */
+ /* generate sine and cosine tables with maximum error less than 1 ULP */
for (i = 1; i < (table_size + 1)/2; i++) {
c[1] = c[0] - ((alpha * c[0]) + (beta * s[0]));
s[1] = s[0] - ((alpha * s[0]) - (beta * c[0]));
@@ -190,13 +190,13 @@ ffts_generate_cosine_sine_pow2_32f(ffts_cpx_32f *const table, int table_size)
FFTS_ASSUME_ALIGNED_32(&cos_sin_pi_table[4 * offset]);
hs = (const double*) &half_secant[2 * offset];
- /* initialize from table */
+ /* initialize from lookup table */
for (i = 0; i <= log_2; i++) {
w[i][0] = ct[i][0];
w[i][1] = ct[i][1];
}
- /* generate sine and cosine table with maximum error less than 0.5 ULP */
+ /* generate sine and cosine tables with maximum error less than 0.5 ULP */
for (i = 1; i < table_size/2; i++) {
/* calculate trailing zeros in index */
log_2 = ffts_ctzl(i);
@@ -218,4 +218,128 @@ mid_point:
exit:
return 0;
+}
+
+int
+ffts_generate_table_1d_real_32f(struct _ffts_plan_t *const p,
+ int sign,
+ int invert)
+{
+ const ffts_cpx_64f *FFTS_RESTRICT ct;
+ const double *FFTS_RESTRICT hs;
+ ffts_cpx_64f FFTS_ALIGN(16) w[32];
+ int i, log_2, offset, N;
+ float *A, *B;
+
+ if (!p) {
+ return -1;
+ }
+
+ A = (float*) FFTS_ASSUME_ALIGNED_32(p->A);
+ B = (float*) FFTS_ASSUME_ALIGNED_32(p->B);
+ N = (int) p->N;
+
+ /* the first */
+ if (sign < 0) {
+ A[0] = 0.5f;
+ A[1] = -0.5f;
+ B[0] = invert ? -0.5f : 0.5f;
+ B[1] = 0.5f;
+ } else {
+ /* peel of the first */
+ A[0] = 1.0f;
+ A[1] = invert ? 1.0f : -1.0f;
+ B[0] = 1.0f;
+ B[1] = 1.0f;
+ }
+
+ if (FFTS_UNLIKELY(N == 4)) {
+ i = 1;
+ goto last;
+ }
+
+ /* calculate table offset */
+ FFTS_ASSUME(N / 4 > 1);
+ log_2 = ffts_ctzl(N);
+ FFTS_ASSUME(log_2 > 2);
+ offset = 34 - log_2;
+ ct = (const ffts_cpx_64f*)
+ FFTS_ASSUME_ALIGNED_32(&cos_sin_pi_table[4 * offset]);
+ hs = (const double*) &half_secant[2 * offset];
+
+ /* initialize from lookup table */
+ for (i = 0; i <= log_2; i++) {
+ w[i][0] = ct[i][0];
+ w[i][1] = ct[i][1];
+ }
+
+ /* generate sine and cosine tables with maximum error less than 0.5 ULP */
+ if (sign < 0) {
+ for (i = 1; i < N/4; i++) {
+ float t0, t1, t2;
+
+ /* calculate trailing zeros in index */
+ log_2 = ffts_ctzl(i);
+
+ t0 = (float) (0.5 * (1.0 - w[log_2][1]));
+ t1 = (float) (0.5 * w[log_2][0]);
+ t2 = (float) (0.5 * (1.0 + w[log_2][1]));
+
+ A[ 2 * i + 0] = t0;
+ A[N - 2 * i + 0] = t0;
+ A[ 2 * i + 1] = -t1;
+ A[N - 2 * i + 1] = t1;
+
+ B[ 2 * i + 0] = invert ? -t2 : t2;
+ B[N - 2 * i + 0] = invert ? -t2 : t2;
+ B[ 2 * i + 1] = t1;
+ B[N - 2 * i + 1] = -t1;
+
+ /* skip and find next trailing zero */
+ offset = (log_2 + 2 + ffts_ctzl(~i >> (log_2 + 2)));
+ w[log_2][0] = hs[log_2] * (w[log_2 + 1][0] + w[offset][0]);
+ w[log_2][1] = hs[log_2] * (w[log_2 + 1][1] + w[offset][1]);
+ }
+ } else {
+ for (i = 1; i < N/4; i++) {
+ float t0, t1, t2;
+
+ /* calculate trailing zeros in index */
+ log_2 = ffts_ctzl(i);
+
+ t0 = (float) (1.0 - w[log_2][1]);
+ t1 = (float) w[log_2][0];
+ t2 = (float) (1.0 + w[log_2][1]);
+
+ A[ 2 * i + 0] = t0;
+ A[N - 2 * i + 0] = t0;
+ A[ 2 * i + 1] = invert ? t1 : -t1;
+ A[N - 2 * i + 1] = invert ? -t1 : t1;
+
+ B[ 2 * i + 0] = t2;
+ B[N - 2 * i + 0] = t2;
+ B[ 2 * i + 1] = t1;
+ B[N - 2 * i + 1] = -t1;
+
+ /* skip and find next trailing zero */
+ offset = (log_2 + 2 + ffts_ctzl(~i >> (log_2 + 2)));
+ w[log_2][0] = hs[log_2] * (w[log_2 + 1][0] + w[offset][0]);
+ w[log_2][1] = hs[log_2] * (w[log_2 + 1][1] + w[offset][1]);
+ }
+ }
+
+last:
+ if (sign < 0) {
+ A[2 * i + 0] = 0.0f;
+ A[2 * i + 1] = 0.0f;
+ B[2 * i + 0] = invert ? -1.0f : 1.0f;
+ B[2 * i + 1] = 0.0f;
+ } else {
+ A[2 * i + 0] = 0.0f;
+ A[2 * i + 1] = 0.0f;
+ B[2 * i + 0] = 2.0f;
+ B[2 * i + 1] = 0.0f;
+ }
+
+ return 0;
} \ No newline at end of file
OpenPOWER on IntegriCloud