int fourier_alloc_data(t_fourier *fourier) {
- int i;
-
- 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;
- }
+ int i,size;
+
+ size=1;
+ for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
+
+ 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");
- for(i=0;i<fourier->dim;i++) {
- free(fourier->data[i]);
- free(fourier->ftdata[i]);
- }
+ free(fourier->data);
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=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];
+ 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=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);
- }
+ 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);
}
}
}
}
int fourier_dft_2d(t_fourier *fourier) {
- return 0;
+ return F_NOT_SUPPORTED;
}
int fourier_dft_3d(t_fourier *fourier) {
- return 0;
+ return F_NOT_SUPPORTED;
}
-
-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;
- }
-}