1 /* fourier.c -- fourier transformation api
3 * author: hackbard@hackdaworld.dyndns.org
9 int fourier_init(t_fourier *fourier,int outfd) {
11 dprintf(outfd,"[fourier] initializing fourier api ...\n");
20 int fourier_alloc_data(t_fourier *fourier) {
25 for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
27 fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
28 fourier->ftdata=&(fourier->data[size]);
30 memset(fourier->data,0,2*size*sizeof(t_complex));
32 if(fourier->data==NULL) {
33 dprintf(fourier->outfd,"[fourier] malloc failed\n");
40 int fourier_shutdown(t_fourier *fourier) {
42 dprintf(fourier->outfd,"[fourier] shutdown\n");
49 int fourier_fft_1d(t_fourier *fourier) {
58 int fourier_dft_1d(t_fourier *fourier) {
63 if(fourier->type&FWD) {
64 for(k=0;k<fourier->data_len[0];k++) {
65 fourier->ftdata[k].r=0;
66 fourier->ftdata[k].i=0;
67 for(i=0;i<fourier->data_len[0];i++) {
68 /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
69 arg=-2.0*k*M_PI*i/fourier->data_len[0];
70 fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
71 fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
73 fourier->ftdata[k].r/=fourier->data_len[0];
74 fourier->ftdata[k].i/=fourier->data_len[0];
78 for(k=0;k<fourier->data_len[0];k++) {
81 for(i=0;i<fourier->data_len[0];i++) {
82 arg=2.0*k*M_PI*i/fourier->data_len[0];
83 fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
84 fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
92 int fourier_dft_2d(t_fourier *fourier) {
93 return F_NOT_SUPPORTED;
96 int fourier_dft_3d(t_fourier *fourier) {
97 return F_NOT_SUPPORTED;