1 /* fourier.c -- fourier transformation api
3 * author: hackbard@hackdaworld.dyndns.org
9 int fourier_init(t_fourier *fourier,int outfd) {
11 fprintf(outfd,"[fourier] initializing fourier api ...\n");
20 int fourier_alloc_data(t_fourier *fourier) {
24 for(i=0;i<fourier->dim;i++) {
25 fourier->data[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
26 fourier->ftdata[i]=(t_complex *)malloc(fourier->data_len[i]*sizeof(t_complex));
27 if((fourier->data[i]==NULL)||(fourier->ftdata[i]==NULL)) {
28 fprintf(fourier->outfd,"[fourier] malloc failed\n");
36 int fourier_shutdown(t_fourier *fourier) {
40 fprintf(fourier->outfd,"[fourier] shutdown\n");
42 for(i=0;i<fourier->dim;i++) {
43 free(fourier->data[dim]);
44 free(fourier->ftdata[dim]);
50 int fourier_dft_1d(t_fourier *fourier) {
55 if(fourier->type&FWD) {
56 for(dim=1;dim<=fourier->dim;dim++) {
57 for(k=0;k<fourier->data_len[dim];k++) {
58 fourier->ftdata[dim][k].r=0;
59 fourier->ftdata[dim][k].i=0;
60 for(i=0;i<fourier->data_len[dim];i++) {
61 /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
62 arg=-1.0*k*M_PI*i/fourier->data_len[dim];
63 fourier->ftdata[dim][k].r+=(cos(arg)*fourier->data[dim][i]-sin(arg)*fourier->data[dim][i].i);
64 fourier->ftdata[dim][k].i+=(sin(arg)*fourier->data[dim][i]+cos(arg)*fourier->data[dim][i].i);
66 fourier->ftdata[dim][k].r/=fourier->data_len[dim];
67 fourier->ftdata[dim][k].i/=fourier->data_len[dim];
72 for(dim=1;dim<=fourier->dim;dim++) {
73 for(k=0;k<fourier->data_len[dim];k++) {
74 fourier->data[dim][k].r=0;
75 fourier->data[dim][k].i=0;
76 for(i=0;i<fourier->data_len[dim];i++) {
77 arg=1.0*k*M_PI*i/fourier->data_len[dim];
78 fourier->data[dim][k].r+=(cos(arg)*fourier->ftdata[dim][i]-sin(arg)*fourier->ftdata[dim][i].i);
79 fourier->data[dim][k].i+=(sin(arg)*fourier->ftdata[dim][i]+cos(arg)*fourier->ftdata[dim][i].i);
88 int fourier_calc(t_fourier *fourier) {
90 fprintf(fourier->outfd,"[fourier] %d dimensional %c-%cft calculation ...\n",
92 ((fourier->type&FWD)?'f':'b'),
93 ((fourier->type&DFT)?'d':'f'));
95 if(fourier->type&DFT) {
96 switch(fourier->dim) {
98 fourier_dft_1d(fourier);
101 fourier_dft_2d(fourier);
104 fourier_dft_3d(fourier);
107 fprintf(fourier->outfd,"[fourier] dimension failure\n");
108 return F_DIM_FAILURE;
113 fprintf(fourier->outfd,"[fourier] fft not supported by now\n");
114 return F_NOT_SUPPORTED;