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 dprintf(fourier->outfd,"[fourier] allocating %d bytes for data ...\n",
28 2*size*sizeof(t_complex));
30 fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
31 fourier->ftdata=&(fourier->data[size]);
33 memset(fourier->data,0,2*size*sizeof(t_complex));
35 if(fourier->data==NULL) {
36 dprintf(fourier->outfd,"[fourier] malloc failed\n");
43 int fourier_shutdown(t_fourier *fourier) {
45 dprintf(fourier->outfd,"[fourier] shutdown\n");
54 int fourier_fft_1d_init(t_fourier *fourier) {
59 /* calculate m = log 2 (N) + check N = 2^m */
62 if(1<<m==fourier->data_len[0]) break;
66 dprintf(fourier->outfd,
67 "[fourier] fft impossible (N > 2^32 or not a power of 2\n");
68 return F_FFT_IMPOSSIBLE;
70 dprintf(fourier->outfd,"[fourier] log2(%d) = %d\n",fourier->data_len[0],m);
71 fourier->log2len[0]=m;
73 /* pointer array which will point to reversed t_complex data */
74 fourier->revdata=(t_complex **)malloc(fourier->data_len[0]*sizeof(t_complex *));
75 if(fourier->revdata==NULL) {
76 dprintf(fourier->outfd,
77 "[fourier] malloc for reversed data pointers failed\n");
81 dprintf(fourier->outfd,"[fourier] fft: allocated data pointers\n");
83 /* swap data (bit reversal) */
84 for(i=0;i<m;i++) mask[i]=1<<i;
85 for(i=0;i<fourier->data_len[0];i++) {
87 for(j=0;j<m;j++) r+=(((i&mask[j])>>j)<<(m-j-1));
88 fourier->revdata[i]=fourier->ftdata+r;
94 int fourier_fft_1d_shutdown(t_fourier *fourier) {
96 dprintf(fourier->outfd,"[fourier] fft shutdown\n");
98 free(fourier->revdata);
104 int fourier_fft_1d(t_fourier *fourier) {
108 /* copy src to destination, destination is modified in place */
109 memcpy(fourier->ftdata,fourier->data,fourier->data_len[0]*sizeof(t_complex));
110 /* from now on access bit reversed data by revdata[i]->... */
112 // for(i=0;i<fourier->log2len[0];i++) {
118 for(i=0;i<fourier->data_len[0];i++)
119 dprintf(fourier->outfd,"%f %f\n",
120 fourier->ftdata[i].r,fourier->revdata[i]->r);
122 return F_NOT_SUPPORTED;
127 int fourier_dft_1d(t_fourier *fourier) {
132 if(fourier->type&FWD) {
133 for(k=0;k<fourier->data_len[0];k++) {
134 fourier->ftdata[k].r=0;
135 fourier->ftdata[k].i=0;
136 for(i=0;i<fourier->data_len[0];i++) {
137 /* f(k) = 1/N sum(n=0-N) f(n) exp(-i*k*2*PI*n/N) */
138 arg=-2.0*k*M_PI*i/fourier->data_len[0];
139 fourier->ftdata[k].r+=(cos(arg)*fourier->data[i].r-sin(arg)*fourier->data[i].i);
140 fourier->ftdata[k].i+=(sin(arg)*fourier->data[i].r+cos(arg)*fourier->data[i].i);
142 fourier->ftdata[k].r/=fourier->data_len[0];
143 fourier->ftdata[k].i/=fourier->data_len[0];
147 for(k=0;k<fourier->data_len[0];k++) {
148 fourier->data[k].r=0;
149 fourier->data[k].i=0;
150 for(i=0;i<fourier->data_len[0];i++) {
151 arg=2.0*k*M_PI*i/fourier->data_len[0];
152 fourier->data[k].r+=(cos(arg)*fourier->ftdata[i].r-sin(arg)*fourier->ftdata[i].i);
153 fourier->data[k].i+=(sin(arg)*fourier->ftdata[i].r+cos(arg)*fourier->ftdata[i].i);
161 int fourier_dft_2d(t_fourier *fourier) {
170 X=fourier->data_len[0];
171 Y=fourier->data_len[1];
175 data=(t_complex *)malloc(size*sizeof(t_complex));
177 dprintf(fourier->outfd,
178 "[fourier] malloc of %d bytes of temp data failed\n",size);
181 memset(data,0,size*sizeof(t_complex));
183 if(fourier->type&BWD) return F_NOT_SUPPORTED;
187 // for(v=0;v<Y;v++) {
189 // for(u=0;u<X;u++) {
190 // dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,v);
191 // fourier->ftdata[off_f+u].r=0;
192 // fourier->ftdata[off_f+u].i=0;
193 // for(y=0;y<Y;y++) {
195 // for(x=0;x<X;x++) {
196 // arg=-2.0*M_PI*((1.0*u*x)/X+(1.0*v*y)/Y);
197 // fourier->ftdata[off_f+u].r+=(cos(arg)*fourier->data[off_r+x].r-sin(arg)*fourier->data[off_r+x].i);
198 // fourier->ftdata[off_f+u].i+=(sin(arg)*fourier->data[off_r+x].r+cos(arg)*fourier->data[off_r+x].i);
201 // fourier->ftdata[off_f+u].r/=(X*Y);
202 // fourier->ftdata[off_f+u].i/=(X*Y);
206 /* the more clever way ... */
211 //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,y);
215 data[off_f+u].r+=(cos(arg)*fourier->data[off_r].r-sin(arg)*fourier->data[off_r].i);
216 data[off_f+u].r+=(sin(arg)*fourier->data[off_r].r+cos(arg)*fourier->data[off_r].i);
223 // dft on index 2 of 'index 1 transformed data'
227 //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
231 fourier->ftdata[off_f+v].r+=(cos(arg)*data[off_r].r-sin(arg)*data[off_r].i);
232 fourier->ftdata[off_f+v].r+=(sin(arg)*data[off_r].r+cos(arg)*data[off_r].i);
234 fourier->ftdata[off_f+v].r/=Y;
235 fourier->ftdata[off_f+v].i/=Y;
245 int fourier_dft_3d(t_fourier *fourier) {
246 return F_NOT_SUPPORTED;