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++) {
- // dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,v);
- // fourier->ftdata[off_f+u].r=0;
- // fourier->ftdata[off_f+u].i=0;
- // 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);
- // }
- // }
+ /*
+ 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+=X;
}
// dft on index 2 of 'index 1 transformed data'
- off_f=0;
for(x=0;x<X;x++) {
+ off_f=0;
for(v=0;v<Y;v++) {
//dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
+ off_r=0;
for(y=0;y<Y;y++) {
arg=-2.0*M_PI*v*y/Y;
- off_r=off_f+y;
- fourier->ftdata[off_f+v].r+=(cos(arg)*data[off_r].r-sin(arg)*data[off_r].i);
- fourier->ftdata[off_f+v].r+=(sin(arg)*data[off_r].r+cos(arg)*data[off_r].i);
+ fourier->ftdata[off_f+x].r+=(cos(arg)*data[off_r+x].r-sin(arg)*data[off_r+x].i);
+ fourier->ftdata[off_f+x].r+=(sin(arg)*data[off_r+x].r+cos(arg)*data[off_r+x].i);
+ off_r+=X;
}
- fourier->ftdata[off_f+v].r/=Y;
- fourier->ftdata[off_f+v].i/=Y;
+ fourier->ftdata[off_f+x].r/=Y;
+ fourier->ftdata[off_f+x].i/=Y;
+ off_f+=X;
}
- off_f+=Y;
}
free(data);