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 | |
13 | struct 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 | |
22 | static int prod(const int *dims, int ndims) |
23 | { |
24 | int x=1; |
25 | while (ndims--) |
26 | x *= *dims++; |
27 | return x; |
28 | } |
29 | |
30 | kiss_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 | |
74 | void 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 | |
100 | void 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 | |