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