+ */
+
+ /* 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);