1 /**/
2 #include "ltfat.h"
3 #include "ltfat/types.h"
4 #include "ltfat/macros.h"
5 #include "ltfat/thirdparty/kiss_fft.h"
6
7 /****** FFT ******/
LTFAT_NAME(fft_plan)8 struct LTFAT_NAME(fft_plan)
9 {
10 ltfat_int L;
11 ltfat_int W;
12 LTFAT_COMPLEX* in;
13 LTFAT_COMPLEX* out;
14 LTFAT_COMPLEX* tmp;
15 LTFAT_KISS(fft_plan)* kiss_plan;
16 };
17
18 LTFAT_API int
LTFAT_NAME(fft)19 LTFAT_NAME(fft)(LTFAT_COMPLEX in[], ltfat_int L, ltfat_int W,
20 LTFAT_COMPLEX out[])
21 {
22 LTFAT_NAME(fft_plan)* p = NULL;
23 int status = LTFATERR_SUCCESS;
24
25 CHECKSTATUS( LTFAT_NAME(fft_init)(L, W, in, out, 0, &p));
26 LTFAT_NAME(fft_execute)(p);
27 LTFAT_NAME(fft_done)(&p);
28 error:
29 return status;
30 }
31
32 static int
LTFAT_NAME(fft_init_common)33 LTFAT_NAME(fft_init_common)(ltfat_int L, ltfat_int W,
34 LTFAT_COMPLEX in[], LTFAT_COMPLEX out[],
35 unsigned inverse, LTFAT_NAME(fft_plan)** p)
36 {
37 LTFAT_NAME(fft_plan)* fftwp = NULL;
38 ltfat_int nextfastL = 0;
39
40 int status = LTFATERR_SUCCESS;
41
42 CHECKNULL(p);
43 CHECK(LTFATERR_NOTPOSARG, L > 0, "L must be positive");
44 CHECK(LTFATERR_NOTPOSARG, W > 0, "W must be positive");
45
46 nextfastL = ltfat_nextfastfft(L);
47
48 if (L != nextfastL)
49 {
50 DEBUG("Warning: L=%td is a \"slow\" FFT lengh. "
51 "Next fast FFT lenght is L=%td. See ltfat_nextfastfft. "
52 "Moreover, some dynamic memory allocation will occur "
53 "during execution.", L, nextfastL);
54 }
55
56 CHECKMEM( fftwp = LTFAT_NEW(LTFAT_NAME(fft_plan)) );
57 fftwp->L = L; fftwp->W = W; fftwp->in = in; fftwp->out = out;
58
59 fftwp->kiss_plan = LTFAT_KISS(fft_alloc)(L, inverse, NULL, NULL);
60 CHECKINIT(fftwp->kiss_plan, "FFTW plan creation failed.");
61
62 if (in == out)
63 CHECKMEM( fftwp->tmp = LTFAT_NAME_COMPLEX(malloc)(L) );
64
65 *p = fftwp;
66 return status;
67 error:
68 if (fftwp)
69 LTFAT_NAME(fft_done)(&fftwp);
70 *p = NULL;
71 return status;
72 }
73
74
75 LTFAT_API int
LTFAT_NAME(fft_init)76 LTFAT_NAME(fft_init)(ltfat_int L, ltfat_int W,
77 LTFAT_COMPLEX in[], LTFAT_COMPLEX out[],
78 unsigned UNUSED(flags), LTFAT_NAME(fft_plan)** p)
79 {
80 return LTFAT_NAME(fft_init_common)(L, W, in, out, 0, p);
81 }
82
83 LTFAT_API int
LTFAT_NAME(fft_execute)84 LTFAT_NAME(fft_execute)(LTFAT_NAME(fft_plan)* p)
85 {
86 return LTFAT_NAME(fft_execute_newarray)( p, p->in, p->out);
87 }
88
89 LTFAT_API int
LTFAT_NAME(fft_execute_newarray)90 LTFAT_NAME(fft_execute_newarray)(LTFAT_NAME(fft_plan)* p,
91 const LTFAT_COMPLEX in[], LTFAT_COMPLEX out[])
92 {
93 int status = LTFATERR_SUCCESS;
94 CHECKNULL(p); CHECKNULL(in); CHECKNULL(out);
95
96 if (in == out)
97 {
98 CHECKNULL(p->tmp);
99
100 for (ltfat_int w = 0; w < p->W; w++)
101 {
102 memcpy(p->tmp, in + w * p->L, p->L * sizeof * p->tmp);
103 LTFAT_KISS(fft)(p->kiss_plan,
104 (const kiss_fft_cpx*) p->tmp,
105 (kiss_fft_cpx*) out + w * p->L);
106 }
107 }
108 else
109 {
110 for (ltfat_int w = 0; w < p->W; w++)
111 LTFAT_KISS(fft)(p->kiss_plan,
112 (const kiss_fft_cpx*) in + w * p->L,
113 (kiss_fft_cpx*) out + w * p->L);
114 }
115
116 error:
117 return status;
118 }
119
120 LTFAT_API int
LTFAT_NAME(fft_done)121 LTFAT_NAME(fft_done)(LTFAT_NAME(fft_plan)** p)
122 {
123 int status = LTFATERR_SUCCESS;
124 LTFAT_NAME(fft_plan)* pp = NULL;
125 CHECKNULL(p); CHECKNULL(*p);
126 pp = *p;
127 if (pp->tmp) ltfat_free(pp->tmp);
128 if (pp->kiss_plan) ltfat_free(pp->kiss_plan);
129 ltfat_free(pp);
130 pp = NULL;
131 error:
132 return status;
133 }
134
135 /******* IFFT ******/
LTFAT_NAME(ifft_plan)136 struct LTFAT_NAME(ifft_plan)
137 {
138 struct LTFAT_NAME(fft_plan) inplan;
139 };
140
141 LTFAT_API int
LTFAT_NAME(ifft)142 LTFAT_NAME(ifft)(LTFAT_COMPLEX in[], ltfat_int L, ltfat_int W,
143 LTFAT_COMPLEX out[])
144 {
145 LTFAT_NAME(ifft_plan)* p = NULL;
146 int status = LTFATERR_SUCCESS;
147
148 CHECKSTATUS( LTFAT_NAME(ifft_init)(L, W, in, out, 0, &p));
149 LTFAT_NAME(ifft_execute)(p);
150 LTFAT_NAME(ifft_done)(&p);
151 error:
152 return status;
153 }
154
155 LTFAT_API int
LTFAT_NAME(ifft_init)156 LTFAT_NAME(ifft_init)(ltfat_int L, ltfat_int W,
157 LTFAT_COMPLEX in[], LTFAT_COMPLEX out[],
158 unsigned UNUSED(flags), LTFAT_NAME(ifft_plan)** p)
159 {
160 return LTFAT_NAME(fft_init_common)(L, W, in, out, 1,
161 (LTFAT_NAME(fft_plan)**) p);
162 }
163
164 LTFAT_API int
LTFAT_NAME(ifft_execute)165 LTFAT_NAME(ifft_execute)(LTFAT_NAME(ifft_plan)* p)
166 {
167 return LTFAT_NAME(fft_execute)((LTFAT_NAME(fft_plan)*) p);
168 }
169
170 LTFAT_API int
LTFAT_NAME(ifft_execute_newarray)171 LTFAT_NAME(ifft_execute_newarray)(LTFAT_NAME(ifft_plan)* p,
172 const LTFAT_COMPLEX in[], LTFAT_COMPLEX out[])
173 {
174 return LTFAT_NAME(fft_execute_newarray)((LTFAT_NAME(fft_plan)*) p, in, out);
175
176 }
177
178 LTFAT_API int
LTFAT_NAME(ifft_done)179 LTFAT_NAME(ifft_done)(LTFAT_NAME(ifft_plan)** p)
180 {
181 return LTFAT_NAME(fft_done)((LTFAT_NAME(fft_plan)**) p);
182 }
183
184 /****** FFTREAL ******/
LTFAT_NAME(fftreal_plan)185 struct LTFAT_NAME(fftreal_plan)
186 {
187 ltfat_int L;
188 ltfat_int W;
189 LTFAT_REAL* in;
190 LTFAT_REAL* out;
191 LTFAT_COMPLEX* tmp;
192 LTFAT_KISS(fft_plan)* kiss_plan_cpx;
193 LTFAT_KISS(fftr_plan)* kiss_plan;
194 };
195
196 LTFAT_API int
LTFAT_NAME(fftreal)197 LTFAT_NAME(fftreal)(LTFAT_REAL in[], ltfat_int L, ltfat_int W,
198 LTFAT_COMPLEX out[])
199 {
200 LTFAT_NAME(fftreal_plan)* p = NULL;
201 int status = LTFATERR_SUCCESS;
202
203 CHECKSTATUS( LTFAT_NAME(fftreal_init)(L, W, in, out, 0, &p));
204 LTFAT_NAME(fftreal_execute)(p);
205 LTFAT_NAME(fftreal_done)(&p);
206 error:
207 return status;
208 }
209
210 static int
LTFAT_NAME(fftreal_init_common)211 LTFAT_NAME(fftreal_init_common)(ltfat_int L, ltfat_int W,
212 LTFAT_REAL in[], LTFAT_REAL out[],
213 unsigned inverse, LTFAT_NAME(fftreal_plan)** p)
214 {
215 LTFAT_NAME(fftreal_plan)* fftwp = NULL;
216 ltfat_int M2;
217 ltfat_int nextfastL;
218
219 int status = LTFATERR_SUCCESS;
220
221 CHECKNULL(p);
222 CHECK(LTFATERR_BADARG, L > 0, "L must be even positive");
223 CHECK(LTFATERR_NOTPOSARG, W > 0, "W must be positive");
224
225 M2 = L / 2 + 1;
226
227 CHECKMEM( fftwp = LTFAT_NEW(LTFAT_NAME(fftreal_plan)) );
228 fftwp->L = L; fftwp->W = W; fftwp->in = in; fftwp->out = out;
229
230 nextfastL = ltfat_nextfastfft(L);
231
232 if (L != nextfastL)
233 {
234 DEBUG("Warning: L=%td is a \"slow\" FFT lengh. "
235 "Next fast FFT lenght is L=%td. See ltfat_nextfastfft. "
236 "Moreover, some dynamic memory allocation will occur "
237 "during execution.", L, nextfastL);
238 }
239
240 if (L % 2)
241 {
242 if (L != nextfastL)
243 {
244 DEBUGNOTE("Warning: Odd L is a \"very slow\" FFT lengh. Full FFT will be performed.");
245 }
246 // Workaround for odd-length transforms
247 fftwp->kiss_plan_cpx = LTFAT_KISS(fft_alloc)(L, inverse, NULL, NULL);
248 CHECKINIT(fftwp->kiss_plan_cpx, "FFTW plan creation failed.");
249 CHECKMEM( fftwp->tmp = LTFAT_NAME_COMPLEX(malloc)(4 * M2 ) );
250 }
251 else
252 {
253 fftwp->kiss_plan = LTFAT_KISS(fftr_alloc)(L, inverse, NULL, NULL);
254 CHECKINIT(fftwp->kiss_plan, "FFTW plan creation failed.");
255
256 if (in == out)
257 CHECKMEM( fftwp->tmp = LTFAT_NAME_COMPLEX(malloc)(L / 2 + 1) );
258 }
259 *p = fftwp;
260 return status;
261 error:
262 if (fftwp)
263 LTFAT_NAME(fftreal_done)(&fftwp);
264 *p = NULL;
265 return status;
266 }
267
268
269 LTFAT_API int
LTFAT_NAME(fftreal_init)270 LTFAT_NAME(fftreal_init)(ltfat_int L, ltfat_int W,
271 LTFAT_REAL in[], LTFAT_COMPLEX out[],
272 unsigned UNUSED(flags), LTFAT_NAME(fftreal_plan)** p)
273 {
274 return LTFAT_NAME(fftreal_init_common)(L, W, in, (LTFAT_REAL*) out, 0, p);
275 }
276
277 LTFAT_API int
LTFAT_NAME(fftreal_execute)278 LTFAT_NAME(fftreal_execute)(LTFAT_NAME(fftreal_plan)* p)
279 {
280 return LTFAT_NAME(fftreal_execute_newarray)( p, p->in,
281 (LTFAT_COMPLEX*) p->out);
282 }
283
284 LTFAT_API int
LTFAT_NAME(fftreal_execute_newarray)285 LTFAT_NAME(fftreal_execute_newarray)(LTFAT_NAME(fftreal_plan)* p,
286 const LTFAT_REAL in[], LTFAT_COMPLEX out[])
287 {
288 int status = LTFATERR_SUCCESS;
289 ltfat_int M2;
290
291 CHECKNULL(p); CHECKNULL(in); CHECKNULL(out);
292
293 M2 = p->L / 2 + 1;
294
295 if (p->L % 2)
296 {
297 ltfat_int step = p->L;
298 if (in == (const LTFAT_REAL*) out)
299 step = 2 * M2;
300
301 for (ltfat_int w = 0; w < p->W; w++)
302 {
303 LTFAT_NAME(real2complex_array)(in + w * step, p->L, p->tmp);
304
305 LTFAT_KISS(fft)(p->kiss_plan_cpx,
306 (const kiss_fft_cpx*) p->tmp,
307 (kiss_fft_cpx*) p->tmp + 2 * M2);
308
309 memcpy(out + w * M2, p->tmp + 2 * M2, M2 * sizeof * out);
310 }
311 }
312 else
313 {
314 if ( in == (const LTFAT_REAL*) out )
315 {
316 CHECKNULL(p->tmp);
317
318 for (ltfat_int w = 0; w < p->W; w++)
319 {
320 memcpy(p->tmp, in + w * 2 * M2, p->L * sizeof * p->in);
321 LTFAT_KISS(fftr)(p->kiss_plan,
322 (const kiss_fft_scalar*) p->tmp,
323 (kiss_fft_cpx*) out + w * M2);
324 }
325 }
326 else
327 {
328 for (ltfat_int w = 0; w < p->W; w++)
329 LTFAT_KISS(fftr)(p->kiss_plan,
330 (const kiss_fft_scalar*) in + w * p->L,
331 (kiss_fft_cpx*) out + w * M2);
332 }
333 }
334 error:
335 return status;
336 }
337
338 LTFAT_API int
LTFAT_NAME(fftreal_done)339 LTFAT_NAME(fftreal_done)(LTFAT_NAME(fftreal_plan)** p)
340 {
341 int status = LTFATERR_SUCCESS;
342 LTFAT_NAME(fftreal_plan)* pp = NULL;
343 CHECKNULL(p); CHECKNULL(*p);
344 pp = *p;
345 if (pp->tmp) ltfat_free(pp->tmp);
346 if (pp->kiss_plan) ltfat_free(pp->kiss_plan);
347 if (pp->kiss_plan_cpx) ltfat_free(pp->kiss_plan_cpx);
348 ltfat_free(pp);
349 pp = NULL;
350 error:
351 return status;
352 }
353
354 /******* IFFTREAL ******/
LTFAT_NAME(ifftreal_plan)355 struct LTFAT_NAME(ifftreal_plan)
356 {
357 struct LTFAT_NAME(fftreal_plan) inplan;
358 };
359
360 LTFAT_API int
LTFAT_NAME(ifftreal)361 LTFAT_NAME(ifftreal)(LTFAT_COMPLEX in[], ltfat_int L, ltfat_int W,
362 LTFAT_REAL out[])
363 {
364 LTFAT_NAME(ifftreal_plan)* p = NULL;
365 int status = LTFATERR_SUCCESS;
366
367 CHECKSTATUS( LTFAT_NAME(ifftreal_init)(L, W, in, out, 0, &p));
368 LTFAT_NAME(ifftreal_execute)(p);
369 LTFAT_NAME(ifftreal_done)(&p);
370 error:
371 return status;
372 }
373
374 LTFAT_API int
LTFAT_NAME(ifftreal_init)375 LTFAT_NAME(ifftreal_init)(ltfat_int L, ltfat_int W,
376 LTFAT_COMPLEX in[], LTFAT_REAL out[],
377 unsigned UNUSED(flags), LTFAT_NAME(ifftreal_plan)** p)
378 {
379 return LTFAT_NAME(fftreal_init_common)(L, W, (LTFAT_REAL*)in, out, 1,
380 (LTFAT_NAME(fftreal_plan)**) p);
381 }
382
383 LTFAT_API int
LTFAT_NAME(ifftreal_execute)384 LTFAT_NAME(ifftreal_execute)(LTFAT_NAME(ifftreal_plan)* pin)
385 {
386 LTFAT_NAME(fftreal_plan)* p = (LTFAT_NAME(fftreal_plan)*) pin;
387 return LTFAT_NAME(ifftreal_execute_newarray)( pin, (const LTFAT_COMPLEX*) p->in,
388 p->out);
389 }
390
391 LTFAT_API int
LTFAT_NAME(ifftreal_execute_newarray)392 LTFAT_NAME(ifftreal_execute_newarray)(LTFAT_NAME(ifftreal_plan)* pin,
393 const LTFAT_COMPLEX in[], LTFAT_REAL out[])
394 {
395 int status = LTFATERR_SUCCESS;
396 ltfat_int M2;
397 LTFAT_NAME(fftreal_plan)* p;
398 CHECKNULL(pin); CHECKNULL(in); CHECKNULL(out);
399 p = (LTFAT_NAME(fftreal_plan)*) pin;
400
401 M2 = p->L / 2 + 1;
402
403 if (p->L % 2)
404 {
405 ltfat_int step = p->L;
406 if (in == (const LTFAT_COMPLEX*) out)
407 step = 2 * M2;
408
409 for (ltfat_int w = 0; w < p->W; w++)
410 {
411 const LTFAT_COMPLEX* inTmp = in + w * M2;
412 memcpy(p->tmp, inTmp, M2 * sizeof * in);
413
414 for (ltfat_int ii = p->L - 1, jj = 1; ii >= M2; ii--, jj++)
415 p->tmp[ii] = conj(inTmp[jj]);
416
417 LTFAT_KISS(fft)(p->kiss_plan_cpx,
418 (const kiss_fft_cpx*) p->tmp,
419 (kiss_fft_cpx*) p->tmp + 2 * M2);
420
421 LTFAT_NAME(complex2real_array)( p->tmp + 2 * M2, p->L, out + w * step);
422 }
423 }
424 else
425 {
426 if (in == (const LTFAT_COMPLEX*) out)
427 {
428 CHECKNULL(p->tmp);
429
430 for (ltfat_int w = 0; w < p->W; w++)
431 {
432 memcpy(p->tmp, in + w * M2, M2 * sizeof * in);
433 LTFAT_KISS(fftri)(p->kiss_plan,
434 (const kiss_fft_cpx*) p->tmp,
435 (kiss_fft_scalar*) out + w * 2 * M2);
436 }
437 }
438 else
439 {
440 for (ltfat_int w = 0; w < p->W; w++)
441 LTFAT_KISS(fftri)(p->kiss_plan,
442 (const kiss_fft_cpx*) in + w * M2,
443 (kiss_fft_scalar*) out + w * p->L);
444 }
445 }
446 error:
447 return status;
448 }
449
450 LTFAT_API int
LTFAT_NAME(ifftreal_done)451 LTFAT_NAME(ifftreal_done)(LTFAT_NAME(ifftreal_plan)** p)
452 {
453 return LTFAT_NAME(fftreal_done)((LTFAT_NAME(fftreal_plan)**) p);
454 }
455