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 }