}
void
-gcp_fft_1d_r2_forward_twiddle(unsigned int log2_fft_length, std::complex<float> *W)
+gcp_fft_1d_r2_twiddle(unsigned int log2_fft_length, std::complex<float> *W)
{
unsigned int n = 1 << log2_fft_length;
W[n/4 - i].imag() = -W[i].real();
}
}
-
-
-void
-gcp_fft_1d_r2_reverse_twiddle(unsigned int log2_fft_length, std::complex<float> *W)
-{
- // FIXME this is wrong/insufficient. inverse is still incorrect
-
- // reverse factors are the conjugate of the forward factors
- gcp_fft_1d_r2_forward_twiddle(log2_fft_length, W);
-
- unsigned int n = 1 << log2_fft_length;
- for (unsigned i=0; i < n/4; i++)
- W[i] = conj(W[i]);
-}
* \brief Submit a job that computes the forward or reverse FFT.
*
* \param mgr is the job manager instance
- * \param log2_fft_length is the log2 of the fft_length (4 <= x <= 13).
+ * \param log2_fft_length is the log2 of the fft_length (4 <= x <= 12).
* \param forward is true to compute the forward xform
* \param out is the fft_length output from FFT (must be 16-byte aligned).
* \param in is the fft_length input to FFT (must be 16-byte aligned).
const std::complex<float> *W);
/*!
- * \brief Compute twiddle factors for forward transform.
+ * \brief Compute twiddle factors
*
* \param log2_fft_length is the log2 of the fft_length.
* \param W is fft_length/4 twiddle factor output (must be 16-byte aligned).
*/
void
-gcp_fft_1d_r2_forward_twiddle(unsigned int log2_fft_length, std::complex<float> *W);
-
-/*!
- * \brief Compute twiddle factors for reverse transform.
- *
- * \param log2_fft_length is the log2 of the fft_length.
- * \param W is fft_length/4 twiddle factor output (must be 16-byte aligned).
- */
-void
-gcp_fft_1d_r2_reverse_twiddle(unsigned int log2_fft_length, std::complex<float> *W);
+gcp_fft_1d_r2_twiddle(unsigned int log2_fft_length, std::complex<float> *W);
#endif /* INCLUDED_GCP_FFT_1D_R2_H */
std::complex<float> *cell_out = (std::complex<float> *) cell_out_void.get();
std::complex<float> *cell_twiddle = (std::complex<float> *) cell_twiddle_void.get();
- if (forward)
- gcp_fft_1d_r2_forward_twiddle(log2_fft_size, cell_twiddle);
- else
- gcp_fft_1d_r2_reverse_twiddle(log2_fft_size, cell_twiddle);
+ gcp_fft_1d_r2_twiddle(log2_fft_size, cell_twiddle);
srandom(1); // we want reproducibility
fprintf(stdout, "%s fft_size = %4d max_rel_error = %e\n",
forward ? "fwd" : "rev", fft_size, max_rel);
- // CPPUNIT_ASSERT(max_rel <= 1e-4);
-
+ CPPUNIT_ASSERT(max_rel <= 5e-3);
}
#include <gc_declare_proc.h>
#include <libfft.h>
+/*
+ * v is really vector complex<float>
+ */
+static void
+conjugate_vector(vector float *v, int nelements)
+{
+ vector float k = {1, -1, 1, -1};
+ int i;
+ for (i = 0; i < nelements; i++)
+ v[i] *= k;
+}
+
static void
gcs_fft_1d_r2(const gc_job_direct_args_t *input,
gc_job_direct_args_t *output __attribute__((unused)),
vector float *in = (vector float *) eaa->arg[1].ls_addr;
vector float *W = (vector float *) eaa->arg[2].ls_addr;
int log2_fft_length = input->arg[0].u32;
- int forward = input->arg[1].u32; // non-zero if forward xform (FIXME use)
+ int forward = input->arg[1].u32; // non-zero if forward xform
- fft_1d_r2(out, in, W, log2_fft_length);
+ if (forward){
+ fft_1d_r2(out, in, W, log2_fft_length);
+ }
+ else {
+ conjugate_vector(in, 1 << (log2_fft_length - 1));
+ fft_1d_r2(out, in, W, log2_fft_length);
+ conjugate_vector(out, 1 << (log2_fft_length - 1));
+ }
}
GC_DECLARE_PROC(gcs_fft_1d_r2, "fft_1d_r2");