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, inverse N-d real FFTs now work to the best of my knowledge

Browse Source
This commit is contained in:
Mark Borgerding 2006-11-14 18:57:44 +00:00
parent b4d5ded242
commit 97ce553a94
4 changed files with 68 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
test/testkiss.py
View File

@ -41,9 +41,9 @@ else:
def dopack(x,cpx=1):
x = Numeric.reshape( x, ( Numeric.size(x),) )
print 'packed=[',
print ' '.join([str(y) for y in x]),
print ']'
#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 ] )
@ -80,9 +80,9 @@ def flatten(x):
def randmat( ndims ):
dims=[]
for i in range( ndims ):
curdim = int( random.uniform(4,6) )
if doreal and i==0:
curdim = int(curdim/2)*2 # force even for first dimension of real
curdim = int( random.uniform(4,7) )
if doreal and i==(ndims-1):
curdim = int(curdim/2)*2 # force even last dimension if real
dims.append( curdim )
return make_random(dims )
@ -94,9 +94,11 @@ def test_fft(ndims):
xver = FFT.real_fftnd(x)
else:
xver = FFT.fftnd(x)
print 'x=',x
print 'xver=',xver
#print 'x=',x
#print 'xver=',xver
open('/tmp/fftexp.dat','w').write(dopack( flatten(xver) , True ) )
x2=dofft(x)
err = xver - x2
errf = flatten(err)
@ -132,22 +134,20 @@ def dofft(x):
p = popen2.Popen3(cmd )
open('/tmp/fftin.dat','w').write(dopack( x , iscomp ) )
p.tochild.write( dopack( x , iscomp ) )
p.tochild.close()
res = dounpack( p.fromchild.read() , 1 )
open('/tmp/fftout.dat','w').write(dopack( flatten(res) , True ) )
if doreal:
dims[-1] = int( dims[-1]/2 ) + 1
res = scale * res
p.wait()
try:
return Numeric.reshape(res,dims)
finally:
print 'len(res)=%d' % len(res)
print 'dims=%s' % str(dims)
return Numeric.reshape(res,dims)
def main():
opts,args = getopt.getopt(sys.argv[1:],'r')
@ -156,11 +156,15 @@ def main():
global doreal
doreal = opts.has_key('-r')
print 'Testing multi-dimensional FFTs'
for dim in range(1,4):
if doreal:
print 'Testing multi-dimensional real FFTs'
else:
print 'Testing multi-dimensional FFTs'
for dim in range(1,5):
test_fft( dim )
if __name__ == "__main__":
random.seed(2);
#random.seed(2);
main()

8
tools/fftutil.c
View File

@ -77,12 +77,12 @@ void fft_filend_real(FILE * fin,FILE * fout,int *dims,int ndims,int isinverse)
for (i=0;i<ndims;++i)
dimprod *= dims[i];
insize = outsize = dimprod;
int rdim = dims[ndims-1];
if (isinverse)
insize = insize*2*(dims[0]/2+1)/dims[0];
insize = insize*2*(rdim/2+1)/rdim;
else
outsize = outsize*2*(dims[0]/2+1)/dims[0];
fprintf(stderr,"insize=%d outsize=%d\n",insize,outsize);
outsize = outsize*2*(rdim/2+1)/rdim;
ibuf = malloc(insize*sizeof(kiss_fft_scalar));
obuf = malloc(outsize*sizeof(kiss_fft_scalar));
@ -144,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;

1
tools/kiss_fftnd.c
View File

@ -177,7 +177,6 @@ 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 );

100
tools/kiss_fftndr.c
View File

@ -14,18 +14,7 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
#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 */
};
#define MAX(x,y) ( ( (x)<(y) )?(y):(x) )
struct kiss_fftndr_state
{
@ -44,22 +33,19 @@ static int prod(const int *dims, int ndims)
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);
int dimReal = dims[ndims-1];
int dimOther = prod(dims,ndims-1);
(void)kiss_fftr_alloc(dims[0],inverse_fft,NULL,&nr);
(void)kiss_fftnd_alloc(dims+1,ndims-1,inverse_fft,NULL,&nd);
(void)kiss_fftr_alloc(dimReal,inverse_fft,NULL,&nr);
(void)kiss_fftnd_alloc(dims,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
MAX( 2*dimOther , 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) {
@ -75,12 +61,8 @@ kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void
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->cfg_r = kiss_fftr_alloc( dimReal,inverse_fft,st+1,&nr);
st->cfg_nd = kiss_fftnd_alloc(dims,ndims-1,inverse_fft, ((char*) st->cfg_r)+nr,&nd);
st->tmpbuf = (char*)st->cfg_nd + nd;
return st;
@ -92,44 +74,50 @@ void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx
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)
kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf;
kiss_fft_cpx * tmp2;
if (dimReal/2+1 > dimOther)
tmp2 = tmp1+dimReal/2+1;
else
tmp2 = tmp1+dimOther;
fprintf(stderr,"dimReal=%d\n",dimReal);
fprintf(stderr,"dimOther=%d\n",dimOther);
fprintf(stderr,"t1=%p\n",t1);
fprintf(stderr,"t2=%p\n",t2);
// timedata is N0 x N1 x ... x Nk real
// take the real data, fft it and transpose
// take a real chunk of data, fft it and place the output at correct intervals
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" );
kiss_fftr( st->cfg_r, timedata + k1*dimReal , tmp1 ); // tmp1 now holds dimReal/2+1 complex points
for (k2=0;k2<dimReal/2+1;++k2)
tmp2[ k2*dimOther+k1 ] = tmp1[k2];
}
// 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" );
for (k2=0;k2<dimReal/2+1;++k2) {
kiss_fftnd(st->cfg_nd, tmp2+k2*dimOther, tmp1); // tmp1 now holds dimOther complex points
for (k1=0;k1<dimOther;++k1)
freqdata[ k1*(dimReal/2+1) + k2] = tmp1[k1];
}
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);
int k1,k2;
int dimReal = st->dimReal;
int dimOther = st->dimOther;
kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf;
kiss_fft_cpx * tmp2;
if (dimReal/2+1 > dimOther)
tmp2 = tmp1+dimReal/2+1;
else
tmp2 = tmp1+dimOther;
for (k2=0;k2<dimReal/2+1;++k2) {
for (k1=0;k1<dimOther;++k1)
tmp1[k1] = freqdata[ k1*(dimReal/2+1) + k2 ];
kiss_fftnd(st->cfg_nd, tmp1, tmp2+k2*dimOther);
}
for (k1=0;k1<dimOther;++k1) {
for (k2=0;k2<dimReal/2+1;++k2)
tmp1[k2] = tmp2[ k2*dimOther+k1 ];
kiss_fftri( st->cfg_r,tmp1,timedata + k1*dimReal);
}
}