improvements
authorhackbard <hackbard>
Fri, 24 Sep 2004 10:47:32 +0000 (10:47 +0000)
committerhackbard <hackbard>
Fri, 24 Sep 2004 10:47:32 +0000 (10:47 +0000)
fourier/fourier.c
fourier/fourier.h

index d6c78f6..e877cd1 100644 (file)
@@ -19,15 +19,19 @@ int fourier_init(t_fourier *fourier,int outfd) {
 
 int fourier_alloc_data(t_fourier *fourier) {
 
-  int i;
-
-  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)) {
-      dprintf(fourier->outfd,"[fourier] malloc failed\n");
-      return F_ALLOC_FAIL;
-    }
+  int i,size;
+
+  size=1;
+  for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
+
+  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 +39,40 @@ int fourier_alloc_data(t_fourier *fourier) {
 
 int fourier_shutdown(t_fourier *fourier) {
 
-  int i;
-
   dprintf(fourier->outfd,"[fourier] shutdown\n");
 
-  for(i=0;i<fourier->dim;i++) {
-    free(fourier->data[i]);
-    free(fourier->ftdata[i]);
-  }
+  free(fourier->data);
 
   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=0;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].r-sin(arg)*fourier->data[dim][i].i);
-          fourier->ftdata[dim][k].i+=(sin(arg)*fourier->data[dim][i].r+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=0;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].r-sin(arg)*fourier->ftdata[dim][i].i);
-          fourier->data[dim][k].i+=(sin(arg)*fourier->ftdata[dim][i].r+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);
       }
     }
   }
@@ -86,39 +81,9 @@ int fourier_dft_1d(t_fourier *fourier) {
 }
 
 int fourier_dft_2d(t_fourier *fourier) {
-  return 0;
+  return F_NOT_SUPPORTED;
 }
 
 int fourier_dft_3d(t_fourier *fourier) {
-  return 0;
+  return F_NOT_SUPPORTED;
 }
-
-int fourier_calc(t_fourier *fourier) {
-
-  dprintf(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:
-        dprintf(fourier->outfd,"[fourier] dimension failure\n");
-        return F_DIM_FAILURE;
-    }
-    return F_SUCCESS;
-  }
-  else {
-    dprintf(fourier->outfd,"[fourier] fft not supported by now\n");
-    return F_NOT_SUPPORTED;
-  }
-} 
index 39aee41..fcaa3e9 100644 (file)
@@ -6,6 +6,7 @@
 /* includes */
 #define _GNU_SOURCE
 #include <stdio.h>
+#include <string.h>
 #include <stdlib.h>
 #include <math.h>
 
@@ -29,18 +30,18 @@ typedef struct s_fourier {
 #define FFT (1<<1)
 #define FWD (1<<2)
 #define BWD (1<<3)
-  int dim;
 #define MAX_DIM 3
-  t_complex *data[MAX_DIM];
-  t_complex *ftdata[MAX_DIM];
+  int dim;
+  t_complex *data;
+  t_complex *ftdata;
   int data_len[MAX_DIM];
 } t_fourier;
 
 /* function prototypes */
 int fourier_init(t_fourier *fourier,int outfd);
+int fourier_alloc_data(t_fourier *fourier);
 int fourier_dft_1d(t_fourier *fourier);
 int fourier_dft_2d(t_fourier *fourier);
 int fourier_dft_3d(t_fourier *fourier);
-int fourier_calc(t_fourier *fourier);
 
 #endif