1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
/*
* Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#include "kiss_fftnd.h"
#include "_kiss_fft_guts.h"
struct kiss_fftnd_state{
int dimprod; /* dimsum would be mighty tasty right now */
int ndims;
int *dims;
kiss_fft_cfg *states; /* cfg states for each dimension */
kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire input */
};
kiss_fftnd_cfg kiss_fftnd_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem)
{
kiss_fftnd_cfg st = NULL;
int i;
int dimprod=1;
size_t memneeded = sizeof(struct kiss_fftnd_state);
char * ptr;
for (i=0;i<ndims;++i) {
size_t sublen=0;
kiss_fft_alloc (dims[i], inverse_fft, NULL, &sublen);
memneeded += sublen; /* st->states[i] */
dimprod *= dims[i];
}
memneeded += sizeof(int) * ndims;/* st->dims */
memneeded += sizeof(void*) * ndims;/* st->states */
memneeded += sizeof(kiss_fft_cpx) * dimprod; /* st->tmpbuf */
if (lenmem == NULL) {/* allocate for the caller*/
st = (kiss_fftnd_cfg) malloc (memneeded);
} else { /* initialize supplied buffer if big enough */
if (*lenmem >= memneeded)
st = (kiss_fftnd_cfg) mem;
*lenmem = memneeded; /*tell caller how big struct is (or would be) */
}
if (!st)
return NULL; /*malloc failed or buffer too small */
st->dimprod = dimprod;
st->ndims = ndims;
ptr=(char*)(st+1);
st->states = (kiss_fft_cfg *)ptr;
ptr += sizeof(void*) * ndims;
st->dims = (int*)ptr;
ptr += sizeof(int) * ndims;
st->tmpbuf = (kiss_fft_cpx*)ptr;
ptr += sizeof(kiss_fft_cpx) * dimprod;
for (i=0;i<ndims;++i) {
size_t len;
st->dims[i] = dims[i];
kiss_fft_alloc (st->dims[i], inverse_fft, NULL, &len);
st->states[i] = kiss_fft_alloc (st->dims[i], inverse_fft, ptr,&len);
ptr += len;
}
/*
Hi there!
If you're looking at this particular code, it probably means you've got a brain-dead bounds checker
that thinks the above code overwrites the end of the array.
It doesn't.
-- Mark
P.S.
The below code might give you some warm fuzzies and help convince you.
*/
if ( ptr - (char*)st != (int)memneeded ) {
fprintf(stderr,
"################################################################################\n"
"Internal error! Memory allocation miscalculation\n"
"################################################################################\n"
);
}
return st;
}
/*
This works by tackling one dimension at a time.
In effect,
Each stage starts out by reshaping the matrix into a DixSi 2d matrix.
A Di-sized fft is taken of each column, transposing the matrix as it goes.
Here's a 3-d example:
Take a 2x3x4 matrix, laid out in memory as a contiguous buffer
[ [ [ a b c d ] [ e f g h ] [ i j k l ] ]
[ [ m n o p ] [ q r s t ] [ u v w x ] ] ]
Stage 0 ( D=2): treat the buffer as a 2x12 matrix
[ [a b ... k l]
[m n ... w x] ]
FFT each column with size 2.
Transpose the matrix at the same time using kiss_fft_stride.
[ [ a+m a-m ]
[ b+n b-n]
...
[ k+w k-w ]
[ l+x l-x ] ]
Note fft([x y]) == [x+y x-y]
Stage 1 ( D=3) treats the buffer (the output of stage D=2) as an 3x8 matrix,
[ [ a+m a-m b+n b-n c+o c-o d+p d-p ]
[ e+q e-q f+r f-r g+s g-s h+t h-t ]
[ i+u i-u j+v j-v k+w k-w l+x l-x ] ]
And perform FFTs (size=3) on each of the columns as above, transposing
the matrix as it goes. The output of stage 1 is
(Legend: ap = [ a+m e+q i+u ]
am = [ a-m e-q i-u ] )
[ [ sum(ap) fft(ap)[0] fft(ap)[1] ]
[ sum(am) fft(am)[0] fft(am)[1] ]
[ sum(bp) fft(bp)[0] fft(bp)[1] ]
[ sum(bm) fft(bm)[0] fft(bm)[1] ]
[ sum(cp) fft(cp)[0] fft(cp)[1] ]
[ sum(cm) fft(cm)[0] fft(cm)[1] ]
[ sum(dp) fft(dp)[0] fft(dp)[1] ]
[ sum(dm) fft(dm)[0] fft(dm)[1] ] ]
Stage 2 ( D=4) treats this buffer as a 4*6 matrix,
[ [ sum(ap) fft(ap)[0] fft(ap)[1] sum(am) fft(am)[0] fft(am)[1] ]
[ sum(bp) fft(bp)[0] fft(bp)[1] sum(bm) fft(bm)[0] fft(bm)[1] ]
[ sum(cp) fft(cp)[0] fft(cp)[1] sum(cm) fft(cm)[0] fft(cm)[1] ]
[ sum(dp) fft(dp)[0] fft(dp)[1] sum(dm) fft(dm)[0] fft(dm)[1] ] ]
Then FFTs each column, transposing as it goes.
The resulting matrix is the 3d FFT of the 2x3x4 input matrix.
Note as a sanity check that the first element of the final
stage's output (DC term) is
sum( [ sum(ap) sum(bp) sum(cp) sum(dp) ] )
, i.e. the summation of all 24 input elements.
*/
void kiss_fftnd(kiss_fftnd_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
{
int i,k;
const kiss_fft_cpx * bufin=fin;
kiss_fft_cpx * bufout;
/*arrange it so the last bufout == fout*/
if ( st->ndims & 1 ) {
bufout = fout;
if (fin==fout) {
memcpy( st->tmpbuf, fin, sizeof(kiss_fft_cpx) * st->dimprod );
bufin = st->tmpbuf;
}
}else
bufout = st->tmpbuf;
for ( k=0; k < st->ndims; ++k) {
int curdim = st->dims[k];
int stride = st->dimprod / curdim;
for ( i=0 ; i<stride ; ++i )
kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride );
/*toggle back and forth between the two buffers*/
if (bufout == st->tmpbuf){
bufout = fout;
bufin = st->tmpbuf;
}else{
bufout = st->tmpbuf;
bufin = fout;
}
}
}
|
