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;
  23
  24   for(i=0;i<fourier->dim;i++) {
  25     fourier->data[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
  26     fourier->ftdata[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
  27     if((fourier->data[i]==NULL)||(fourier->ftdata[i]==NULL)) {
  28       dprintf(fourier->outfd,"[fourier] malloc failed\n");
  29       return F_ALLOC_FAIL;
  30     }
  31   }
  32
  33   return F_SUCCESS;
  34 }
  35
  36 int fourier_shutdown(t_fourier *fourier) {
  37
  38   int i;
  39
  40   dprintf(fourier->outfd,"[fourier] shutdown\n");
  41
  42   for(i=0;i<fourier->dim;i++) {
  43     free(fourier->data[i]);
  44     free(fourier->ftdata[i]);
  45   }
  46
  47   return F_SUCCESS;
  48 }
  49
  50 int fourier_dft_1d(t_fourier *fourier) {
  51
  52   int i,k,dim;
  53   double arg;
  54
  55   if(fourier->type&FWD) {
  56     for(dim=0;dim<fourier->dim;dim++) {
  57       for(k=0;k<fourier->data_len[dim];k++) {
  58         fourier->ftdata[dim][k].r=0;
  59         fourier->ftdata[dim][k].i=0;
  60         for(i=0;i<fourier->data_len[dim];i++) {
  61           /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
  62           arg=-1.0*k*M_PI*i/fourier->data_len[dim];
  63           fourier->ftdata[dim][k].r+=(cos(arg)*fourier->data[dim][i].r-sin(arg)*fourier->data[dim][i].i);
  64           fourier->ftdata[dim][k].i+=(sin(arg)*fourier->data[dim][i].r+cos(arg)*fourier->data[dim][i].i);
  65         }
  66         fourier->ftdata[dim][k].r/=fourier->data_len[dim];
  67         fourier->ftdata[dim][k].i/=fourier->data_len[dim];
  68       }
  69     }
  70   }
  71   else {
  72     for(dim=0;dim<fourier->dim;dim++) {
  73       for(k=0;k<fourier->data_len[dim];k++) {
  74         fourier->data[dim][k].r=0;
  75         fourier->data[dim][k].i=0;
  76         for(i=0;i<fourier->data_len[dim];i++) {
  77           arg=1.0*k*M_PI*i/fourier->data_len[dim];
  78           fourier->data[dim][k].r+=(cos(arg)*fourier->ftdata[dim][i].r-sin(arg)*fourier->ftdata[dim][i].i);
  79           fourier->data[dim][k].i+=(sin(arg)*fourier->ftdata[dim][i].r+cos(arg)*fourier->ftdata[dim][i].i);
  80         }
  81       }
  82     }
  83   }
  84
  85   return F_SUCCESS;
  86 }
  87
  88 int fourier_dft_2d(t_fourier *fourier) {
  89   return 0;
  90 }
  91
  92 int fourier_dft_3d(t_fourier *fourier) {
  93   return 0;
  94 }
  95
  96 int fourier_calc(t_fourier *fourier) {
  97
  98   dprintf(fourier->outfd,"[fourier] %d dimensional %c-%cft calculation ...\n",
  99           fourier->dim,
 100           ((fourier->type&FWD)?'f':'b'),
 101           ((fourier->type&DFT)?'d':'f'));
 102
 103   if(fourier->type&DFT) {
 104     switch(fourier->dim) {
 105       case 1:
 106         fourier_dft_1d(fourier);
 107         break;
 108       case 2:
 109         fourier_dft_2d(fourier);
 110         break;
 111       case 3:
 112         fourier_dft_3d(fourier);
 113         break;
 114       default:
 115         dprintf(fourier->outfd,"[fourier] dimension failure\n");
 116         return F_DIM_FAILURE;
 117     }
 118     return F_SUCCESS;
 119   }
 120   else {
 121     dprintf(fourier->outfd,"[fourier] fft not supported by now\n");
 122     return F_NOT_SUPPORTED;
 123   }
 124 }