X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=fourier%2Ffourier.c;h=89f42c24ce8b132c9d22cbb23646dad95eeac827;hb=HEAD;hp=d322bc03ae0d807ea5e5a4da1a52c801dd704a52;hpb=f8382861ee7bff6d129bd149579d26bf9b81ada0;p=my-code%2Fapi.git

diff --git a/fourier/fourier.c b/fourier/fourier.c
index d322bc0..89f42c2 100644
--- a/fourier/fourier.c
+++ b/fourier/fourier.c
@@ -103,7 +103,7 @@ int fourier_fft_1d_shutdown(t_fourier *fourier) {
 
 int fourier_fft_1d(t_fourier *fourier) {
 
-  int i,j;
+  int i;
 
   /* copy src to destination, destination is modified in place */
   memcpy(fourier->ftdata,fourier->data,fourier->data_len[0]*sizeof(t_complex));
@@ -117,7 +117,7 @@ int fourier_fft_1d(t_fourier *fourier) {
 
   for(i=0;i<fourier->data_len[0];i++) 
     dprintf(fourier->outfd,"%f %f\n",
-            fourier->ftdata[i]->r,fourier->revdata[i]->i);
+            fourier->ftdata[i].r,fourier->revdata[i]->r);
 
   return F_NOT_SUPPORTED;
 
@@ -163,31 +163,80 @@ int fourier_dft_2d(t_fourier *fourier) {
   int off_f,off_r;
   int u,v,x,y;
   int X,Y;
+  int size;
   double arg;
+  t_complex *data;
 
   X=fourier->data_len[0];
   Y=fourier->data_len[1];
 
+  size=X*Y;
+
+  data=(t_complex *)malloc(size*sizeof(t_complex));
+  if(data==NULL) {
+    dprintf(fourier->outfd,
+            "[fourier] malloc of %d bytes of temp data failed\n",size);
+    return F_ALLOC_FAIL;
+  }
+  memset(data,0,size*sizeof(t_complex));
+
   if(fourier->type&BWD) return F_NOT_SUPPORTED;
 
-  /* stupid way - actually you would do: ft_on_index_1(ft_on_index_2(f(x,y))) */
+  /* stupid way */
+  /*
   for(v=0;v<Y;v++) {
-    off=v*X;
+    off_f=v*X;
     for(u=0;u<X;u++) {
-      fourier->ftdata[off+u].r=0;
-      fourier->ftdata[off+u].i=0;
       for(y=0;y<Y;y++) {
         off_r=y*X;
         for(x=0;x<X;x++) {
-          arg=-2.0*M_PI*((1.0*u*x)/X+(1.0*v*y)/Y);
-          fourier->ftdata[off+u].r+=(cos(arg)*fourier->data[off_r+x].r-sin(arg)*fourier->data[off_r+x].i);
-          fourier->ftdata[off+u].i+=(sin(arg)*fourier->data[off_r+x].r+cos(arg)*fourier->data[off_r+x].i);
+         arg=-2.0*M_PI*((1.0*u*x)/X+(1.0*v*y)/Y);
+          fourier->ftdata[off_f+u].r+=(cos(arg)*fourier->data[off_r+x].r-sin(arg)*fourier->data[off_r+x].i);
+          fourier->ftdata[off_f+u].i+=(sin(arg)*fourier->data[off_r+x].r+cos(arg)*fourier->data[off_r+x].i);
         }
       }
       fourier->ftdata[off_f+u].r/=(X*Y);
       fourier->ftdata[off_f+u].i/=(X*Y);
     }
   }
+  */
+
+  /* the more clever way ... */
+  // dft on index 1
+  off_f=0;
+  for(y=0;y<Y;y++) {
+    for(u=0;u<X;u++) {
+      //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,y);
+      for(x=0;x<X;x++) {
+        arg=-2.0*M_PI*u*x/X;
+        off_r=off_f+x;
+        data[off_f+u].r+=(cos(arg)*fourier->data[off_r].r-sin(arg)*fourier->data[off_r].i);
+        data[off_f+u].r+=(sin(arg)*fourier->data[off_r].r+cos(arg)*fourier->data[off_r].i);
+      }
+      data[off_f+u].r/=X;
+      data[off_f+u].i/=X;
+    }
+    off_f+=X;
+  }
+  //  dft on index 2 of 'index 1 transformed data'
+  for(x=0;x<X;x++) {
+    off_f=x;
+    for(v=0;v<Y;v++) {
+      //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
+      off_r=x;
+      for(y=0;y<Y;y++) {
+        arg=-2.0*M_PI*v*y/Y;
+        fourier->ftdata[off_f].r+=(cos(arg)*data[off_r].r-sin(arg)*data[off_r].i);
+        fourier->ftdata[off_f].r+=(sin(arg)*data[off_r].r+cos(arg)*data[off_r].i);
+        off_r+=X;
+      }
+      fourier->ftdata[off_f].r/=Y;
+      fourier->ftdata[off_f].i/=Y;
+      off_f+=X;
+    }
+  }
+ 
+  free(data);
 
   return F_SUCCESS;
 }