File indexing completed on 2024-06-16 04:35:42
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_fftndr.h" 0016 #include "../_kiss_fft_guts.h" 0017 #define MAX(x,y) ( ( (x)<(y) )?(y):(x) ) 0018 0019 struct kiss_fftndr_state 0020 { 0021 int dimReal; 0022 int dimOther; 0023 kiss_fftr_cfg cfg_r; 0024 kiss_fftnd_cfg cfg_nd; 0025 void * tmpbuf; 0026 }; 0027 0028 static int prod(const int *dims, int ndims) 0029 { 0030 int x=1; 0031 while (ndims--) 0032 x *= *dims++; 0033 return x; 0034 } 0035 0036 kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem) 0037 { 0038 kiss_fftndr_cfg st = NULL; 0039 size_t nr=0 , nd=0,ntmp=0; 0040 int dimReal = dims[ndims-1]; 0041 int dimOther = prod(dims,ndims-1); 0042 size_t memneeded; 0043 0044 (void)kiss_fftr_alloc(dimReal,inverse_fft,NULL,&nr); 0045 (void)kiss_fftnd_alloc(dims,ndims-1,inverse_fft,NULL,&nd); 0046 ntmp = 0047 MAX( 2*dimOther , dimReal+2) * sizeof(kiss_fft_scalar) // freq buffer for one pass 0048 + dimOther*(dimReal+2) * sizeof(kiss_fft_scalar); // large enough to hold entire input in case of in-place 0049 0050 memneeded = sizeof( struct kiss_fftndr_state ) + nr + nd + ntmp; 0051 0052 if (lenmem==NULL) { 0053 st = (kiss_fftndr_cfg) malloc(memneeded); 0054 }else{ 0055 if (*lenmem >= memneeded) 0056 st = (kiss_fftndr_cfg)mem; 0057 *lenmem = memneeded; 0058 } 0059 if (st==NULL) 0060 return NULL; 0061 memset( st , 0 , memneeded); 0062 0063 st->dimReal = dimReal; 0064 st->dimOther = dimOther; 0065 st->cfg_r = kiss_fftr_alloc( dimReal,inverse_fft,st+1,&nr); 0066 st->cfg_nd = kiss_fftnd_alloc(dims,ndims-1,inverse_fft, ((char*) st->cfg_r)+nr,&nd); 0067 st->tmpbuf = (char*)st->cfg_nd + nd; 0068 0069 return st; 0070 } 0071 0072 void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) 0073 { 0074 int k1,k2; 0075 int dimReal = st->dimReal; 0076 int dimOther = st->dimOther; 0077 int nrbins = dimReal/2+1; 0078 0079 kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf; 0080 kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther); 0081 0082 // timedata is N0 x N1 x ... x Nk real 0083 0084 // take a real chunk of data, fft it and place the output at correct intervals 0085 for (k1=0;k1<dimOther;++k1) { 0086 kiss_fftr( st->cfg_r, timedata + k1*dimReal , tmp1 ); // tmp1 now holds nrbins complex points 0087 for (k2=0;k2<nrbins;++k2) 0088 tmp2[ k2*dimOther+k1 ] = tmp1[k2]; 0089 } 0090 0091 for (k2=0;k2<nrbins;++k2) { 0092 kiss_fftnd(st->cfg_nd, tmp2+k2*dimOther, tmp1); // tmp1 now holds dimOther complex points 0093 for (k1=0;k1<dimOther;++k1) 0094 freqdata[ k1*(nrbins) + k2] = tmp1[k1]; 0095 } 0096 } 0097 0098 void kiss_fftndri(kiss_fftndr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata) 0099 { 0100 int k1,k2; 0101 int dimReal = st->dimReal; 0102 int dimOther = st->dimOther; 0103 int nrbins = dimReal/2+1; 0104 kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf; 0105 kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther); 0106 0107 for (k2=0;k2<nrbins;++k2) { 0108 for (k1=0;k1<dimOther;++k1) 0109 tmp1[k1] = freqdata[ k1*(nrbins) + k2 ]; 0110 kiss_fftnd(st->cfg_nd, tmp1, tmp2+k2*dimOther); 0111 } 0112 0113 for (k1=0;k1<dimOther;++k1) { 0114 for (k2=0;k2<nrbins;++k2) 0115 tmp1[k2] = tmp2[ k2*dimOther+k1 ]; 0116 kiss_fftri( st->cfg_r,tmp1,timedata + k1*dimReal); 0117 } 0118 }