fourier/fourier.c

   1 /* fourier.c -- fourier transformation api
   2  *
   3  * author: hackbard@hackdaworld.dyndns.org
   4  *
   5  */
   6
   7 #include "fourier.h"
   8
   9 int fourier_init(t_fourier *fourier,int outfd) {
  10
  11  dprintf(outfd,"[fourier] initializing fourier api ...\n");
  12
  13  fourier->outfd=outfd;
  14  fourier->type=0;
  15  fourier->dim=0;
  16
  17  return F_SUCCESS;
  18 }
  19
  20 int fourier_alloc_data(t_fourier *fourier) {
  21
  22   int i,size;
  23
  24   size=1;
  25   for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
  26
  27   fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
  28   fourier->ftdata=&(fourier->data[size]);
  29
  30   memset(fourier->data,0,2*size*sizeof(t_complex));
  31
  32   if(fourier->data==NULL) {
  33     dprintf(fourier->outfd,"[fourier] malloc failed\n");
  34     return F_ALLOC_FAIL;
  35   }
  36
  37   return F_SUCCESS;
  38 }
  39
  40 int fourier_shutdown(t_fourier *fourier) {
  41
  42   dprintf(fourier->outfd,"[fourier] shutdown\n");
  43
  44   free(fourier->data);
  45
  46   return F_SUCCESS;
  47 }
  48
  49 int fourier_fft_1d(t_fourier *fourier) {
  50
  51   int i,j,k;
  52
  53
  54
  55   return F_SUCCESS;
  56 }
  57
  58 int fourier_dft_1d(t_fourier *fourier) {
  59
  60   int i,k;
  61   double arg;
  62
  63   if(fourier->type&FWD) {
  64     for(k=0;k<fourier->data_len[0];k++) {
  65       fourier->ftdata[k].r=0;
  66       fourier->ftdata[k].i=0;
  67       for(i=0;i<fourier->data_len[0];i++) {
  68         /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
  69         arg=-2.0*k*M_PI*i/fourier->data_len[0];
  70         fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
  71         fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
  72       }
  73       fourier->ftdata[k].r/=fourier->data_len[0];
  74       fourier->ftdata[k].i/=fourier->data_len[0];
  75     }
  76   }
  77   else {
  78     for(k=0;k<fourier->data_len[0];k++) {
  79       fourier->data[k].r=0;
  80       fourier->data[k].i=0;
  81       for(i=0;i<fourier->data_len[0];i++) {
  82         arg=2.0*k*M_PI*i/fourier->data_len[0];
  83         fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
  84         fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
  85       }
  86     }
  87   }
  88
  89   return F_SUCCESS;
  90 }
  91
  92 int fourier_dft_2d(t_fourier *fourier) {
  93   return F_NOT_SUPPORTED;
  94 }
  95
  96 int fourier_dft_3d(t_fourier *fourier) {
  97   return F_NOT_SUPPORTED;
  98 }