summaryrefslogtreecommitdiff
path: root/gsl-1.9/fft/test_complex_source.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/fft/test_complex_source.c')
-rw-r--r--gsl-1.9/fft/test_complex_source.c448
1 files changed, 448 insertions, 0 deletions
diff --git a/gsl-1.9/fft/test_complex_source.c b/gsl-1.9/fft/test_complex_source.c
new file mode 100644
index 0000000..a4c0d27
--- /dev/null
+++ b/gsl-1.9/fft/test_complex_source.c
@@ -0,0 +1,448 @@
+/* fft/test_complex.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "bitreverse.h"
+#include "signals.h"
+#include "compare.h"
+
+void FUNCTION(test_complex,func) (size_t stride, size_t n);
+int FUNCTION(test,offset) (const BASE data[], size_t stride,
+ size_t n, size_t offset);
+void FUNCTION(test_complex,bitreverse_order) (size_t stride, size_t n) ;
+void FUNCTION(test_complex,radix2) (size_t stride, size_t n);
+
+int FUNCTION(test,offset) (const BASE data[], size_t stride,
+ size_t n, size_t offset)
+{
+ int status = 0 ;
+ size_t i, j, k = 0 ;
+
+ for (i = 0; i < n; i++)
+ {
+ k += 2 ;
+
+ for (j = 1; j < stride; j++)
+ {
+ status |= data[k] != k + offset ;
+ k++ ;
+ status |= data[k] != k + offset ;
+ k++ ;
+ }
+ }
+ return status ;
+}
+
+
+void FUNCTION(test_complex,func) (size_t stride, size_t n)
+{
+ size_t i ;
+ int status ;
+
+ TYPE(gsl_fft_complex_wavetable) * cw ;
+ TYPE(gsl_fft_complex_workspace) * cwork ;
+
+ BASE * complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * fft_complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * fft_complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+
+ for (i = 0 ; i < 2 * n * stride ; i++)
+ {
+ complex_data[i] = (BASE)i ;
+ complex_tmp[i] = (BASE)(i + 1000.0) ;
+ fft_complex_data[i] = (BASE)(i + 2000.0) ;
+ fft_complex_tmp[i] = (BASE)(i + 3000.0) ;
+ }
+
+ gsl_set_error_handler (NULL); /* abort on any errors */
+
+ /* Test allocation */
+
+ {
+ cw = FUNCTION(gsl_fft_complex_wavetable,alloc) (n);
+ gsl_test (cw == 0, NAME(gsl_fft_complex_wavetable)
+ "_alloc, n = %d, stride = %d", n, stride);
+ }
+
+ {
+ cwork = FUNCTION(gsl_fft_complex_workspace,alloc) (n);
+ gsl_test (cwork == 0, NAME(gsl_fft_complex_workspace)
+ "_alloc, n = %d", n);
+ }
+
+
+ /* Test mixed radix fft with noise */
+
+ {
+ FUNCTION(fft_signal,complex_noise) (n, stride, complex_data, fft_complex_data);
+ for (i = 0 ; i < n ; i++)
+ {
+ REAL(complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
+ IMAG(complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
+ }
+
+ FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
+
+ for (i = 0 ; i < n ; i++)
+ {
+ REAL(fft_complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
+ IMAG(fft_complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
+ }
+
+ status = FUNCTION(compare_complex,results) ("dft", fft_complex_data,
+ "fft of noise", complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_forward with signal_noise, n = %d, stride = %d", n, stride);
+
+ if (stride > 1)
+ {
+ status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_forward avoids unstrided data, n = %d, stride = %d",
+ n, stride);
+ }
+ }
+
+ /* Test the inverse fft */
+
+ {
+ status = FUNCTION(gsl_fft_complex,inverse) (complex_data, stride, n, cw, cwork);
+ status = FUNCTION(compare_complex,results) ("orig", complex_tmp,
+ "fft inverse", complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_inverse with signal_noise, n = %d, stride = %d", n, stride);
+
+ if (stride > 1)
+ {
+ status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_inverse other data untouched, n = %d, stride = %d",
+ n, stride);
+ }
+
+ }
+
+ /* Test the backward fft */
+
+ {
+ status = FUNCTION(gsl_fft_complex,backward) (fft_complex_tmp, stride, n, cw, cwork);
+
+ for (i = 0; i < n; i++)
+ {
+ REAL(complex_tmp,stride,i) *= n;
+ IMAG(complex_tmp,stride,i) *= n;
+ }
+ status = FUNCTION(compare_complex,results) ("orig",
+ complex_tmp,
+ "fft backward",
+ fft_complex_tmp,
+ stride, n, 1e6);
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_backward with signal_noise, n = %d, stride = %d", n, stride);
+
+ if (stride > 1)
+ {
+ status = FUNCTION(test, offset) (fft_complex_tmp, stride, n, 3000) ;
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_backward avoids unstrided data, n = %d, stride = %d",
+ n, stride);
+ }
+
+ }
+
+ /* Test a pulse signal */
+
+ {
+ FUNCTION(fft_signal,complex_pulse) (1, n, stride, 1.0, 0.0, complex_data,
+ fft_complex_data);
+ FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
+ status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
+ "fft of pulse", complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_forward with signal_pulse, n = %d, stride = %d", n, stride);
+
+ }
+
+
+ /* Test a constant signal */
+
+ {
+ FUNCTION(fft_signal,complex_constant) (n, stride, 1.0, 0.0, complex_data,
+ fft_complex_data);
+ FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
+ status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
+ "fft of constant",
+ complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_forward with signal_constant, n = %d, stride = %d", n, stride);
+ }
+
+ /* Test an exponential (cos/sin) signal */
+
+ {
+ status = 0;
+ for (i = 0; i < n; i++)
+ {
+ FUNCTION(fft_signal,complex_exp) ((int)i, n, stride, 1.0, 0.0, complex_data,
+ fft_complex_data);
+ FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
+ status |= FUNCTION(compare_complex,results) ("analytic",
+ fft_complex_data,
+ "fft of exp",
+ complex_data,
+ stride, n, 1e6);
+ }
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_forward with signal_exp, n = %d, stride = %d", n, stride);
+ }
+
+ FUNCTION(gsl_fft_complex_wavetable,free) (cw);
+ FUNCTION(gsl_fft_complex_workspace,free) (cwork);
+
+ free (complex_data);
+ free (complex_tmp);
+ free (fft_complex_data);
+ free (fft_complex_tmp);
+}
+
+
+void
+FUNCTION(test_complex,bitreverse_order) (size_t stride, size_t n)
+{
+ int status ;
+ size_t logn, i ;
+
+ BASE * tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * reversed_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+
+ for (i = 0; i < 2 * stride * n; i++)
+ {
+ data[i] = (BASE)i ;
+ }
+
+ memcpy (tmp, data, 2 * n * stride * sizeof(BASE)) ;
+
+ logn = 0 ; while (n > (1U<<logn)) {logn++ ; } ;
+
+ /* do a naive bit reversal as a baseline for testing the other routines */
+
+ for (i = 0; i < n; i++)
+ {
+ size_t i_tmp = i ;
+ size_t j = 0 ;
+ size_t bit ;
+
+ for (bit = 0; bit < logn; bit++)
+ {
+ j <<= 1; /* reverse shift i into j */
+ j |= i_tmp & 1;
+ i_tmp >>= 1;
+ }
+
+ reversed_data[2*j*stride] = data[2*i*stride] ;
+ reversed_data[2*j*stride+1] = data[2*i*stride+1] ;
+ }
+
+ FUNCTION(fft_complex,bitreverse_order) (data, stride, n, logn);
+
+ status = FUNCTION(compare_complex,results) ("naive bit reverse",
+ reversed_data,
+ "fft_complex_bitreverse_order",
+ data,
+ stride, n, 1e6);
+
+ gsl_test (status, "fft_complex_bitreverse_order, n = %d", n);
+
+ free (reversed_data) ;
+ free (data) ;
+ free (tmp) ;
+}
+
+void FUNCTION(test_complex,radix2) (size_t stride, size_t n)
+{
+ size_t i ;
+ int status ;
+
+ BASE * complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * fft_complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+ BASE * fft_complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
+
+ for (i = 0 ; i < 2 * n * stride ; i++)
+ {
+ complex_data[i] = (BASE)i ;
+ complex_tmp[i] = (BASE)(i + 1000.0) ;
+ fft_complex_data[i] = (BASE)(i + 2000.0) ;
+ fft_complex_tmp[i] = (BASE)(i + 3000.0) ;
+ }
+
+ gsl_set_error_handler (NULL); /* abort on any errors */
+
+ /* Test radix-2 fft with noise */
+
+ {
+ FUNCTION(fft_signal,complex_noise) (n, stride, complex_data,
+ fft_complex_data);
+ for (i = 0 ; i < n ; i++)
+ {
+ REAL(complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
+ IMAG(complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
+ }
+
+ FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
+
+ for (i = 0 ; i < n ; i++)
+ {
+ REAL(fft_complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
+ IMAG(fft_complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
+ }
+
+ status = FUNCTION(compare_complex,results) ("dft", fft_complex_data,
+ "fft of noise", complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_forward with signal_noise, n = %d, stride = %d",
+ n, stride);
+
+ if (stride > 1)
+ {
+ status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_forward avoids unstrided data, n = %d, stride = %d",
+ n, stride);
+ }
+ }
+
+ /* Test the inverse fft */
+
+ {
+ status = FUNCTION(gsl_fft_complex,radix2_inverse) (complex_data, stride, n);
+ status = FUNCTION(compare_complex,results) ("orig", complex_tmp,
+ "fft inverse", complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_inverse with signal_noise, n = %d, stride = %d", n, stride);
+
+ if (stride > 1)
+ {
+ status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_inverse other data untouched, n = %d, stride = %d",
+ n, stride);
+ }
+
+ }
+
+ /* Test the backward fft */
+
+ {
+ status = FUNCTION(gsl_fft_complex,radix2_backward) (fft_complex_tmp, stride, n);
+
+ for (i = 0; i < n; i++)
+ {
+ REAL(complex_tmp,stride,i) *= n;
+ IMAG(complex_tmp,stride,i) *= n;
+ }
+ status = FUNCTION(compare_complex,results) ("orig",
+ complex_tmp,
+ "fft backward",
+ fft_complex_tmp,
+ stride, n, 1e6);
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_backward with signal_noise, n = %d, stride = %d", n, stride);
+
+ if (stride > 1)
+ {
+ status = FUNCTION(test, offset) (fft_complex_tmp, stride, n, 3000) ;
+
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_backward avoids unstrided data, n = %d, stride = %d",
+ n, stride);
+ }
+
+ }
+
+ /* Test a pulse signal */
+
+ {
+ FUNCTION(fft_signal,complex_pulse) (1, n, stride, 1.0, 0.0, complex_data,
+ fft_complex_data);
+ FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
+ status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
+ "fft of pulse", complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_forward with signal_pulse, n = %d, stride = %d", n, stride);
+
+ }
+
+
+ /* Test a constant signal */
+
+ {
+ FUNCTION(fft_signal,complex_constant) (n, stride, 1.0, 0.0, complex_data,
+ fft_complex_data);
+ FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
+ status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
+ "fft of constant",
+ complex_data,
+ stride, n, 1e6);
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_forward with signal_constant, n = %d, stride = %d",
+ n, stride);
+ }
+
+ /* Test an exponential (cos/sin) signal */
+
+ {
+ status = 0;
+ for (i = 0; i < n; i++)
+ {
+ FUNCTION(fft_signal,complex_exp) ((int)i, n, stride, 1.0, 0.0, complex_data,
+ fft_complex_data);
+ FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
+ status |= FUNCTION(compare_complex,results) ("analytic",
+ fft_complex_data,
+ "fft of exp",
+ complex_data,
+ stride, n, 1e6);
+ }
+ gsl_test (status, NAME(gsl_fft_complex)
+ "_radix2_forward with signal_exp, n = %d, stride = %d", n, stride);
+ }
+
+
+ free (complex_data);
+ free (complex_tmp);
+ free (fft_complex_data);
+ free (fft_complex_tmp);
+}
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 07:25:41 +0000