1 /*
  2   Radix-2 DIF recursion, iteration, FFTW3, and non-FFT DFT time comparison
  3   John Bryan 2017
  4   -lrt -lm -lfftw3 -O3   
  5   gcc 4.6.3
  6 */
  7 
  8 
  9 #include <stdio.h>
 10 #include <stdlib.h>
 11 #include <time.h>
 12 #include <fftw3.h>
 13 #include <complex.h>
 14 #include <math.h>
 15 #include <float.h>
 16 #include <assert.h>
 17 #include <stdint.h>
 18 #include <inttypes.h>
 19 
 20 
 21 #define PI acos(-1)
 22 #define MILLION 1000000
 23 #define BILLION 1000000000L
 24 #define NUMBER_OF_LENGTHS 7
 25 #define NUMBER_OF_ITERATIONS 50
 26 #define OFFSET 6
 27 
 28 
 29 double randf(void)
 30 {
 31     double low=-1.0, high=1.0;
 32     return (rand()/(double)(RAND_MAX))*abs(low-high)+low;
 33 }
 34 
 35 
 36 double get_op_time()
 37 {
 38     complex double a,b;
 39     struct timespec start,end;
 40     double ns, sum = 0.;
 41     int i;
 42     for (i=0; i<MILLION; i++)
 43     {
 44          a = randf()+randf()*I;
 45          b = randf()+randf()*I;
 46          clock_gettime(CLOCK_MONOTONIC, &start);
 47          a*b;
 48          clock_gettime(CLOCK_MONOTONIC, &end);
 49          ns = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
 50 
 51          sum = sum+ns;
 52 
 53      }
 54      return (sum/MILLION);
 55 }
 56 
 57 
 58 int n_close(double A, double B, double maxDiff, double maxRelDiff)
 59 {
 60      float diff = fabs(A - B);
 61      if (diff <= maxDiff)
 62          {return 0;}
 63      A = fabs(A);
 64      B = fabs(B);
 65      double largest = (double)(B > A) ? B : A;
 66      if (diff <= largest * maxRelDiff)
 67          {return 0;}
 68      assert(0);
 69 }
 70 
 71 
 72 void bracewell_buneman(complex double *x, int m, int N)
 73 {
 74      int M=1, i,  T, k, r[N], muplus, q, p, upper_range, nu;
 75      int nprime, rprime, rprimeprime;
 76      complex double temp;
 77      muplus = (m+1)>>1;
 78      for (i=0; i < N; i++)
 79          {r[i] = 0;}
 80      upper_range = muplus + 1;
 81      for (nu=1; nu < upper_range; nu++)
 82      {
 83          for (k=0; k<M; k++)
 84          {
 85              T = r[k] << 1;
 86              r[k] = T;
 87              r[k+M] = T + 1;
 88           }
 89           M = M << 1;
 90      }
 91      if (m & 0x01) {M = M >> 1;};
 92      for (q=1; q < M; q++)
 93      {
 94           nprime = q - M;
 95           rprimeprime = r[q] * M;
 96           for (p=0; p < r[q]; p++)
 97           {
 98               nprime = nprime + M;
 99               rprime = rprimeprime + r[p];
100               temp = x[nprime];
101               x[nprime] = x[rprime];
102               x[rprime] = temp;
103           }
104       }
105       return;
106 }
107 
108 
109 void recursion_dif (int BaseE, int N, complex double *x, complex double *twiddle, int tss)
110 {
111     int Nprime=0, BaseO=0, n=0, twiddle_factor;
112     complex double e, o;
113     if (N==1) return;
114     else
115     {
116        Nprime = N/2;
117        BaseO = BaseE + Nprime;
118        for (n=0; n < Nprime; n++)
119        {
120            twiddle_factor = n*tss;
121            e = x[BaseE+n] + x[BaseO+n];
122            o = (x[BaseE+n] - x[BaseO+n]) * twiddle[twiddle_factor];
123            x[BaseE+n] = e;
124            x[BaseO+n] = o;
125         }
126         tss = 2*tss;
127         recursion_dif(BaseE, Nprime, x, twiddle, tss);
128         recursion_dif(BaseO, Nprime, x, twiddle, tss);
129     }
130 }
131 
132 
133 void iteration_dif(int N, complex double *x2, complex double *twiddle, int p)
134 {
135     int Nprime = 0, BaseE = 0, BaseO = 0, n = 0;
136     int twiddle_factor = 0, Bp = 1, tss = 1, P=0, b=0;
137     complex double e, o;
138     for (P=0; P<p; P++)
139     {
140        Nprime = N/2;
141        BaseE = 0;
142        for (b=0; b<Bp; b++)
143        {
144            BaseO = BaseE+Nprime;
145            for (n=0; n<Nprime; n++)
146            {
147                twiddle_factor = n*tss;
148                e = x2[BaseE+n]+x2[BaseO+n];
149                o = (x2[BaseE+n]-x2[BaseO+n])*twiddle[twiddle_factor];
150                x2[BaseE+n] = e;
151                x2[BaseO+n] = o;
152             }
153             BaseE = BaseE+N;
154        }
155        Bp = Bp*2;
156        N = N/2;
157        tss = 2*tss;
158     }
159 }
160 
161 
162 int main()
163 {
164     int p;
165     int N;
166     int n, u, c, k, i, nk, tss=1, baseE=0;
167     double dft_mean[NUMBER_OF_LENGTHS], difr_mean[NUMBER_OF_LENGTHS];
168     double difi_mean[NUMBER_OF_LENGTHS], fftw3_mean[NUMBER_OF_LENGTHS];
169     uint64_t dft_mean_int[NUMBER_OF_LENGTHS], difr_mean_int[NUMBER_OF_LENGTHS];
170     uint64_t difi_mean_int[NUMBER_OF_LENGTHS], fftw3_mean_int[NUMBER_OF_LENGTHS];
171     double dft_sum, difr_sum, difi_sum, fftw3_sum;
172     double dft_diff, difr_diff, difi_diff, fftw3_diff;
173     struct timespec start, end;
174     FILE *f = fopen("file.txt","w");
175     double complex_mult_time = get_op_time();
176     double iterations_complex_mult_time_product = NUMBER_OF_ITERATIONS * complex_mult_time;
177     complex double *x, *dft, *twiddle, *complex_exponential;
178     fftw_complex *in, *out;
179     fftw_plan p2;
180     for (u=0; u<NUMBER_OF_LENGTHS; u++)
181     {
182         printf(" \n%d sequence length(s) left to compute performance...\n",NUMBER_OF_LENGTHS-u);
183         dft_sum = 0.;
184         difr_sum = 0.;
185         difi_sum = 0.;
186         fftw3_sum = 0.;
187         dft_diff = 0.;
188         difr_diff = 0.;
189         difi_diff = 0.;
190         fftw3_diff = 0.;
191         p = u + OFFSET;
192         N = (int)powf(2.0, (float)(p));
193         dft = malloc(sizeof(complex double) * N);
194         x = malloc(sizeof(complex double) * N);
195         complex_exponential = malloc(sizeof(complex double) * N);
196         twiddle = malloc(sizeof(complex double) * N/2);
197         for ( nk = 0; nk < N; nk++ )
198             complex_exponential[nk]=cexp(-2*PI*I*nk/N);
199         for (c=0; c<NUMBER_OF_ITERATIONS; c++)
200         {
201             srand(1);
202             for ( i = 0; i < N; i++ )
203                 x[i] = randf()+randf()*I;
204             for ( i = 0; i < N; i++ )
205                 dft[i] = 0.0 + 0.0*I;
206             clock_gettime(CLOCK_MONOTONIC, &start);
207             for (k=0; k<N; k++)
208                for (n=0; n<N; n++)
209                    dft[k] = dft[k] + (x[n] * complex_exponential[(n*k)%N]);
210             clock_gettime(CLOCK_MONOTONIC, &end);
211             dft_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
212             dft_sum = dft_sum + dft_diff;
213             for ( i = 0; i < N/2; i++ )
214                 twiddle[i]=cexp(-2*PI*I*i/N);
215             clock_gettime(CLOCK_MONOTONIC, &start);
216             iteration_dif(N,x,twiddle,p);
217             bracewell_buneman(x,p,N);
218             clock_gettime(CLOCK_MONOTONIC, &end);
219             difi_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
220             difi_sum = difi_sum + difi_diff;
221             for (i=0;i<N;i++)
222             {
223                 n_close(creal(x[i])/N,creal(dft[i])/N,0.01,DBL_EPSILON);
224                 n_close(cimag(x[i])/N,cimag(dft[i])/N,0.01,DBL_EPSILON);
225             }
226             srand(1);
227             for ( i = 0; i < N; i++ )
228                 x[i] = randf()+randf()*I;
229             clock_gettime(CLOCK_MONOTONIC, &start);
230             recursion_dif(baseE,N,x,twiddle,tss);
231             bracewell_buneman(x,p,N);
232             clock_gettime(CLOCK_MONOTONIC, &end);
233             difr_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
234             difr_sum = difr_sum + difr_diff;
235             for (i=0;i<N;i++)
236             {
237                 n_close(creal(x[i])/N,creal(dft[i])/N,0.01,DBL_EPSILON);
238                 n_close(cimag(x[i])/N,cimag(dft[i])/N,0.01,DBL_EPSILON);
239             }
240             in  = fftw_malloc(sizeof(fftw_complex)*N);
241             out = fftw_malloc(sizeof(fftw_complex)*N );
242             srand(1);
243             for ( i = 0; i < N; i++ )
244             {
245                 x[i] = randf()+randf()*I;
246                 in[i][0] = creal(x[i]);
247                 in[i][1] = cimag(x[i]);
248              }
249              p2 = fftw_plan_dft_1d ( N, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
250             clock_gettime(CLOCK_MONOTONIC, &start);
251             fftw_execute(p2);
252             clock_gettime(CLOCK_MONOTONIC,&end);
253             fftw3_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
254             fftw3_sum = fftw3_sum + fftw3_diff;
255             fftw_destroy_plan(p2);
256         }
257         free(x);
258         free(dft);
259         free(twiddle);
260         free(complex_exponential);
261         fftw_free(in);
262         fftw_free(out);
263         dft_mean[u] = dft_sum/iterations_complex_mult_time_product;
264         difr_mean[u] = difr_sum/iterations_complex_mult_time_product;
265         difi_mean[u] = difi_sum/iterations_complex_mult_time_product;
266         fftw3_mean[u] = fftw3_sum/iterations_complex_mult_time_product;
267         dft_mean_int[u] = (uint64_t)round(dft_mean[u]);
268         difr_mean_int[u] = (uint64_t)round(difr_mean[u]);
269         difi_mean_int[u] = (uint64_t)round(difi_mean[u]);
270         fftw3_mean_int[u] = (uint64_t)round(fftw3_mean[u]);
271 
272      }
273      for ( i = 0; i < NUMBER_OF_LENGTHS; i++ )
274         fprintf (f, "  %3d   "  "%"PRIu64  "\n", (int)powf(2.0, (float)(i+OFFSET)) , difr_mean_int[i] );
275      for ( i = 0; i< NUMBER_OF_LENGTHS; i++ )
276         fprintf (f, "  %3d   "  "%"PRIu64  "\n", (int)powf(2.0, (float)(i+OFFSET)) , difi_mean_int[i] );
277      for ( i = 0; i < NUMBER_OF_LENGTHS; i++ )
278         fprintf (f, "  %3d   "  "%"PRIu64  "\n", (int)powf(2.0, (float)(i+OFFSET)) , dft_mean_int[i] );
279      for ( i = 0; i < NUMBER_OF_LENGTHS; i++ )
280         fprintf (f, "  %3d   "  "%"PRIu64  "\n", (int)powf(2.0, (float)(i+OFFSET)) , fftw3_mean_int[i] );
281 
282      fclose(f);
283      system("python plot11.py");
284      exit(0);
285 }