From b4d5ded242c481f736ce588df536ca3010b14976 Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Tue, 14 Nov 2006 15:32:22 +0000 Subject: [PATCH] fwd kissfftndr sorta works, but the packing is wrong, or at least different from FFTW --- test/testkiss.py | 38 +++++++------ tools/Makefile | 2 +- tools/fftutil.c | 67 +++++++++++++++++++--- tools/kiss_fftnd.c | 3 +- tools/kiss_fftndr.c | 135 ++++++++++++++++++++++++++++++++++++++++++++ tools/kiss_fftndr.h | 47 +++++++++++++++ 6 files changed, 264 insertions(+), 28 deletions(-) create mode 100644 tools/kiss_fftndr.c create mode 100644 tools/kiss_fftndr.h diff --git a/test/testkiss.py b/test/testkiss.py index c6a2c9f..87edce2 100755 --- a/test/testkiss.py +++ b/test/testkiss.py @@ -40,6 +40,11 @@ else: def dopack(x,cpx=1): x = Numeric.reshape( x, ( Numeric.size(x),) ) + + print 'packed=[', + print ' '.join([str(y) for y in x]), + print ']' + if cpx: s = ''.join( [ struct.pack(fmt*2,c.real,c.imag) for c in x ] ) else: @@ -75,25 +80,22 @@ def flatten(x): def randmat( ndims ): dims=[] for i in range( ndims ): - curdim = int( random.uniform(2,4) ) + curdim = int( random.uniform(4,6) ) + if doreal and i==0: + curdim = int(curdim/2)*2 # force even for first dimension of real dims.append( curdim ) return make_random(dims ) def test_fft(ndims): - if ndims == 1: - nfft = int(random.uniform(50,100)) - if doreal: - nfft = int(nfft/2)*2 - - x = Numeric.array(make_random( [ nfft ] ) ) - else: - x=randmat( ndims ) + x=randmat( ndims ) print 'dimensions=%s' % str( Numeric.shape(x) ), if doreal: xver = FFT.real_fftnd(x) else: xver = FFT.fftnd(x) + print 'x=',x + print 'xver=',xver x2=dofft(x) err = xver - x2 @@ -130,6 +132,7 @@ def dofft(x): p = popen2.Popen3(cmd ) + p.tochild.write( dopack( x , iscomp ) ) p.tochild.close() @@ -140,7 +143,11 @@ def dofft(x): res = scale * res p.wait() - return Numeric.reshape(res,dims) + try: + return Numeric.reshape(res,dims) + finally: + print 'len(res)=%d' % len(res) + print 'dims=%s' % str(dims) def main(): opts,args = getopt.getopt(sys.argv[1:],'r') @@ -149,14 +156,11 @@ def main(): global doreal doreal = opts.has_key('-r') - if doreal: - print 'Note: Real optimization not yet done for odd length ffts and multi-D' - test_fft(1) - else: - print 'Testing multi-dimensional FFTs' - for dim in range(1,4): - test_fft( dim ) + print 'Testing multi-dimensional FFTs' + for dim in range(1,4): + test_fft( dim ) if __name__ == "__main__": + random.seed(2); main() diff --git a/tools/Makefile b/tools/Makefile index 86c4d6d..28e8cbd 100644 --- a/tools/Makefile +++ b/tools/Makefile @@ -43,7 +43,7 @@ $(FASTFILTREAL): ../kiss_fft.c kiss_fastfir.c kiss_fftr.c $(FASTFILT): ../kiss_fft.c kiss_fastfir.c $(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+ -DFAST_FILT_UTIL -$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fftnd.c kiss_fftr.c +$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fftnd.c kiss_fftr.c kiss_fftndr.c $(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+ $(PSDPNG): ../kiss_fft.c psdpng.c kiss_fftr.c diff --git a/tools/fftutil.c b/tools/fftutil.c index c80ff72..73bb670 100644 --- a/tools/fftutil.c +++ b/tools/fftutil.c @@ -19,8 +19,7 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND #include #include "kiss_fft.h" -#include "kiss_fftnd.h" -#include "kiss_fftr.h" +#include "kiss_fftndr.h" static void fft_file(FILE * fin,FILE * fout,int nfft,int isinverse) @@ -62,6 +61,51 @@ void fft_filend(FILE * fin,FILE * fout,int *dims,int ndims,int isinverse) free (buf); } + + +#define CHK fprintf(stderr,"LINE=%d\t",__LINE__) + +static +void fft_filend_real(FILE * fin,FILE * fout,int *dims,int ndims,int isinverse) +{ + int dimprod=1,i; + kiss_fftndr_cfg st; + void *ibuf; + void *obuf; + int insize,outsize; // size in bytes + + for (i=0;i 0) { + if (isinverse) { + kiss_fftndri(st, + (kiss_fft_cpx*)ibuf, + (kiss_fft_scalar*)obuf); + }else{ + kiss_fftndr(st, + (kiss_fft_scalar*)ibuf, + (kiss_fft_cpx*)obuf); + } + fwrite (obuf, sizeof(kiss_fft_scalar), outsize,fout); + } + free(st); + free(ibuf); + free(obuf); +} + static void fft_file_real(FILE * fin,FILE * fout,int nfft,int isinverse) { @@ -100,7 +144,7 @@ int get_dims(char * arg,int * dims) if (p0) *p0++ = '\0'; dims[ndims++] = atoi(arg); - /* fprintf(stderr,"dims[%d] = %d\n",ndims-1,dims[ndims-1]); */ + fprintf(stderr,"dims[%d] = %d\n",ndims-1,dims[ndims-1]); arg = p0; }while (p0); return ndims; @@ -147,12 +191,17 @@ int main(int argc,char ** argv) ++optind; } - if (ndims>1 && !isreal) - fft_filend(fin,fout,dims,ndims,isinverse); - else if (ndims==1 && !isreal) - fft_file(fin,fout,dims[0],isinverse); - else if (ndims==1 && isreal) - fft_file_real(fin,fout,dims[0],isinverse); + if (ndims==1) { + if (isreal) + fft_file_real(fin,fout,dims[0],isinverse); + else + fft_file(fin,fout,dims[0],isinverse); + }else{ + if (isreal) + fft_filend_real(fin,fout,dims,ndims,isinverse); + else + fft_filend(fin,fout,dims,ndims,isinverse); + } if (fout!=stdout) fclose(fout); if (fin!=stdin) fclose(fin); diff --git a/tools/kiss_fftnd.c b/tools/kiss_fftnd.c index bf9e6b1..00ec783 100644 --- a/tools/kiss_fftnd.c +++ b/tools/kiss_fftnd.c @@ -22,7 +22,7 @@ struct kiss_fftnd_state{ int ndims; int *dims; kiss_fft_cfg *states; /* cfg states for each dimension */ - kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire buffer */ + 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) @@ -177,6 +177,7 @@ void kiss_fftnd(kiss_fftnd_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) for ( k=0; k < st->ndims; ++k) { int curdim = st->dims[k]; int stride = st->dimprod / curdim; + fprintf(stderr,"curdim = %d, stride = %d\n",curdim,stride); for ( i=0 ; istates[k], bufin+i , bufout+i*curdim, stride ); diff --git a/tools/kiss_fftndr.c b/tools/kiss_fftndr.c new file mode 100644 index 0000000..7a638e5 --- /dev/null +++ b/tools/kiss_fftndr.c @@ -0,0 +1,135 @@ +/* +Copyright (c) 2003-2004, Mark Borgerding + +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * 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. + * 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. + +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. +*/ + +#include "kiss_fftndr.h" +#include "_kiss_fft_guts.h" + + +// copied from kiss_fftnd.c +struct kiss_fftnd_state{ + int dimprod; + int ndims; + int *dims; + kiss_fft_cfg *states; /* cfg states for each dimension */ + kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire input */ +}; + + + +struct kiss_fftndr_state +{ + int dimReal; + int dimOther; + kiss_fftr_cfg cfg_r; + kiss_fftnd_cfg cfg_nd; + void * tmpbuf; +}; + +static int prod(const int *dims, int ndims) +{ + int x=1; + while (ndims--) + x *= *dims++; + return x; +} + +#define CHK fprintf(stderr,"line=%d\t",__LINE__) +kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem) +{ + kiss_fftndr_cfg st = NULL; + size_t nr=0 , nd=0,ntmp=0; + int dimReal = dims[0]; + int dimOther = prod(dims+1,ndims-1); + + (void)kiss_fftr_alloc(dims[0],inverse_fft,NULL,&nr); + (void)kiss_fftnd_alloc(dims+1,ndims-1,inverse_fft,NULL,&nd); + ntmp = + dimReal * sizeof(kiss_fft_scalar) // time buffer for one pass + + (dimReal+2) * sizeof(kiss_fft_scalar) // freq buffer for one pass + + dimOther*(dimReal+2) * sizeof(kiss_fft_scalar); // large enough to hold entire input in case of in-place + + + size_t memneeded = sizeof( struct kiss_fftndr_state ) + nr + nd + ntmp; + + if (lenmem==NULL) { + st = (kiss_fftndr_cfg) malloc(memneeded); + }else{ + if (*lenmem >= memneeded) + st = (kiss_fftndr_cfg)mem; + *lenmem = memneeded; + } + if (st==NULL) + return NULL; + memset( st , 0 , memneeded); + + st->dimReal = dimReal; + st->dimOther = dimOther; + st->cfg_r = kiss_fftr_alloc( dims[0],inverse_fft,st+1,&nr); + st->cfg_nd = kiss_fftnd_alloc(dims+1,ndims-1,inverse_fft, ((char*) st->cfg_r)+nr,&nd); + + // tell the N-D complex FFT to stride the input + st->cfg_nd->dimprod *= (dims[0]/2+1); + + st->tmpbuf = (char*)st->cfg_nd + nd; + + return st; +} + +void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) +{ + int k1,k2; + int dimReal = st->dimReal; + int dimOther = st->dimOther; + + kiss_fft_scalar * t1 = (kiss_fft_scalar*)st->tmpbuf; // enough to hold dimReal scalars + kiss_fft_cpx * t2 = (kiss_fft_cpx*)(t1+dimReal); // enough to hold the entire input (only used if in==out) + + fprintf(stderr,"dimReal=%d\n",dimReal); + fprintf(stderr,"dimOther=%d\n",dimOther); + fprintf(stderr,"t1=%p\n",t1); + fprintf(stderr,"t2=%p\n",t2); + + + // take the real data, fft it and transpose + for (k1=0;k1cfg_r, t1, t2 + k1*(dimReal/2+1) ); + + fprintf(stderr,"t1=fft([ " ); + for (k2=0;k2cfg_nd, t2,freqdata); + + fprintf(stderr,"out=[ " ); + for (k1=0;k1dimReal * st->dimOther); +} diff --git a/tools/kiss_fftndr.h b/tools/kiss_fftndr.h new file mode 100644 index 0000000..38ec3ab --- /dev/null +++ b/tools/kiss_fftndr.h @@ -0,0 +1,47 @@ +#ifndef KISS_NDR_H +#define KISS_NDR_H + +#include "kiss_fft.h" +#include "kiss_fftr.h" +#include "kiss_fftnd.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct kiss_fftndr_state *kiss_fftndr_cfg; + + +kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem); +/* + dims[0] must be even + + If you don't care to allocate space, use mem = lenmem = NULL +*/ + + +void kiss_fftndr( + kiss_fftndr_cfg cfg, + const kiss_fft_scalar *timedata, + kiss_fft_cpx *freqdata); +/* + input timedata has dims[0] X dims[1] X ... X dims[ndims-1] scalar points + output freqdata has dims[0] X dims[1] X ... X dims[ndims-1]/2+1 complex points +*/ + +void kiss_fftndri( + kiss_fftndr_cfg cfg, + const kiss_fft_cpx *freqdata, + kiss_fft_scalar *timedata); +/* + input and output dimensions are the exact opposite of kiss_fftndr +*/ + + +#define kiss_fftr_free free + +#ifdef __cplusplus +} +#endif + +#endif