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_dft_1d(t_fourier *fourier) {
  50
  51   int i,k;
  52   double arg;
  53
  54   if(fourier->type&FWD) {
  55     for(k=0;k<fourier->data_len[0];k++) {
  56       fourier->ftdata[k].r=0;
  57       fourier->ftdata[k].i=0;
  58       for(i=0;i<fourier->data_len[0];i++) {
  59         /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
  60         arg=-2.0*k*M_PI*i/fourier->data_len[0];
  61         fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
  62         fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
  63       }
  64       fourier->ftdata[k].r/=fourier->data_len[0];
  65       fourier->ftdata[k].i/=fourier->data_len[0];
  66     }
  67   }
  68   else {
  69     for(k=0;k<fourier->data_len[0];k++) {
  70       fourier->data[k].r=0;
  71       fourier->data[k].i=0;
  72       for(i=0;i<fourier->data_len[0];i++) {
  73         arg=2.0*k*M_PI*i/fourier->data_len[0];
  74         fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
  75         fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
  76       }
  77     }
  78   }
  79
  80   return F_SUCCESS;
  81 }
  82
  83 int fourier_dft_2d(t_fourier *fourier) {
  84   return F_NOT_SUPPORTED;
  85 }
  86
  87 int fourier_dft_3d(t_fourier *fourier) {
  88   return F_NOT_SUPPORTED;
  89 }