int fourier_init(t_fourier *fourier,int outfd) {
- fprintf(outfd,"[fourier] initializing fourier api ...\n");
+ dprintf(outfd,"[fourier] initializing fourier api ...\n");
fourier->outfd=outfd;
fourier->type=0;
int fourier_alloc_data(t_fourier *fourier) {
- int i;
+ int i,size;
- for(i=0;i<fourier->dim;i++) {
- fourier->data[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
- fourier->ftdata[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
- if((fourier->data[i]==NULL)||(fourier->ftdata[i]==NULL)) {
- fprintf(fourier->outfd,"[fourier] malloc failed\n");
- return F_ALLOC_FAIL;
- }
+ size=1;
+ for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
+
+ dprintf(fourier->outfd,"[fourier] allocating %d bytes for data ...\n",
+ 2*size*sizeof(t_complex));
+
+ fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
+ fourier->ftdata=&(fourier->data[size]);
+
+ memset(fourier->data,0,2*size*sizeof(t_complex));
+
+ if(fourier->data==NULL) {
+ dprintf(fourier->outfd,"[fourier] malloc failed\n");
+ return F_ALLOC_FAIL;
}
return F_SUCCESS;
int fourier_shutdown(t_fourier *fourier) {
- int i;
+ dprintf(fourier->outfd,"[fourier] shutdown\n");
+
+ free(fourier->data);
+
+ return F_SUCCESS;
+}
+
+/* fft functions */
+
+int fourier_fft_1d_init(t_fourier *fourier) {
+
+ unsigned int i,j,m,r;
+ int mask[32];
+
+ /* calculate m = log 2 (N) + check N = 2^m */
+ m=0;
+ while(m<32) {
+ if(1<<m==fourier->data_len[0]) break;
+ else m++;
+ }
+ if(m==32) {
+ dprintf(fourier->outfd,
+ "[fourier] fft impossible (N > 2^32 or not a power of 2\n");
+ return F_FFT_IMPOSSIBLE;
+ }
+ dprintf(fourier->outfd,"[fourier] log2(%d) = %d\n",fourier->data_len[0],m);
+ fourier->log2len[0]=m;
+
+ /* pointer array which will point to reversed t_complex data */
+ fourier->revdata=(t_complex **)malloc(fourier->data_len[0]*sizeof(t_complex *));
+ if(fourier->revdata==NULL) {
+ dprintf(fourier->outfd,
+ "[fourier] malloc for reversed data pointers failed\n");
+ return F_ALLOC_FAIL;
+ }
- fprintf(fourier->outfd,"[fourier] shutdown\n");
+ dprintf(fourier->outfd,"[fourier] fft: allocated data pointers\n");
- for(i=0;i<fourier->dim;i++) {
- free(fourier->data[dim]);
- free(fourier->ftdata[dim]);
+ /* swap data (bit reversal) */
+ for(i=0;i<m;i++) mask[i]=1<<i;
+ for(i=0;i<fourier->data_len[0];i++) {
+ r=0;
+ for(j=0;j<m;j++) r+=(((i&mask[j])>>j)<<(m-j-1));
+ fourier->revdata[i]=fourier->ftdata+r;
}
return F_SUCCESS;
}
+int fourier_fft_1d_shutdown(t_fourier *fourier) {
+
+ dprintf(fourier->outfd,"[fourier] fft shutdown\n");
+
+ free(fourier->revdata);
+
+ return F_SUCCESS;
+}
+
+
+int fourier_fft_1d(t_fourier *fourier) {
+
+ int i;
+
+ /* copy src to destination, destination is modified in place */
+ memcpy(fourier->ftdata,fourier->data,fourier->data_len[0]*sizeof(t_complex));
+ /* from now on access bit reversed data by revdata[i]->... */
+
+ // for(i=0;i<fourier->log2len[0];i++) {
+
+
+ // }
+
+
+ for(i=0;i<fourier->data_len[0];i++)
+ dprintf(fourier->outfd,"%f %f\n",
+ fourier->ftdata[i].r,fourier->revdata[i]->r);
+
+ return F_NOT_SUPPORTED;
+
+ return F_SUCCESS;
+}
+
int fourier_dft_1d(t_fourier *fourier) {
- int i,k,dim;
+ int i,k;
double arg;
if(fourier->type&FWD) {
- for(dim=1;dim<=fourier->dim;dim++) {
- for(k=0;k<fourier->data_len[dim];k++) {
- fourier->ftdata[dim][k].r=0;
- fourier->ftdata[dim][k].i=0;
- for(i=0;i<fourier->data_len[dim];i++) {
- /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
- arg=-1.0*k*M_PI*i/fourier->data_len[dim];
- fourier->ftdata[dim][k].r+=(cos(arg)*fourier->data[dim][i]-sin(arg)*fourier->data[dim][i].i);
- fourier->ftdata[dim][k].i+=(sin(arg)*fourier->data[dim][i]+cos(arg)*fourier->data[dim][i].i);
- }
- fourier->ftdata[dim][k].r/=fourier->data_len[dim];
- fourier->ftdata[dim][k].i/=fourier->data_len[dim];
+ for(k=0;k<fourier->data_len[0];k++) {
+ fourier->ftdata[k].r=0;
+ fourier->ftdata[k].i=0;
+ for(i=0;i<fourier->data_len[0];i++) {
+ /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
+ arg=-2.0*k*M_PI*i/fourier->data_len[0];
+ fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
+ fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
}
+ fourier->ftdata[k].r/=fourier->data_len[0];
+ fourier->ftdata[k].i/=fourier->data_len[0];
}
}
else {
- for(dim=1;dim<=fourier->dim;dim++) {
- for(k=0;k<fourier->data_len[dim];k++) {
- fourier->data[dim][k].r=0;
- fourier->data[dim][k].i=0;
- for(i=0;i<fourier->data_len[dim];i++) {
- arg=1.0*k*M_PI*i/fourier->data_len[dim];
- fourier->data[dim][k].r+=(cos(arg)*fourier->ftdata[dim][i]-sin(arg)*fourier->ftdata[dim][i].i);
- fourier->data[dim][k].i+=(sin(arg)*fourier->ftdata[dim][i]+cos(arg)*fourier->ftdata[dim][i].i);
- }
+ for(k=0;k<fourier->data_len[0];k++) {
+ fourier->data[k].r=0;
+ fourier->data[k].i=0;
+ for(i=0;i<fourier->data_len[0];i++) {
+ arg=2.0*k*M_PI*i/fourier->data_len[0];
+ fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
+ fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
}
}
}
return F_SUCCESS;
}
-int fourier_calc(t_fourier *fourier) {
-
- fprintf(fourier->outfd,"[fourier] %d dimensional %c-%cft calculation ...\n",
- fourier->dim,
- ((fourier->type&FWD)?'f':'b'),
- ((fourier->type&DFT)?'d':'f'));
-
- if(fourier->type&DFT) {
- switch(fourier->dim) {
- case 1:
- fourier_dft_1d(fourier);
- break;
- case 2:
- fourier_dft_2d(fourier);
- break;
- case 3:
- fourier_dft_3d(fourier);
- break;
- default:
- fprintf(fourier->outfd,"[fourier] dimension failure\n");
- return F_DIM_FAILURE;
+int fourier_dft_2d(t_fourier *fourier) {
+
+ 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);
}
- return F_SUCCESS;
}
- else {
- fprintf(fourier->outfd,"[fourier] fft not supported by now\n");
- return F_NOT_SUPPORTED;
+ */
+
+ /* 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;
+}
+
+int fourier_dft_3d(t_fourier *fourier) {
+ return F_NOT_SUPPORTED;
+}