1
2 #include "pffft_double.h"
3
4 #include <stdio.h>
5 #include <stdlib.h>
6
7
c_forward_complex_double(const int transformLen)8 void c_forward_complex_double(const int transformLen)
9 {
10 printf("running %s()\n", __FUNCTION__);
11
12 /* first check - might be skipped */
13 if (transformLen < pffftd_min_fft_size(PFFFT_COMPLEX))
14 {
15 fprintf(stderr, "Error: minimum FFT transformation length is %d\n", pffftd_min_fft_size(PFFFT_COMPLEX));
16 return;
17 }
18
19 /* instantiate FFT and prepare transformation for length N */
20 PFFFTD_Setup *ffts = pffftd_new_setup(transformLen, PFFFT_COMPLEX);
21
22 /* one more check */
23 if (!ffts)
24 {
25 fprintf(stderr,
26 "Error: transformation length %d is not decomposable into small prime factors. "
27 "Next valid transform size is: %d ; next power of 2 is: %d\n",
28 transformLen,
29 pffftd_nearest_transform_size(transformLen, PFFFT_COMPLEX, 1),
30 pffftd_next_power_of_two(transformLen) );
31 return;
32 }
33
34 /* allocate aligned vectors for input X and output Y */
35 double *X = (double*)pffftd_aligned_malloc(transformLen * 2 * sizeof(double)); /* complex: re/im interleaved */
36 double *Y = (double*)pffftd_aligned_malloc(transformLen * 2 * sizeof(double)); /* complex: re/im interleaved */
37 double *W = (double*)pffftd_aligned_malloc(transformLen * 2 * sizeof(double));
38
39 /* prepare some input data */
40 for (int k = 0; k < 2 * transformLen; k += 4)
41 {
42 X[k] = k / 2; /* real */
43 X[k+1] = (k / 2) & 1; /* imag */
44
45 X[k+2] = -1 - k / 2; /* real */
46 X[k+3] = (k / 2) & 1; /* imag */
47 }
48
49 /* do the forward transform; write complex spectrum result into Y */
50 pffftd_transform_ordered(ffts, X, Y, W, PFFFT_FORWARD);
51
52 /* print spectral output */
53 printf("output should be complex spectrum with %d complex bins\n", transformLen);
54 for (int k = 0; k < 2 * transformLen; k += 2)
55 printf("Y[%d] = %f + i * %f\n", k/2, Y[k], Y[k+1]);
56
57 pffftd_aligned_free(W);
58 pffftd_aligned_free(Y);
59 pffftd_aligned_free(X);
60 pffftd_destroy_setup(ffts);
61 }
62
63
main(int argc,char * argv[])64 int main(int argc, char *argv[])
65 {
66 int N = (1 < argc) ? atoi(argv[1]) : 16;
67 c_forward_complex_double(N);
68 return 0;
69 }
70