X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=fourier%2Ffourier.c;h=9f810ecef37084d7d0bd4e17291688236b611752;hb=0e73e99a5415a252a15883b7f3594411781d2e9e;hp=fca6d6a2fdda3072aad543e6c85165102925c0b1;hpb=81e8529f1b659c27e65c66d12cb122ad19fcfba3;p=my-code%2Fapi.git

diff --git a/fourier/fourier.c b/fourier/fourier.c
index fca6d6a..9f810ec 100644
--- a/fourier/fourier.c
+++ b/fourier/fourier.c
@@ -8,7 +8,7 @@
 
 int fourier_init(t_fourier *fourier,int outfd) {
 
- fprintf(outfd,"[fourier] initializing fourier api ...\n");
+ dprintf(outfd,"[fourier] initializing fourier api ...\n");
 
  fourier->outfd=outfd;
  fourier->type=0;
@@ -19,15 +19,22 @@ int fourier_init(t_fourier *fourier,int outfd) {
 
 int fourier_alloc_data(t_fourier *fourier) {
 
-  int i;
+  int i,size;
 
-  for(i=0;i<fourier->dim;i++) {
-    fourier->data[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
-    fourier->ftdata[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
-    if((fourier->data[i]==NULL)||(fourier->ftdata[i]==NULL)) {
-      fprintf(fourier->outfd,"[fourier] malloc failed\n");
-      return F_ALLOC_FAIL;
-    }
+  size=1;
+  for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
+
+  dprintf(fourier->outfd,"[fourier] allocating %d bytes for data ...\n",
+          2*size*sizeof(t_complex));
+
+  fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
+  fourier->ftdata=&(fourier->data[size]);
+
+  memset(fourier->data,0,2*size*sizeof(t_complex));
+
+  if(fourier->data==NULL) {
+    dprintf(fourier->outfd,"[fourier] malloc failed\n");
+    return F_ALLOC_FAIL;
   }
 
   return F_SUCCESS;
@@ -35,49 +42,115 @@ int fourier_alloc_data(t_fourier *fourier) {
 
 int fourier_shutdown(t_fourier *fourier) {
 
-  int i;
+  dprintf(fourier->outfd,"[fourier] shutdown\n");
+
+  free(fourier->data);
+
+  return F_SUCCESS;
+}
+
+/* fft functions */
+
+int fourier_fft_1d_init(t_fourier *fourier) {
+
+  unsigned int i,j,m,r;
+  int mask[32];
+
+  /* calculate m = log 2 (N) + check N = 2^m */
+  m=0;
+  while(m<32) {
+    if(1<<m==fourier->data_len[0]) break;
+    else m++;
+  }
+  if(m==32) {
+    dprintf(fourier->outfd,
+            "[fourier] fft impossible (N > 2^32 or not a power of 2\n");
+    return F_FFT_IMPOSSIBLE;
+  }
+  dprintf(fourier->outfd,"[fourier] log2(%d) = %d\n",fourier->data_len[0],m);
+  fourier->log2len[0]=m;
+
+  /* pointer array which will point to reversed t_complex data */
+  fourier->revdata=(t_complex **)malloc(fourier->data_len[0]*sizeof(t_complex *));
+  if(fourier->revdata==NULL) {
+    dprintf(fourier->outfd,
+            "[fourier] malloc for reversed data pointers failed\n");
+    return F_ALLOC_FAIL;
+  }
 
-  fprintf(fourier->outfd,"[fourier] shutdown\n");
+  dprintf(fourier->outfd,"[fourier] fft: allocated data pointers\n");
 
-  for(i=0;i<fourier->dim;i++) {
-    free(fourier->data[dim]);
-    free(fourier->ftdata[dim]);
+  /* swap data (bit reversal) */
+  for(i=0;i<m;i++) mask[i]=1<<i;
+  for(i=0;i<fourier->data_len[0];i++) {
+    r=0;
+    for(j=0;j<m;j++) r+=(((i&mask[j])>>j)<<(m-j-1));
+    fourier->revdata[i]=fourier->ftdata+r;
   }
 
   return F_SUCCESS;
 }
 
+int fourier_fft_1d_shutdown(t_fourier *fourier) {
+
+  dprintf(fourier->outfd,"[fourier] fft shutdown\n");
+
+  free(fourier->revdata);
+
+  return F_SUCCESS;
+}
+
+
+int fourier_fft_1d(t_fourier *fourier) {
+
+  int i;
+
+  /* copy src to destination, destination is modified in place */
+  memcpy(fourier->ftdata,fourier->data,fourier->data_len[0]*sizeof(t_complex));
+  /* from now on access bit reversed data by revdata[i]->... */
+
+  // for(i=0;i<fourier->log2len[0];i++) {
+
+
+  // }
+
+
+  for(i=0;i<fourier->data_len[0];i++) 
+    dprintf(fourier->outfd,"%f %f\n",
+            fourier->ftdata[i].r,fourier->revdata[i]->r);
+
+  return F_NOT_SUPPORTED;
+
+  return F_SUCCESS;
+}
+
 int fourier_dft_1d(t_fourier *fourier) {
 
-  int i,k,dim;
+  int i,k;
   double arg;
 
   if(fourier->type&FWD) {
-    for(dim=1;dim<=fourier->dim;dim++) {
-      for(k=0;k<fourier->data_len[dim];k++) {
-        fourier->ftdata[dim][k].r=0;
-        fourier->ftdata[dim][k].i=0;
-        for(i=0;i<fourier->data_len[dim];i++) {
-    	  /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
-          arg=-1.0*k*M_PI*i/fourier->data_len[dim];
-          fourier->ftdata[dim][k].r+=(cos(arg)*fourier->data[dim][i]-sin(arg)*fourier->data[dim][i].i);
-          fourier->ftdata[dim][k].i+=(sin(arg)*fourier->data[dim][i]+cos(arg)*fourier->data[dim][i].i);
-        }
-        fourier->ftdata[dim][k].r/=fourier->data_len[dim];
-        fourier->ftdata[dim][k].i/=fourier->data_len[dim];
+    for(k=0;k<fourier->data_len[0];k++) {
+      fourier->ftdata[k].r=0;
+      fourier->ftdata[k].i=0;
+      for(i=0;i<fourier->data_len[0];i++) {
+        /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
+        arg=-2.0*k*M_PI*i/fourier->data_len[0];
+        fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
+        fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
       }
+      fourier->ftdata[k].r/=fourier->data_len[0];
+      fourier->ftdata[k].i/=fourier->data_len[0];
     }
   }
   else {
-    for(dim=1;dim<=fourier->dim;dim++) {
-      for(k=0;k<fourier->data_len[dim];k++) {
-        fourier->data[dim][k].r=0;
-        fourier->data[dim][k].i=0;
-        for(i=0;i<fourier->data_len[dim];i++) {
-          arg=1.0*k*M_PI*i/fourier->data_len[dim];
-          fourier->data[dim][k].r+=(cos(arg)*fourier->ftdata[dim][i]-sin(arg)*fourier->ftdata[dim][i].i);
-          fourier->data[dim][k].i+=(sin(arg)*fourier->ftdata[dim][i]+cos(arg)*fourier->ftdata[dim][i].i);
-        }
+    for(k=0;k<fourier->data_len[0];k++) {
+      fourier->data[k].r=0;
+      fourier->data[k].i=0;
+      for(i=0;i<fourier->data_len[0];i++) {
+        arg=2.0*k*M_PI*i/fourier->data_len[0];
+        fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
+        fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
       }
     }
   }
@@ -85,32 +158,89 @@ int fourier_dft_1d(t_fourier *fourier) {
   return F_SUCCESS;
 }
 
-int fourier_calc(t_fourier *fourier) {
-
-  fprintf(fourier->outfd,"[fourier] %d dimensional %c-%cft calculation ...\n",
-          fourier->dim,
-          ((fourier->type&FWD)?'f':'b'),
-          ((fourier->type&DFT)?'d':'f'));
-
-  if(fourier->type&DFT) {
-    switch(fourier->dim) {
-      case 1:
-        fourier_dft_1d(fourier);
-        break;
-      case 2:
-        fourier_dft_2d(fourier);
-        break;
-      case 3:
-        fourier_dft_3d(fourier);
-        break;
-      default:
-        fprintf(fourier->outfd,"[fourier] dimension failure\n");
-        return F_DIM_FAILURE;
+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 */
+  /*
+  for(v=0;v<Y;v++) {
+    off_f=v*X;
+    for(u=0;u<X;u++) {
+      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_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);
     }
-    return F_SUCCESS;
   }
-  else {
-    fprintf(fourier->outfd,"[fourier] fft not supported by now\n");
-    return F_NOT_SUPPORTED;
+  */
+
+  /* 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=0;
+    for(v=0;v<Y;v++) {
+      //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
+      off_r=0;
+      for(y=0;y<Y;y++) {
+        arg=-2.0*M_PI*v*y/Y;
+        fourier->ftdata[off_f+x].r+=(cos(arg)*data[off_r+x].r-sin(arg)*data[off_r+x].i);
+        fourier->ftdata[off_f+x].r+=(sin(arg)*data[off_r+x].r+cos(arg)*data[off_r+x].i);
+        off_r+=X;
+      }
+      fourier->ftdata[off_f+x].r/=Y;
+      fourier->ftdata[off_f+x].i/=Y;
+      off_f+=X;
+    }
+  }
+
+  free(data);
+
+  return F_SUCCESS;
+}
+
+int fourier_dft_3d(t_fourier *fourier) {
+  return F_NOT_SUPPORTED;
+}