added 2d dft + bmp api (TESTING - may break!)

[my-code/api.git] / fourier / fourier.c

diff --git a/fourier/fourier.c b/fourier/fourier.c
index 045bf65..d322bc0 100644 (file)
--- a/fourier/fourier.c
+++ b/fourier/fourier.c
@@ -24,6 +24,9 @@ int fourier_alloc_data(t_fourier *fourier) {
   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]);
 
@@ -46,11 +49,77 @@ int fourier_shutdown(t_fourier *fourier) {
   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;
+  }
+
+  dprintf(fourier->outfd,"[fourier] fft: allocated data pointers\n");
+
+  /* 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,j,k;
+  int i,j;
+
+  /* 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]->i);
+
+  return F_NOT_SUPPORTED;
 
   return F_SUCCESS;
 }
@@ -90,7 +159,37 @@ int fourier_dft_1d(t_fourier *fourier) {
 }
 
 int fourier_dft_2d(t_fourier *fourier) {
-  return F_NOT_SUPPORTED;
+
+  int off_f,off_r;
+  int u,v,x,y;
+  int X,Y;
+  double arg;
+
+  X=fourier->data_len[0];
+  Y=fourier->data_len[1];
+
+  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))) */
+  for(v=0;v<Y;v++) {
+    off=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);
+        }
+      }
+      fourier->ftdata[off_f+u].r/=(X*Y);
+      fourier->ftdata[off_f+u].i/=(X*Y);
+    }
+  }
+
+  return F_SUCCESS;
 }
 
 int fourier_dft_3d(t_fourier *fourier) {