mirror/kissfft
1
0
Fork 0
mirror of https://github.com/mborgerding/kissfft.git synced 2025-05-27 21:20:27 -04:00
Code Issues Packages Projects Releases Wiki Activity

fwd kissfftndr sorta works, but the packing is wrong, or at least different from FFTW

Browse Source
This commit is contained in:
Mark Borgerding 2006-11-14 15:32:22 +00:00
parent 7f68a2916b
commit b4d5ded242
6 changed files with 264 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
test/testkiss.py
View File

@ -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()

2
tools/Makefile
View File

@ -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

67
tools/fftutil.c
View File

@ -19,8 +19,7 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
#include <unistd.h>
#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<ndims;++i)
dimprod *= dims[i];
insize = outsize = dimprod;
if (isinverse)
insize = insize*2*(dims[0]/2+1)/dims[0];
else
outsize = outsize*2*(dims[0]/2+1)/dims[0];
fprintf(stderr,"insize=%d outsize=%d\n",insize,outsize);
ibuf = malloc(insize*sizeof(kiss_fft_scalar));
obuf = malloc(outsize*sizeof(kiss_fft_scalar));
st = kiss_fftndr_alloc(dims, ndims, isinverse, 0, 0);
while ( fread (ibuf, sizeof(kiss_fft_scalar), insize, fin) > 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);

3
tools/kiss_fftnd.c
View File

@ -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 ; i<stride ; ++i )
kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride );

135
tools/kiss_fftndr.c Normal file
View File

@ -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;k1<dimOther;++k1) {
for (k2=0;k2<dimReal;++k2)
t1[k2] = timedata[k1 + k2*dimOther]; // select the appropriate samples into a contiguous buffer
//then fft the N real samples to create N/2+1 complex points
kiss_fftr( st->cfg_r, t1, t2 + k1*(dimReal/2+1) );
fprintf(stderr,"t1=fft([ " );
for (k2=0;k2<dimReal;++k2)
fprintf(stderr,"%f ",t1[k2] );
fprintf(stderr,"])\n" );
}
// fft the remaining dimensions just like the N-D complex case
kiss_fftnd(st->cfg_nd, t2,freqdata);
fprintf(stderr,"out=[ " );
for (k1=0;k1<dimOther;++k1) {
for (k2=0;k2<dimReal/2+1;++k2) {
kiss_fft_cpx c = freqdata[k1*(dimReal/2+1) +k2];
fprintf(stderr,"%f%+fi ",c.r , c.i );
}
fprintf(stderr,"\n" );
}
fprintf(stderr,"]\n" );
}
void kiss_fftndri(kiss_fftndr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
{
memset(timedata,0,sizeof(kiss_fft_scalar) * st->dimReal * st->dimOther);
}

47
tools/kiss_fftndr.h Normal file
View File

@ -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