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