]> hackdaworld.org Git - my-code/api.git/blob - fourier/fourier.c
modifications i forgot to commit
[my-code/api.git] / fourier / fourier.c
1 /* fourier.c -- fourier transformation api
2  *
3  * author: hackbard@hackdaworld.dyndns.org
4  *
5  */
6
7 #include "fourier.h"
8
9 int fourier_init(t_fourier *fourier,int outfd) {
10
11  dprintf(outfd,"[fourier] initializing fourier api ...\n");
12
13  fourier->outfd=outfd;
14  fourier->type=0;
15  fourier->dim=0;
16
17  return F_SUCCESS;
18 }
19
20 int fourier_alloc_data(t_fourier *fourier) {
21
22   int i,size;
23
24   size=1;
25   for(i=0;i<fourier->dim;i++) size*=fourier->data_len[i];
26
27   dprintf(fourier->outfd,"[fourier] allocating %d bytes for data ...\n",
28           2*size*sizeof(t_complex));
29
30   fourier->data=(t_complex *)malloc(2*size*sizeof(t_complex));
31   fourier->ftdata=&(fourier->data[size]);
32
33   memset(fourier->data,0,2*size*sizeof(t_complex));
34
35   if(fourier->data==NULL) {
36     dprintf(fourier->outfd,"[fourier] malloc failed\n");
37     return F_ALLOC_FAIL;
38   }
39
40   return F_SUCCESS;
41 }
42
43 int fourier_shutdown(t_fourier *fourier) {
44
45   dprintf(fourier->outfd,"[fourier] shutdown\n");
46
47   free(fourier->data);
48
49   return F_SUCCESS;
50 }
51
52 /* fft functions */
53
54 int fourier_fft_1d_init(t_fourier *fourier) {
55
56   unsigned int i,j,m,r;
57   int mask[32];
58
59   /* calculate m = log 2 (N) + check N = 2^m */
60   m=0;
61   while(m<32) {
62     if(1<<m==fourier->data_len[0]) break;
63     else m++;
64   }
65   if(m==32) {
66     dprintf(fourier->outfd,
67             "[fourier] fft impossible (N > 2^32 or not a power of 2\n");
68     return F_FFT_IMPOSSIBLE;
69   }
70   dprintf(fourier->outfd,"[fourier] log2(%d) = %d\n",fourier->data_len[0],m);
71   fourier->log2len[0]=m;
72
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");
78     return F_ALLOC_FAIL;
79   }
80
81   dprintf(fourier->outfd,"[fourier] fft: allocated data pointers\n");
82
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++) {
86     r=0;
87     for(j=0;j<m;j++) r+=(((i&mask[j])>>j)<<(m-j-1));
88     fourier->revdata[i]=fourier->ftdata+r;
89   }
90
91   return F_SUCCESS;
92 }
93
94 int fourier_fft_1d_shutdown(t_fourier *fourier) {
95
96   dprintf(fourier->outfd,"[fourier] fft shutdown\n");
97
98   free(fourier->revdata);
99
100   return F_SUCCESS;
101 }
102
103
104 int fourier_fft_1d(t_fourier *fourier) {
105
106   int i;
107
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]->... */
111
112   // for(i=0;i<fourier->log2len[0];i++) {
113
114
115   // }
116
117
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);
121
122   return F_NOT_SUPPORTED;
123
124   return F_SUCCESS;
125 }
126
127 int fourier_dft_1d(t_fourier *fourier) {
128
129   int i,k;
130   double arg;
131
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);
141       }
142       fourier->ftdata[k].r/=fourier->data_len[0];
143       fourier->ftdata[k].i/=fourier->data_len[0];
144     }
145   }
146   else {
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);
154       }
155     }
156   }
157
158   return F_SUCCESS;
159 }
160
161 int fourier_dft_2d(t_fourier *fourier) {
162
163   int off_f,off_r;
164   int u,v,x,y;
165   int X,Y;
166   int size;
167   double arg;
168   t_complex *data;
169
170   X=fourier->data_len[0];
171   Y=fourier->data_len[1];
172
173   size=X*Y;
174
175   data=(t_complex *)malloc(size*sizeof(t_complex));
176   if(data==NULL) {
177     dprintf(fourier->outfd,
178             "[fourier] malloc of %d bytes of temp data failed\n",size);
179     return F_ALLOC_FAIL;
180   }
181   memset(data,0,size*sizeof(t_complex));
182
183   if(fourier->type&BWD) return F_NOT_SUPPORTED;
184
185   /* stupid way */
186   /*
187   for(v=0;v<Y;v++) {
188     off_f=v*X;
189     for(u=0;u<X;u++) {
190       for(y=0;y<Y;y++) {
191         off_r=y*X;
192         for(x=0;x<X;x++) {
193          arg=-2.0*M_PI*((1.0*u*x)/X+(1.0*v*y)/Y);
194           fourier->ftdata[off_f+u].r+=(cos(arg)*fourier->data[off_r+x].r-sin(arg)*fourier->data[off_r+x].i);
195           fourier->ftdata[off_f+u].i+=(sin(arg)*fourier->data[off_r+x].r+cos(arg)*fourier->data[off_r+x].i);
196         }
197       }
198       fourier->ftdata[off_f+u].r/=(X*Y);
199       fourier->ftdata[off_f+u].i/=(X*Y);
200     }
201   }
202   */
203
204   /* the more clever way ... */
205   // dft on index 1
206   off_f=0;
207   for(y=0;y<Y;y++) {
208     for(u=0;u<X;u++) {
209       //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",u,y);
210       for(x=0;x<X;x++) {
211         arg=-2.0*M_PI*u*x/X;
212         off_r=off_f+x;
213         data[off_f+u].r+=(cos(arg)*fourier->data[off_r].r-sin(arg)*fourier->data[off_r].i);
214         data[off_f+u].r+=(sin(arg)*fourier->data[off_r].r+cos(arg)*fourier->data[off_r].i);
215       }
216       data[off_f+u].r/=X;
217       data[off_f+u].i/=X;
218     }
219     off_f+=X;
220   }
221   //  dft on index 2 of 'index 1 transformed data'
222   for(x=0;x<X;x++) {
223     off_f=x;
224     for(v=0;v<Y;v++) {
225       //dprintf(fourier->outfd,"[fourier] (u=%d,v=%d)\n",v,x);
226       off_r=x;
227       for(y=0;y<Y;y++) {
228         arg=-2.0*M_PI*v*y/Y;
229         fourier->ftdata[off_f].r+=(cos(arg)*data[off_r].r-sin(arg)*data[off_r].i);
230         fourier->ftdata[off_f].r+=(sin(arg)*data[off_r].r+cos(arg)*data[off_r].i);
231         off_r+=X;
232       }
233       fourier->ftdata[off_f].r/=Y;
234       fourier->ftdata[off_f].i/=Y;
235       off_f+=X;
236     }
237   }
238  
239   free(data);
240
241   return F_SUCCESS;
242 }
243
244 int fourier_dft_3d(t_fourier *fourier) {
245   return F_NOT_SUPPORTED;
246 }