1/*
2 * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
3 * This file is part of KISS FFT - https://github.com/mborgerding/kissfft
4 *
5 * SPDX-License-Identifier: BSD-3-Clause
6 * See COPYING file for more information.
7 */
8
9#include "kiss_fftndr.h"
10#include "_kiss_fft_guts.h"
11#define MAX(x,y) ( ( (x)<(y) )?(y):(x) )
12
13struct kiss_fftndr_state
14{
15 int dimReal;
16 int dimOther;
17 kiss_fftr_cfg cfg_r;
18 kiss_fftnd_cfg cfg_nd;
19 void * tmpbuf;
20};
21
22static int prod(const int *dims, int ndims)
23{
24 int x=1;
25 while (ndims--)
26 x *= *dims++;
27 return x;
28}
29
30kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem)
31{
32 KISS_FFT_ALIGN_CHECK(mem)
33
34 kiss_fftndr_cfg st = NULL;
35 size_t nr=0 , nd=0,ntmp=0;
36 int dimReal = dims[ndims-1];
37 int dimOther = prod(dims,ndims-1);
38 size_t memneeded;
39 char * ptr = NULL;
40
41 (void)kiss_fftr_alloc(dimReal,inverse_fft,NULL,&nr);
42 (void)kiss_fftnd_alloc(dims,ndims-1,inverse_fft,NULL,&nd);
43 ntmp =
44 MAX( 2*dimOther , dimReal+2) * sizeof(kiss_fft_scalar) // freq buffer for one pass
45 + dimOther*(dimReal+2) * sizeof(kiss_fft_scalar); // large enough to hold entire input in case of in-place
46
47 memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof( struct kiss_fftndr_state )) + KISS_FFT_ALIGN_SIZE_UP(nr) + KISS_FFT_ALIGN_SIZE_UP(nd) + KISS_FFT_ALIGN_SIZE_UP(ntmp);
48
49 if (lenmem==NULL) {
50 ptr = (char*) malloc(memneeded);
51 }else{
52 if (*lenmem >= memneeded)
53 ptr = (char *)mem;
54 *lenmem = memneeded;
55 }
56 if (ptr==NULL)
57 return NULL;
58
59 st = (kiss_fftndr_cfg) ptr;
60 memset( st , 0 , memneeded);
61 ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftndr_state));
62
63 st->dimReal = dimReal;
64 st->dimOther = dimOther;
65 st->cfg_r = kiss_fftr_alloc( dimReal,inverse_fft,ptr,&nr);
66 ptr += KISS_FFT_ALIGN_SIZE_UP(nr);
67 st->cfg_nd = kiss_fftnd_alloc(dims,ndims-1,inverse_fft, ptr,&nd);
68 ptr += KISS_FFT_ALIGN_SIZE_UP(nd);
69 st->tmpbuf = ptr;
70
71 return st;
72}
73
74void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
75{
76 int k1,k2;
77 int dimReal = st->dimReal;
78 int dimOther = st->dimOther;
79 int nrbins = dimReal/2+1;
80
81 kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf;
82 kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther);
83
84 // timedata is N0 x N1 x ... x Nk real
85
86 // take a real chunk of data, fft it and place the output at correct intervals
87 for (k1=0;k1<dimOther;++k1) {
88 kiss_fftr( st->cfg_r, timedata + k1*dimReal , tmp1 ); // tmp1 now holds nrbins complex points
89 for (k2=0;k2<nrbins;++k2)
90 tmp2[ k2*dimOther+k1 ] = tmp1[k2];
91 }
92
93 for (k2=0;k2<nrbins;++k2) {
94 kiss_fftnd(st->cfg_nd, tmp2+k2*dimOther, tmp1); // tmp1 now holds dimOther complex points
95 for (k1=0;k1<dimOther;++k1)
96 freqdata[ k1*(nrbins) + k2] = tmp1[k1];
97 }
98}
99
100void kiss_fftndri(kiss_fftndr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
101{
102 int k1,k2;
103 int dimReal = st->dimReal;
104 int dimOther = st->dimOther;
105 int nrbins = dimReal/2+1;
106 kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf;
107 kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther);
108
109 for (k2=0;k2<nrbins;++k2) {
110 for (k1=0;k1<dimOther;++k1)
111 tmp1[k1] = freqdata[ k1*(nrbins) + k2 ];
112 kiss_fftnd(st->cfg_nd, tmp1, tmp2+k2*dimOther);
113 }
114
115 for (k1=0;k1<dimOther;++k1) {
116 for (k2=0;k2<nrbins;++k2)
117 tmp1[k2] = tmp2[ k2*dimOther+k1 ];
118 kiss_fftri( st->cfg_r,tmp1,timedata + k1*dimReal);
119 }
120}
121