70 lines
2.1 KiB
C
70 lines
2.1 KiB
C
|
|
#include "pffft_double.h"
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
|
|
|
|
void c_forward_complex_double(const int transformLen)
|
|
{
|
|
printf("running %s()\n", __FUNCTION__);
|
|
|
|
/* first check - might be skipped */
|
|
if (transformLen < pffftd_min_fft_size(PFFFT_COMPLEX))
|
|
{
|
|
fprintf(stderr, "Error: minimum FFT transformation length is %d\n", pffftd_min_fft_size(PFFFT_COMPLEX));
|
|
return;
|
|
}
|
|
|
|
/* instantiate FFT and prepare transformation for length N */
|
|
PFFFTD_Setup *ffts = pffftd_new_setup(transformLen, PFFFT_COMPLEX);
|
|
|
|
/* one more check */
|
|
if (!ffts)
|
|
{
|
|
fprintf(stderr,
|
|
"Error: transformation length %d is not decomposable into small prime factors. "
|
|
"Next valid transform size is: %d ; next power of 2 is: %d\n",
|
|
transformLen,
|
|
pffftd_nearest_transform_size(transformLen, PFFFT_COMPLEX, 1),
|
|
pffftd_next_power_of_two(transformLen) );
|
|
return;
|
|
}
|
|
|
|
/* allocate aligned vectors for input X and output Y */
|
|
double *X = (double*)pffftd_aligned_malloc(transformLen * 2 * sizeof(double)); /* complex: re/im interleaved */
|
|
double *Y = (double*)pffftd_aligned_malloc(transformLen * 2 * sizeof(double)); /* complex: re/im interleaved */
|
|
double *W = (double*)pffftd_aligned_malloc(transformLen * 2 * sizeof(double));
|
|
|
|
/* prepare some input data */
|
|
for (int k = 0; k < 2 * transformLen; k += 4)
|
|
{
|
|
X[k] = k / 2; /* real */
|
|
X[k+1] = (k / 2) & 1; /* imag */
|
|
|
|
X[k+2] = -1 - k / 2; /* real */
|
|
X[k+3] = (k / 2) & 1; /* imag */
|
|
}
|
|
|
|
/* do the forward transform; write complex spectrum result into Y */
|
|
pffftd_transform_ordered(ffts, X, Y, W, PFFFT_FORWARD);
|
|
|
|
/* print spectral output */
|
|
printf("output should be complex spectrum with %d complex bins\n", transformLen);
|
|
for (int k = 0; k < 2 * transformLen; k += 2)
|
|
printf("Y[%d] = %f + i * %f\n", k/2, Y[k], Y[k+1]);
|
|
|
|
pffftd_aligned_free(W);
|
|
pffftd_aligned_free(Y);
|
|
pffftd_aligned_free(X);
|
|
pffftd_destroy_setup(ffts);
|
|
}
|
|
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
int N = (1 < argc) ? atoi(argv[1]) : 16;
|
|
c_forward_complex_double(N);
|
|
return 0;
|
|
}
|