1 
  2 // FHT 
  3 // John Bryan 2017
  4 
  5 // -lrt -lm -lfftw3 
  6 
  7 #include <stdlib.h>
  8 #include <stdio.h>
  9 #include <time.h>
 10 #include <fftw3.h>
 11 #include <complex.h>
 12 #include <math.h>
 13 #include <stdint.h>
 14 #include <assert.h>
 15 #include <float.h>
 16 
 17 
 18 
 19 #define nu 1000000
 20 #define BILLION 1000000000L
 21 // rows is the number of lengths
 22 #define rows 10   
 23 // columns is the number of iterations averaged
 24 #define columns 300 
 25 
 26 
 27 double randf(void)
 28 {
 29     //return random value between -1 and 1
 30     double low=-1.0, high=1.0;
 31     return (rand()/(double)(RAND_MAX))*abs(low-high)+low;
 32 }
 33 
 34 
 35 double Avg(long int array[]){
 36     // compute average
 37     size_t c;
 38     double sum = 0;
 39     for (c = 0; c < nu; c++){ sum += array[c]; }
 40     return sum / nu;
 41 }
 42 
 43 
 44 
 45 int n_close(long double A, long double B, long double maxDiff, long double maxRelDiff)
 46 {
 47        float diff = fabs(A - B);
 48        if (diff <= maxDiff)
 49            {return 0;}
 50        A = fabs(A);
 51        B = fabs(B);
 52        long double largest = (long double)(B > A) ? B : A;
 53        if (diff <= largest * maxRelDiff)
 54            {return 0;}
 55        assert(0);
 56 }
 57 
 58 
 59 void swap (complex long double *x, int i, int j)
 60 {
 61     // used in bitreversal 
 62     complex long double temp;
 63     temp=x[i];
 64     x[i]=x[j];
 65     x[j]=temp;
 66 }
 67 
 68 
 69 void swap2 (long double *x, int i, int j)
 70 {
 71     // used in bitreversal2
 72     long double temp;
 73     temp=x[i];
 74     x[i]=x[j];
 75     x[j]=temp;
 76 }
 77 
 78 
 79 complex long double * bitreversal(complex long double *x, int N)
 80 {
 81      // used just before dit fft
 82      int M=1,T,k,n,h,r[N];
 83      for (h=0;h<N;h++)
 84          r[h]=0;
 85      while (M<N)
 86      {
 87          k=0;
 88          while (k<M)
 89          {
 90              T= (int)(2*r[k]);
 91              r[k]=T;
 92              r[k+M]=T+1;
 93              k=k+1;
 94           }
 95           M = (int)(2*M);
 96       }
 97       n=1;
 98       while (n<(N-1)) {
 99           if (r[n]>n) swap(x,n,r[n]);
100           n=n+1; }
101       return x;
102 }
103 
104 
105 long double * bitreversal2(long double *x, int N)
106 {
107      // used just before fht 
108      int M=1,T,k,n,h,r[N];
109      for (h=0;h<N;h++)
110          r[h]=0;
111      while (M<N)
112      {
113          k=0;
114          while (k<M)
115          {
116              T= (int)(2*r[k]);
117              r[k]=T;
118              r[k+M]=T+1;
119              k=k+1;
120           }
121           M = (int)(2*M);
122       }
123       n=1;
124       while (n<(N-1)) {
125           if (r[n]>n) swap2(x,n,r[n]);
126           n=n+1; }
127       return x;
128 }
129 
130 complex long double * fft(int N, complex long double *x2, complex long double *t,int p)
131 {
132     // DIT fft
133     int Np=2,Npp=1,BaseT=0,BaseB=0,n=0,tf=0,Bp=N/2,tss=N/2,P=0,b=0;
134     complex long double bot,top;
135     for (P=0; P<p; P++)
136     {
137        Npp=Np/2;
138        BaseT=0;
139        for (b=0; b<Bp; b++)
140        {
141            BaseB=BaseT+Npp;
142            for (n=0;n<Npp;n++)
143            {
144                if (n==0)
145                   bot=x2[BaseB+n];
146                else
147                {
148                   tf=n*tss;
149                   bot=x2[BaseB+n]*t[tf];
150                }
151                top=x2[BaseT+n];
152                x2[BaseT+n]=top+bot;
153                x2[BaseB+n]=top-bot;
154             }
155             BaseT=BaseT+Np;
156        }
157        Bp=Bp/2;
158        Np=Np*2;
159        tss=tss/2;
160     }
161     return x2;
162 }
163 
164 
165 long double *  fht(long double *x, long double *xo, long double *c, long double *s,int N,int p)
166 {
167     // fast hartley tranform 
168     int Np=2,Npp=1,baseT=0,baseB=0,baseBB=0,Nmn=0,n=0,tf=0,Bp=N/2,tss=N/2,P=0,b=0;
169     long double xcs;
170     for (P=0; P<p; P++)
171     {
172        Npp=Np/2;
173        baseT=0;
174        for (b=0; b<Bp; b++)
175        {
176            baseB=baseT+Npp;
177            baseBB=baseT+Np;
178            for (n=0;n<Npp;n++)
179            {
180                tf=n*tss;
181                Nmn=(baseBB-n)%baseBB;
182                if (P%2==0)
183                {
184                    if (n==0)
185                        xcs=x[baseB+n];
186                    else
187                        xcs=(x[baseB+n]*c[tf])+(x[Nmn]*s[tf]);
188                    xo[baseT+n]=x[baseT+n]+xcs;
189                    xo[baseB+n]=x[baseT+n]-xcs;
190                }
191                else
192                {
193                    if (n==0)
194                        xcs=xo[baseB+n];
195                    else
196                         xcs=(xo[baseB+n]*c[tf])+(xo[Nmn]*s[tf]);
197                    x[baseT+n]=xo[baseT+n]+xcs;
198                    x[baseB+n]=xo[baseT+n]-xcs;
199                 }
200            }
201            baseT=baseT+Np;
202        }
203        Bp=Bp/2;
204        Np=Np*2;
205        tss=tss/2;
206     }
207     if (p%2==0)
208         return x;
209     else
210         return xo;
211 }
212 
213 
214 double dAvg(int array[], size_t length){
215     // compute average
216     size_t c;
217     double sum = 0;
218     for (c = 0; c < length; c++){ sum += array[c]; }
219     return sum / (double) length;
220 }
221 
222 
223 int main()
224 {
225    int u,i,N,cc=0,z=0,flag=0;
226    int uN=rows;
227    int it=columns;
228    double mean3[uN],mean4[uN];
229    float PI=acos(-1);
230    complex long double *x2,*t;
231    long double *x,*xo,*x3,*c,*s,*check;
232    long int stime;
233    int ts, times3[uN][it], times4[uN][it];
234    struct timespec start,end;
235    FILE *f=fopen("file.txt","w");
236    fftw_complex  *in,*out;
237    fftw_plan p;
238    for (u=0;u<uN;u++)
239    {
240       for (cc=0;cc<it;cc++)
241       {
242           N=(int)powf(2.0,(float)(u+4));
243           c=malloc(sizeof(long double)*N/2);
244           s=malloc(sizeof(long double)*N/2);
245           t=malloc(sizeof(complex long double)*N/2);
246           for ( i = 0; i < N/2; i++ ) t[i]=cexp(-2*PI*I*i/N);
247           for ( i = 0; i < N/2; i++ ) c[i]=cos(2*PI*i/N);
248           for ( i = 0; i < N/2; i++ ) s[i]=sin(2*PI*i/N);
249           x=malloc(sizeof(long double)*N);
250           xo=malloc(sizeof(long double)*N);
251           x3=malloc(sizeof(long double)*N);
252           check=malloc(sizeof(long double)*N);
253           x2=malloc(sizeof(complex long double)*N);
254           for ( i = 0; i < N; i++ ) x[i]=i;
255           for ( i = 0; i < N; i++ )x2[i]=x[i] +(0.*I);
256           in=fftw_malloc(sizeof(fftw_complex)*N);
257           out=fftw_malloc(sizeof(fftw_complex)*N );
258           for ( i = 0; i < N; i++ )
259           {
260               in[i][0]=creal(x[i]);
261               in[i][1]=cimag(x[i]);
262           }
263           x2=bitreversal(x2,N);
264           clock_gettime(CLOCK_MONOTONIC,&start);
265           x2=fft(N,x2,t,u+4);
266           clock_gettime(CLOCK_MONOTONIC,&end);
267           ts = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
268           times4[u][cc]=ts;
269           p = fftw_plan_dft_1d ( N, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
270           clock_gettime(CLOCK_MONOTONIC,&start);
271           fftw_execute(p);
272           clock_gettime(CLOCK_MONOTONIC,&end);
273           ts = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
274           times3[u][cc]=ts;
275           fftw_destroy_plan(p);
276           fftw_free(in);
277           x=bitreversal2(x,N);
278           clock_gettime(CLOCK_MONOTONIC,&start);
279           x3=fht(x,xo,c,s,N,u+4);
280           clock_gettime(CLOCK_MONOTONIC,&end);
281           ts = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
282           times3[u][cc]=ts;
283            for (i=0;i<N;i++)
284            {
285               x3[i]=x3[i]/N;
286               x2[i]=x2[i]/N;
287            }
288            for (i=0;i<N;i++)
289            {
290               out[i][0]=out[i][0]/N;
291               out[i][1]=out[i][1]/N;
292               check[i]=out[i][0]-out[i][1];
293            }
294            for (i=0;i<N;i++)
295            {
296                n_close(x3[i],check[i],0.01,DBL_EPSILON);
297                n_close(creal(x2[i]),out[i][0],0.01,DBL_EPSILON);
298                n_close(cimag(x2[i]),out[i][1],0.01,DBL_EPSILON);
299            }
300       }
301    }
302         printf("\n\nfht and fft are correct\n\n");
303    fftw_free(out);
304    free(x);
305    free(x2);
306    free(xo);
307    free(t);
308    free(c);
309    free(s);
310    for (z=0; z<rows; z++)
311    {
312         mean3[z]=dAvg(times3[z], columns);
313         mean4[z]=dAvg(times4[z], columns);
314    }
315    //loop fft time results
316    for ( i = 0; i < uN; i++ )
317       fprintf (f, "  %3d   %llu \n",(int)powf(2.0,(float)(i+4)) , (long long unsigned int)mean4[i] );
318    //fht time results
319    for ( i = 0; i < uN; i++ )
320       fprintf (f, "  %3d   %llu \n",(int)powf(2.0,(float)(i+4)) , (long long unsigned int)mean3[i] );
321    fclose(f);
322    exit(0);
323 }