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;
return F_SUCCESS;
}
-int fourier_set_datalen(t_fourier *fourier,int dim,int len) {
+int fourier_alloc_data(t_fourier *fourier) {
int i;
- if(dim>MAX_DIM) {
- fprintf(fourier->outfd,"[fourier] dimension %d not supported\n",dim);
- return F_SET_DLEN_ERROR;
+ 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)) {
+ dprintf(fourier->outfd,"[fourier] malloc failed\n");
+ return F_ALLOC_FAIL;
+ }
}
- fourier->data_len[dim]=len;
-
return F_SUCCESS;
}
-int fourier_malloc(t_fourier *fourier) {
+int fourier_shutdown(t_fourier *fourier) {
int i;
- if(fourier->dim==0) {
- fprintf(fourier->outfd,"[fourier] dimension not specified\n");
- return F_DIM_ERROR;
- }
+ dprintf(fourier->outfd,"[fourier] shutdown\n");
+
for(i=0;i<fourier->dim;i++) {
- if(fourier->data_len[i]==0) {
- fprintf(fourier->outfd,"[fourier] invalid data len for dim %d\n",i);
- return F_DLEN_ERROR;
- }
+ free(fourier->data[i]);
+ free(fourier->ftdata[i]);
}
- for(i=0;i<fourier->dim;i++) {
- if((fourier->data[i]=(t_complex *)malloc(sizeof(t_complex)*fourier->data_len[i]))==NULL) {
- fprintf(fourier->outfd,"[fourier] malloc for data dim %d failed\n",i);
- return F_MALLOC_ERROR;
+ return F_SUCCESS;
+}
+
+int fourier_dft_1d(t_fourier *fourier) {
+
+ int i,k,dim;
+ double arg;
+
+ if(fourier->type&FWD) {
+ for(dim=0;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].r-sin(arg)*fourier->data[dim][i].i);
+ fourier->ftdata[dim][k].i+=(sin(arg)*fourier->data[dim][i].r+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];
+ }
+ }
+ }
+ else {
+ for(dim=0;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].r-sin(arg)*fourier->ftdata[dim][i].i);
+ fourier->data[dim][k].i+=(sin(arg)*fourier->ftdata[dim][i].r+cos(arg)*fourier->ftdata[dim][i].i);
+ }
+ }
}
}
return F_SUCCESS;
}
+
+int fourier_dft_2d(t_fourier *fourier) {
+ return 0;
+}
+
+int fourier_dft_3d(t_fourier *fourier) {
+ return 0;
+}
+
+int fourier_calc(t_fourier *fourier) {
+
+ dprintf(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:
+ dprintf(fourier->outfd,"[fourier] dimension failure\n");
+ return F_DIM_FAILURE;
+ }
+ return F_SUCCESS;
+ }
+ else {
+ dprintf(fourier->outfd,"[fourier] fft not supported by now\n");
+ return F_NOT_SUPPORTED;
+ }
+}