+
+ int off_f,off_r;
+ int u,v,x,y;
+ int X,Y;
+ int size;
+ double arg;
+ t_complex *data;
+
+ X=fourier->data_len[0];
+ Y=fourier->data_len[1];
+
+ size=X*Y;
+
+ data=(t_complex *)malloc(size*sizeof(t_complex));
+ if(data==NULL) {
+ dprintf(fourier->outfd,
+ "[fourier] malloc of %d bytes of temp data failed\n",size);
+ return F_ALLOC_FAIL;
+ }
+ memset(data,0,size*sizeof(t_complex));
+
+ if(fourier->type&BWD) return F_NOT_SUPPORTED;
+
+ /* stupid way */
+ /*
+ for(v=0;v<Y;v++) {
+ off_f=v*X;
+ for(u=0;u<X;u++) {
+ for(y=0;y<Y;y++) {
+ off_r=y*X;
+ for(x=0;x<X;x++) {
+ arg=-2.0*M_PI*((1.0*u*x)/X+(1.0*v*y)/Y);
+ fourier->ftdata[off_f+u].r+=(cos(arg)*fourier->data[off_r+x].r-sin(arg)*fourier->data[off_r+x].i);
+ fourier->ftdata[off_f+u].i+=(sin(arg)*fourier->data[off_r+x].r+cos(arg)*fourier->data[off_r+x].i);
+ }
+ }
+ fourier->ftdata[off_f+u].r/=(X*Y);
+ fourier->ftdata[off_f+u].i/=(X*Y);
+ }
+ }
+ */
+
+ /* the more clever way ... */
+ // dft on index 1
+ off_f=0;
+ for(y=0;y<Y;y++) {
+ for(u=0;u<X;u++) {
+ //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,y);
+ for(x=0;x<X;x++) {
+ arg=-2.0*M_PI*u*x/X;
+ off_r=off_f+x;
+ data[off_f+u].r+=(cos(arg)*fourier->data[off_r].r-sin(arg)*fourier->data[off_r].i);
+ data[off_f+u].r+=(sin(arg)*fourier->data[off_r].r+cos(arg)*fourier->data[off_r].i);
+ }
+ data[off_f+u].r/=X;
+ data[off_f+u].i/=X;
+ }
+ off_f+=X;
+ }
+ // dft on index 2 of 'index 1 transformed data'
+ for(x=0;x<X;x++) {
+ off_f=x;
+ for(v=0;v<Y;v++) {
+ //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
+ off_r=x;
+ for(y=0;y<Y;y++) {
+ arg=-2.0*M_PI*v*y/Y;
+ fourier->ftdata[off_f].r+=(cos(arg)*data[off_r].r-sin(arg)*data[off_r].i);
+ fourier->ftdata[off_f].r+=(sin(arg)*data[off_r].r+cos(arg)*data[off_r].i);
+ off_r+=X;
+ }
+ fourier->ftdata[off_f].r/=Y;
+ fourier->ftdata[off_f].i/=Y;
+ off_f+=X;
+ }
+ }
+
+ free(data);
+
+ return F_SUCCESS;