From 81e8529f1b659c27e65c66d12cb122ad19fcfba3 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Thu, 23 Sep 2004 10:15:04 +0000
Subject: [PATCH] 1d fourier trafo function added

---
 fourier/fourier.c | 97 +++++++++++++++++++++++++++++++++++++----------
 fourier/fourier.h | 20 +++++++---
 2 files changed, 92 insertions(+), 25 deletions(-)

diff --git a/fourier/fourier.c b/fourier/fourier.c
index 14878e2..fca6d6a 100644
--- a/fourier/fourier.c
+++ b/fourier/fourier.c
@@ -17,41 +17,100 @@ int fourier_init(t_fourier *fourier,int outfd) {
  return F_SUCCESS;
 }
 
-int fourier_set_datalen(t_fourier *fourier,int dim,int len) {
+int fourier_alloc_data(t_fourier *fourier) {
 
   int i;
 
-  if(dim>MAX_DIM) {
-    fprintf(fourier->outfd,"[fourier] dimension %d not supported\n",dim);
-    return F_SET_DLEN_ERROR;
+  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;
+    }
   }
 
-  fourier->data_len[dim]=len;
-
   return F_SUCCESS;
 }
 
-int fourier_malloc(t_fourier *fourier) {
+int fourier_shutdown(t_fourier *fourier) {
 
   int i;
 
-  if(fourier->dim==0) {
-    fprintf(fourier->outfd,"[fourier] dimension not specified\n");
-    return F_DIM_ERROR;
-  }
+  fprintf(fourier->outfd,"[fourier] shutdown\n");
+
   for(i=0;i<fourier->dim;i++) {
-    if(fourier->data_len[i]==0) {
-      fprintf(fourier->outfd,"[fourier] invalid data len for dim %d\n",i);
-      return F_DLEN_ERROR;
-    }
+    free(fourier->data[dim]);
+    free(fourier->ftdata[dim]);
   }
 
-  for(i=0;i<fourier->dim;i++) {
-    if((fourier->data[i]=(t_complex *)malloc(sizeof(t_complex)*fourier->data_len[i]))==NULL) {
-      fprintf(fourier->outfd,"[fourier] malloc for data dim %d failed\n",i);
-      return F_MALLOC_ERROR;
+  return F_SUCCESS;
+}
+
+int fourier_dft_1d(t_fourier *fourier) {
+
+  int i,k,dim;
+  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];
+      }
+    }
+  }
+  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);
+        }
+      }
     }
   }
 
   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;
+    }
+    return F_SUCCESS;
+  }
+  else {
+    fprintf(fourier->outfd,"[fourier] fft not supported by now\n");
+    return F_NOT_SUPPORTED;
+  }
+} 
diff --git a/fourier/fourier.h b/fourier/fourier.h
index 5a8ddb3..f3cab6b 100644
--- a/fourier/fourier.h
+++ b/fourier/fourier.h
@@ -9,12 +9,9 @@
 /* defines */
 #define F_SUCCESS 1
 #define F_ERROR -1
-#define F_DIM_ERROR -2
-#define F_DLEN_ERROR -3
-#define F_MALLOC_ERROR -4
-#define F_SET_DLEN_ERROR -5
-
-#define MAX_DIM 3
+#define F_NOT_SUPPORTED -2
+#define F_DIM_FAILURE -3
+#define F_ALLOC_FAIL -4
 
 /* fourier specific variables */
 typedef struct s_complex {
@@ -25,9 +22,20 @@ typedef struct s_complex {
 typedef struct s_fourier {
   int outfd;
   unsigned char type;
+#define DFT (1<<0)
+#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 data_len[MAX_DIM];
 } t_fourier;
 
+/* function prototypes */
+int fourier_init(t_fourier *fourier,int outfd);
+int fourier_dft_1d(t_fourier *fourier);
+int fourier_calc(t_fourier *fourier);
+
 #endif
-- 
2.20.1