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   dprintf(fourier->outfd,"[fourier] allocating %d bytes for data ...\n",
  28           2*size*sizeof(t_complex));
  29
  30   fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
  31   fourier->ftdata=&(fourier->data[size]);
  32
  33   memset(fourier->data,0,2*size*sizeof(t_complex));
  34
  35   if(fourier->data==NULL) {
  36     dprintf(fourier->outfd,"[fourier] malloc failed\n");
  37     return F_ALLOC_FAIL;
  38   }
  39
  40   return F_SUCCESS;
  41 }
  42
  43 int fourier_shutdown(t_fourier *fourier) {
  44
  45   dprintf(fourier->outfd,"[fourier] shutdown\n");
  46
  47   free(fourier->data);
  48
  49   return F_SUCCESS;
  50 }
  51
  52 /* fft functions */
  53
  54 int fourier_fft_1d_init(t_fourier *fourier) {
  55
  56   unsigned int i,j,m,r;
  57   int mask[32];
  58
  59   /* calculate m = log 2 (N) + check N = 2^m */
  60   m=0;
  61   while(m<32) {
  62     if(1<<m==fourier->data_len[0]) break;
  63     else m++;
  64   }
  65   if(m==32) {
  66     dprintf(fourier->outfd,
  67             "[fourier] fft impossible (N > 2^32 or not a power of 2\n");
  68     return F_FFT_IMPOSSIBLE;
  69   }
  70   dprintf(fourier->outfd,"[fourier] log2(%d) = %d\n",fourier->data_len[0],m);
  71   fourier->log2len[0]=m;
  72
  73   /* pointer array which will point to reversed t_complex data */
  74   fourier->revdata=(t_complex **)malloc(fourier->data_len[0]*sizeof(t_complex *));
  75   if(fourier->revdata==NULL) {
  76     dprintf(fourier->outfd,
  77             "[fourier] malloc for reversed data pointers failed\n");
  78     return F_ALLOC_FAIL;
  79   }
  80
  81   dprintf(fourier->outfd,"[fourier] fft: allocated data pointers\n");
  82
  83   /* swap data (bit reversal) */
  84   for(i=0;i<m;i++) mask[i]=1<<i;
  85   for(i=0;i<fourier->data_len[0];i++) {
  86     r=0;
  87     for(j=0;j<m;j++) r+=(((i&mask[j])>>j)<<(m-j-1));
  88     fourier->revdata[i]=fourier->ftdata+r;
  89   }
  90
  91   return F_SUCCESS;
  92 }
  93
  94 int fourier_fft_1d_shutdown(t_fourier *fourier) {
  95
  96   dprintf(fourier->outfd,"[fourier] fft shutdown\n");
  97
  98   free(fourier->revdata);
  99
 100   return F_SUCCESS;
 101 }
 102
 103
 104 int fourier_fft_1d(t_fourier *fourier) {
 105
 106   int i;
 107
 108   /* copy src to destination, destination is modified in place */
 109   memcpy(fourier->ftdata,fourier->data,fourier->data_len[0]*sizeof(t_complex));
 110   /* from now on access bit reversed data by revdata[i]->... */
 111
 112   // for(i=0;i<fourier->log2len[0];i++) {
 113
 114
 115   // }
 116
 117
 118   for(i=0;i<fourier->data_len[0];i++)
 119     dprintf(fourier->outfd,"%f %f\n",
 120             fourier->ftdata[i].r,fourier->revdata[i]->r);
 121
 122   return F_NOT_SUPPORTED;
 123
 124   return F_SUCCESS;
 125 }
 126
 127 int fourier_dft_1d(t_fourier *fourier) {
 128
 129   int i,k;
 130   double arg;
 131
 132   if(fourier->type&FWD) {
 133     for(k=0;k<fourier->data_len[0];k++) {
 134       fourier->ftdata[k].r=0;
 135       fourier->ftdata[k].i=0;
 136       for(i=0;i<fourier->data_len[0];i++) {
 137         /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
 138         arg=-2.0*k*M_PI*i/fourier->data_len[0];
 139         fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
 140         fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
 141       }
 142       fourier->ftdata[k].r/=fourier->data_len[0];
 143       fourier->ftdata[k].i/=fourier->data_len[0];
 144     }
 145   }
 146   else {
 147     for(k=0;k<fourier->data_len[0];k++) {
 148       fourier->data[k].r=0;
 149       fourier->data[k].i=0;
 150       for(i=0;i<fourier->data_len[0];i++) {
 151         arg=2.0*k*M_PI*i/fourier->data_len[0];
 152         fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
 153         fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
 154       }
 155     }
 156   }
 157
 158   return F_SUCCESS;
 159 }
 160
 161 int fourier_dft_2d(t_fourier *fourier) {
 162
 163   int off_f,off_r;
 164   int u,v,x,y;
 165   int X,Y;
 166   int size;
 167   double arg;
 168   t_complex *data;
 169
 170   X=fourier->data_len[0];
 171   Y=fourier->data_len[1];
 172
 173   size=X*Y;
 174
 175   data=(t_complex *)malloc(size*sizeof(t_complex));
 176   if(data==NULL) {
 177     dprintf(fourier->outfd,
 178             "[fourier] malloc of %d bytes of temp data failed\n",size);
 179     return F_ALLOC_FAIL;
 180   }
 181   memset(data,0,size*sizeof(t_complex));
 182
 183   if(fourier->type&BWD) return F_NOT_SUPPORTED;
 184
 185   /* stupid way */
 186   /*
 187   for(v=0;v<Y;v++) {
 188     off_f=v*X;
 189     for(u=0;u<X;u++) {
 190       for(y=0;y<Y;y++) {
 191         off_r=y*X;
 192         for(x=0;x<X;x++) {
 193          arg=-2.0*M_PI*((1.0*u*x)/X+(1.0*v*y)/Y);
 194           fourier->ftdata[off_f+u].r+=(cos(arg)*fourier->data[off_r+x].r-sin(arg)*fourier->data[off_r+x].i);
 195           fourier->ftdata[off_f+u].i+=(sin(arg)*fourier->data[off_r+x].r+cos(arg)*fourier->data[off_r+x].i);
 196         }
 197       }
 198       fourier->ftdata[off_f+u].r/=(X*Y);
 199       fourier->ftdata[off_f+u].i/=(X*Y);
 200     }
 201   }
 202   */
 203
 204   /* the more clever way ... */
 205   // dft on index 1
 206   off_f=0;
 207   for(y=0;y<Y;y++) {
 208     for(u=0;u<X;u++) {
 209       //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,y);
 210       for(x=0;x<X;x++) {
 211         arg=-2.0*M_PI*u*x/X;
 212         off_r=off_f+x;
 213         data[off_f+u].r+=(cos(arg)*fourier->data[off_r].r-sin(arg)*fourier->data[off_r].i);
 214         data[off_f+u].r+=(sin(arg)*fourier->data[off_r].r+cos(arg)*fourier->data[off_r].i);
 215       }
 216       data[off_f+u].r/=X;
 217       data[off_f+u].i/=X;
 218     }
 219     off_f+=X;
 220   }
 221   //  dft on index 2 of 'index 1 transformed data'
 222   for(x=0;x<X;x++) {
 223     off_f=0;
 224     for(v=0;v<Y;v++) {
 225       //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
 226       off_r=0;
 227       for(y=0;y<Y;y++) {
 228         arg=-2.0*M_PI*v*y/Y;
 229         fourier->ftdata[off_f+x].r+=(cos(arg)*data[off_r+x].r-sin(arg)*data[off_r+x].i);
 230         fourier->ftdata[off_f+x].r+=(sin(arg)*data[off_r+x].r+cos(arg)*data[off_r+x].i);
 231         off_r+=X;
 232       }
 233       fourier->ftdata[off_f+x].r/=Y;
 234       fourier->ftdata[off_f+x].i/=Y;
 235       off_f+=X;
 236     }
 237   }
 238
 239   free(data);
 240
 241   return F_SUCCESS;
 242 }
 243
 244 int fourier_dft_3d(t_fourier *fourier) {
 245   return F_NOT_SUPPORTED;
 246 }