From b468fd96d71a151c1a7a2abfd1703fc3caf51e74 Mon Sep 17 00:00:00 2001
From: Mark Borgerding <mark@borgerding.net>
Date: Tue, 30 Dec 2003 15:18:41 +0000
Subject: [PATCH] Added mucho comments.

---
 tools/kfc.c        |  4 +--
 tools/kiss_fftnd.c | 68 ++++++++++++++++++++++++++++++++++++++++++++--
 2 files changed, 68 insertions(+), 4 deletions(-)

diff --git a/tools/kfc.c b/tools/kfc.c
index 02e2f5c..f3e4056 100644
--- a/tools/kfc.c
+++ b/tools/kfc.c
@@ -32,12 +32,12 @@ static const void * find_cached_fft(int nfft,int inverse)
     cached_fft * prev=NULL;
     while ( cur ) {
         if ( cur->nfft == nfft && inverse == cur->inverse )
-            break;//found the right node
+            break;/*found the right node*/
         prev = cur;
         cur = (cached_fft*)prev->next;
     }
     if (cur== NULL) {
-        // no cached node found, need to create a new one
+        /* no cached node found, need to create a new one*/
         kiss_fft_alloc(nfft,inverse,0,&len);
         cur = (cached_fft*)malloc(sizeof(cached_fft) + len );
         if (cur == NULL)
diff --git a/tools/kiss_fftnd.c b/tools/kiss_fftnd.c
index 4a3a292..81cb3de 100644
--- a/tools/kiss_fftnd.c
+++ b/tools/kiss_fftnd.c
@@ -1,3 +1,5 @@
+
+
 /*
 Copyright (c) 2003, Mark Borgerding
 
@@ -76,6 +78,68 @@ void * kiss_fftnd_alloc(int *dims,int ndims,int inverse_fft,void*mem,size_t*lenm
     return st;
 }
 
+/*
+ This works by tackling one dimension at a time.
+
+ In effect,
+ Each stage starts out by reshaping the matrix into a DixSi 2d matrix.
+ A Di-sized fft is taken of each column, transposing the matrix as it goes.
+
+Here's a 3-d example:
+Take a 2x3x4 matrix, laid out in memory as a contiguous buffer
+ [ [ [ a b c d ] [ e f g h ] [ i j k l ] ]
+   [ [ m n o p ] [ q r s t ] [ u v w x ] ] ]
+
+Stage 0 ( D=2): treat the buffer as a 2x12 matrix
+   [ [a b ... k l]
+     [m n ... w x] ]
+
+   FFT each column with size 2.
+   Transpose the matrix at the same time using kiss_fft_stride.
+
+   [ [ a+m a-m ]
+     [ b+n b-n]
+     ...
+     [ k+w k-w ]
+     [ l+x l-x ] ]
+
+   Note fft([x y]) == [x+y x-y]
+
+Stage 1 ( D=3) treats the buffer (the output of stage D=2) as an 3x8 matrix,
+   [ [ a+m a-m b+n b-n c+o c-o d+p d-p ] 
+     [ e+q e-q f+r f-r g+s g-s h+t h-t ]
+     [ i+u i-u j+v j-v k+w k-w l+x l-x ] ]
+
+   And perform FFTs (size=3) on each of the columns as above, transposing 
+   the matrix as it goes.  The output of stage 1 is 
+       (Legend: ap = [ a+m e+q i+u ]
+                am = [ a-m e-q i-u ] )
+   
+   [ [ sum(ap) fft(ap)[0] fft(ap)[1] ]
+     [ sum(am) fft(am)[0] fft(am)[1] ]
+     [ sum(bp) fft(bp)[0] fft(bp)[1] ]
+     [ sum(bm) fft(bm)[0] fft(bm)[1] ]
+     [ sum(cp) fft(cp)[0] fft(cp)[1] ]
+     [ sum(cm) fft(cm)[0] fft(cm)[1] ]
+     [ sum(dp) fft(dp)[0] fft(dp)[1] ]
+     [ sum(dm) fft(dm)[0] fft(dm)[1] ]  ]
+
+Stage 2 ( D=4) treats this buffer as a 4*6 matrix,
+   [ [ sum(ap) fft(ap)[0] fft(ap)[1] sum(am) fft(am)[0] fft(am)[1] ]
+     [ sum(bp) fft(bp)[0] fft(bp)[1] sum(bm) fft(bm)[0] fft(bm)[1] ]
+     [ sum(cp) fft(cp)[0] fft(cp)[1] sum(cm) fft(cm)[0] fft(cm)[1] ]
+     [ sum(dp) fft(dp)[0] fft(dp)[1] sum(dm) fft(dm)[0] fft(dm)[1] ]  ]
+
+   Then FFTs each column, transposing as it goes.
+
+   The resulting matrix is the 3d FFT of the 2x3x4 input matrix.
+
+   Note as a sanity check that the first element of the final 
+   stage's output (DC term) is 
+   sum( [ sum(ap) sum(bp) sum(cp) sum(dp) ] )
+   , i.e. the summation of all 24 input elements. 
+
+*/
 void kiss_fftnd(const void * cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
 {
     int i,k;
@@ -92,8 +156,8 @@ void kiss_fftnd(const void * cfg,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;
-        for (i=0;i<stride;++i) 
-            kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim,stride );
+        for ( i=0 ; i<stride ; ++i ) 
+            kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride );
 
         /*toggle back and forth between the two buffers*/
         if (bufout == st->tmpbuf){