From a005d5248c2634e62b181d8a7be870863e9a11b0 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Wed, 14 Sep 2005 17:31:46 +0000
Subject: [PATCH] fftw3 support

---
 dft.c    | 64 ++++++++++++++++++++++++++++++++++++++++++++------------
 dft.h    |  8 ++++++-
 makefile |  2 +-
 3 files changed, 59 insertions(+), 15 deletions(-)

diff --git a/dft.c b/dft.c
index 1d26cee..848e6e5 100644
--- a/dft.c
+++ b/dft.c
@@ -5,17 +5,25 @@
  *
  */
 
+//#define USE_FFTW3
+
 #include "dft.h"
 
 int main(int argc,char **argv) {
 
+#ifdef USE_FFTW3
+  fftw_plan plan;
+  fftw_complex *in,*out;
+#else
   t_fourier fourier;
+#endif
   t_bmp src;
   t_bmp dst;
   t_bmp cut;
   int x,y;
   int offy,offt;
-  int scale;
+  double scale;
+  int size;
   double max;
   double *mag;
 
@@ -27,9 +35,11 @@ int main(int argc,char **argv) {
   strcpy(src.file,argv[1]);
   strcpy(dst.file,argv[2]);
 
+#ifndef USE_FFTW3
   fourier_init(&fourier,1);
   fourier.type=DFT|FWD;
   fourier.dim=2;
+#endif
 
   bmp_read_file(&src);
 
@@ -39,37 +49,60 @@ int main(int argc,char **argv) {
   dst.height=cut.height;
   bmp_alloc_map(&dst);
 
+  size=cut.width*cut.height;
+
+#ifndef USE_FFTW3
   fourier.data_len[0]=cut.width;
   fourier.data_len[1]=cut.height;
   fourier_alloc_data(&fourier);
+#else
+  in=fftw_malloc(sizeof(fftw_complex)*size);
+  memset(in,0,size*sizeof(fftw_complex));
+  out=fftw_malloc(sizeof(fftw_complex)*size);
+  memset(out,0,size*sizeof(fftw_complex));
+  if(in==NULL||out==NULL) {
+    printf("fftw alloc error!\n");
+    return -23;
+  }
+  plan=fftw_plan_dft_2d(cut.width,cut.height,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
+#endif
 
-  mag=(double *)malloc(fourier.data_len[0]*fourier.data_len[1]*sizeof(double));
+  mag=(double *)malloc(size*sizeof(double));
   if(mag==NULL) {
     printf("unable to alloc mag memory\n");
     return -1;
   }
 
   // copy the data
-  offy=0;
-  for(y=0;y<fourier.data_len[1];y++) {
-    for(x=0;x<fourier.data_len[0];x++) {
-      offt=offy+x;
-      fourier.data[offt].r=(src.map[offt].r+src.map[offt].g+src.map[offt].b)/3;
-    }
-    offy+=fourier.data_len[0];
+  for(x=0;x<size;x++) {
+#ifdef USE_FFTW3
+    in[x][0]=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
+    in[x][1]=0;
+#else
+    fourier.data[x].r=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
+#endif
   }
 
   // do the fourier trafo
+#ifdef USE_FFTW3
+  fftw_execute(plan);
+#else
   fourier_dft_2d(&fourier);
+#endif
+
 
   printf("fourier done!\n");
 
   // copy back the data, intensity = sqrt(r^2+i^2)
   max=0;
   offt=0;
-  for(y=0;y<fourier.data_len[1];y++) {
-    for(x=0;x<fourier.data_len[0];x++) {
+  for(y=0;y<cut.height;y++) {
+    for(x=0;x<cut.width;x++) {
+#ifdef USE_FFTW3
+      mag[offt]=sqrt(out[offt][0]*out[offt][0]+out[offt][1]*out[offt][1]);
+#else
       mag[offt]=sqrt(fourier.ftdata[offt].r*fourier.ftdata[offt].r+fourier.ftdata[offt].i*fourier.ftdata[offt].i);
+#endif
       if(mag[offt]>max) max=mag[offt];
       offt+=1;
     }
@@ -77,11 +110,11 @@ int main(int argc,char **argv) {
 
   printf("found max: %f\n",max);
 
-  if(max!=0) scale=(int)255/max;
+  if(max!=0) scale=255.0/max;
   else scale=0;
 
   if(scale!=0) {
-    printf("scaling image intensity: %d\n",scale);
+    printf("scaling image intensity: %f\n",scale);
     offt=0;
     for(y=0;y<dst.height;y++) {
       for(x=0;x<dst.width;x++) {
@@ -102,7 +135,12 @@ int main(int argc,char **argv) {
 
   bmp_write_file(&dst);
 
+#ifdef USE_FFTW3
+  fftw_destroy_plan(plan);
+  fftw_free(in); fftw_free(out);
+#else
   fourier_shutdown(&fourier);
+#endif
 
   bmp_shutdown(&src);
   bmp_shutdown(&dst);
diff --git a/dft.h b/dft.h
index fe0e66c..5b78dc8 100644
--- a/dft.h
+++ b/dft.h
@@ -10,5 +10,11 @@
 #include <string.h>
 #include <string.h>
 
-#include "fourier.h"
+#ifdef USE_FFTW3
+  #include <fftw3.h>
+  #include <math.h>
+#else
+  #include "fourier.h"
+#endif
+
 #include "bmp.h"
diff --git a/makefile b/makefile
index 3f5b43d..9ecc4f4 100644
--- a/makefile
+++ b/makefile
@@ -3,7 +3,7 @@
 INCLUDEDIR = /usr/include
 
 CFLAGS = -O3 -Wall
-LIBS = -lncurses -lm
+LIBS = -lncurses -lm -lfftw3
 
 OBJS = network.o event.o input.o display.o audio.o fourier.o bmp.o
 
-- 
2.20.1