File indexing completed on 2024-06-09 04:49:45
0001 /* 0002 SPDX-FileCopyrightText: 2003-2004 Mark Borgerding <Mark@Borgerding.net> 0003 SPDX-License-Identifier: BSD-3-Clause 0004 All rights reserved. 0005 0006 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 0007 0008 * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 0009 * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 0010 * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. 0011 0012 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 0013 */ 0014 0015 #include "kiss_fftnd.h" 0016 #include "../_kiss_fft_guts.h" 0017 0018 struct kiss_fftnd_state{ 0019 int dimprod; /* dimsum would be mighty tasty right now */ 0020 int ndims; 0021 int *dims; 0022 kiss_fft_cfg *states; /* cfg states for each dimension */ 0023 kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire input */ 0024 }; 0025 0026 kiss_fftnd_cfg kiss_fftnd_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem) 0027 { 0028 kiss_fftnd_cfg st = NULL; 0029 int i; 0030 int dimprod=1; 0031 size_t memneeded = sizeof(struct kiss_fftnd_state); 0032 char * ptr; 0033 0034 for (i=0;i<ndims;++i) { 0035 size_t sublen=0; 0036 kiss_fft_alloc (dims[i], inverse_fft, NULL, &sublen); 0037 memneeded += sublen; /* st->states[i] */ 0038 dimprod *= dims[i]; 0039 } 0040 memneeded += sizeof(int) * ndims;/* st->dims */ 0041 memneeded += sizeof(void*) * ndims;/* st->states */ 0042 memneeded += sizeof(kiss_fft_cpx) * dimprod; /* st->tmpbuf */ 0043 0044 if (lenmem == NULL) {/* allocate for the caller*/ 0045 st = (kiss_fftnd_cfg) malloc (memneeded); 0046 } else { /* initialize supplied buffer if big enough */ 0047 if (*lenmem >= memneeded) 0048 st = (kiss_fftnd_cfg) mem; 0049 *lenmem = memneeded; /*tell caller how big struct is (or would be) */ 0050 } 0051 if (!st) 0052 return NULL; /*malloc failed or buffer too small */ 0053 0054 st->dimprod = dimprod; 0055 st->ndims = ndims; 0056 ptr=(char*)(st+1); 0057 0058 st->states = (kiss_fft_cfg *)ptr; 0059 ptr += sizeof(void*) * ndims; 0060 0061 st->dims = (int*)ptr; 0062 ptr += sizeof(int) * ndims; 0063 0064 st->tmpbuf = (kiss_fft_cpx*)ptr; 0065 ptr += sizeof(kiss_fft_cpx) * dimprod; 0066 0067 for (i=0;i<ndims;++i) { 0068 size_t len; 0069 st->dims[i] = dims[i]; 0070 kiss_fft_alloc (st->dims[i], inverse_fft, NULL, &len); 0071 st->states[i] = kiss_fft_alloc (st->dims[i], inverse_fft, ptr,&len); 0072 ptr += len; 0073 } 0074 /* 0075 Hi there! 0076 0077 If you're looking at this particular code, it probably means you've got a brain-dead bounds checker 0078 that thinks the above code overwrites the end of the array. 0079 0080 It doesn't. 0081 0082 -- Mark 0083 0084 P.S. 0085 The below code might give you some warm fuzzies and help convince you. 0086 */ 0087 if ( ptr - (char*)st != (int)memneeded ) { 0088 fprintf(stderr, 0089 "################################################################################\n" 0090 "Internal error! Memory allocation miscalculation\n" 0091 "################################################################################\n" 0092 ); 0093 } 0094 return st; 0095 } 0096 0097 /* 0098 This works by tackling one dimension at a time. 0099 0100 In effect, 0101 Each stage starts out by reshaping the matrix into a DixSi 2d matrix. 0102 A Di-sized fft is taken of each column, transposing the matrix as it goes. 0103 0104 Here's a 3-d example: 0105 Take a 2x3x4 matrix, laid out in memory as a contiguous buffer 0106 [ [ [ a b c d ] [ e f g h ] [ i j k l ] ] 0107 [ [ m n o p ] [ q r s t ] [ u v w x ] ] ] 0108 0109 Stage 0 ( D=2): treat the buffer as a 2x12 matrix 0110 [ [a b ... k l] 0111 [m n ... w x] ] 0112 0113 FFT each column with size 2. 0114 Transpose the matrix at the same time using kiss_fft_stride. 0115 0116 [ [ a+m a-m ] 0117 [ b+n b-n] 0118 ... 0119 [ k+w k-w ] 0120 [ l+x l-x ] ] 0121 0122 Note fft([x y]) == [x+y x-y] 0123 0124 Stage 1 ( D=3) treats the buffer (the output of stage D=2) as an 3x8 matrix, 0125 [ [ a+m a-m b+n b-n c+o c-o d+p d-p ] 0126 [ e+q e-q f+r f-r g+s g-s h+t h-t ] 0127 [ i+u i-u j+v j-v k+w k-w l+x l-x ] ] 0128 0129 And perform FFTs (size=3) on each of the columns as above, transposing 0130 the matrix as it goes. The output of stage 1 is 0131 (Legend: ap = [ a+m e+q i+u ] 0132 am = [ a-m e-q i-u ] ) 0133 0134 [ [ sum(ap) fft(ap)[0] fft(ap)[1] ] 0135 [ sum(am) fft(am)[0] fft(am)[1] ] 0136 [ sum(bp) fft(bp)[0] fft(bp)[1] ] 0137 [ sum(bm) fft(bm)[0] fft(bm)[1] ] 0138 [ sum(cp) fft(cp)[0] fft(cp)[1] ] 0139 [ sum(cm) fft(cm)[0] fft(cm)[1] ] 0140 [ sum(dp) fft(dp)[0] fft(dp)[1] ] 0141 [ sum(dm) fft(dm)[0] fft(dm)[1] ] ] 0142 0143 Stage 2 ( D=4) treats this buffer as a 4*6 matrix, 0144 [ [ sum(ap) fft(ap)[0] fft(ap)[1] sum(am) fft(am)[0] fft(am)[1] ] 0145 [ sum(bp) fft(bp)[0] fft(bp)[1] sum(bm) fft(bm)[0] fft(bm)[1] ] 0146 [ sum(cp) fft(cp)[0] fft(cp)[1] sum(cm) fft(cm)[0] fft(cm)[1] ] 0147 [ sum(dp) fft(dp)[0] fft(dp)[1] sum(dm) fft(dm)[0] fft(dm)[1] ] ] 0148 0149 Then FFTs each column, transposing as it goes. 0150 0151 The resulting matrix is the 3d FFT of the 2x3x4 input matrix. 0152 0153 Note as a sanity check that the first element of the final 0154 stage's output (DC term) is 0155 sum( [ sum(ap) sum(bp) sum(cp) sum(dp) ] ) 0156 , i.e. the summation of all 24 input elements. 0157 0158 */ 0159 void kiss_fftnd(kiss_fftnd_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) 0160 { 0161 int i,k; 0162 const kiss_fft_cpx * bufin=fin; 0163 kiss_fft_cpx * bufout; 0164 0165 /*arrange it so the last bufout == fout*/ 0166 if ( st->ndims & 1 ) { 0167 bufout = fout; 0168 if (fin==fout) { 0169 memcpy( st->tmpbuf, fin, sizeof(kiss_fft_cpx) * st->dimprod ); 0170 bufin = st->tmpbuf; 0171 } 0172 }else 0173 bufout = st->tmpbuf; 0174 0175 for ( k=0; k < st->ndims; ++k) { 0176 int curdim = st->dims[k]; 0177 int stride = st->dimprod / curdim; 0178 0179 for ( i=0 ; i<stride ; ++i ) 0180 kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride ); 0181 0182 /*toggle back and forth between the two buffers*/ 0183 if (bufout == st->tmpbuf){ 0184 bufout = fout; 0185 bufin = st->tmpbuf; 0186 }else{ 0187 bufout = st->tmpbuf; 0188 bufin = fout; 0189 } 0190 } 0191 }