From cb5312efdc48425d0d4fecf1a947a49933b12db8 Mon Sep 17 00:00:00 2001
From: Mark Borgerding <mark@borgerding.net>
Date: Thu, 6 Nov 2003 03:59:31 +0000
Subject: [PATCH] 2d fft seems to work

---
 kiss_fft.c      | 45 ++++++++++++++++++++++++++++-----------------
 kiss_fft.h      |  6 ++++++
 tools/fftutil.c | 26 ++++++++++++++++----------
 3 files changed, 50 insertions(+), 27 deletions(-)

diff --git a/kiss_fft.c b/kiss_fft.c
index c4cc0fd..152c1c5 100644
--- a/kiss_fft.c
+++ b/kiss_fft.c
@@ -36,6 +36,7 @@ typedef struct {
     int minus2; /*signify a 2-d transform*/
     kiss_fft_state * rowst;
     kiss_fft_state * colst;
+    kiss_fft_cpx * tmpbuf;
 }kiss_fft2d_state;
 
 
@@ -417,36 +418,46 @@ void * kiss_fft_alloc(int nfft,int inverse_fft)
 void * kiss_fft2d_alloc(int nrows,int ncols,int inverse_fft)
 {
     kiss_fft2d_state *st = NULL;
-    int size1,size2;
+    int size1,size2,sizetmp;
     size1 = allocsize(ncols);
     size2 = allocsize(nrows);
+    sizetmp = ncols > nrows ? ncols : nrows;
 
-    st = (kiss_fft2d_state *) malloc ( sizeof(kiss_fft2d_state) + size1 + size2 );
+    st = (kiss_fft2d_state *) malloc ( sizeof(kiss_fft2d_state) + size1 + size2 + sizetmp );
     if (!st)
         return NULL;
     
+    st->minus2 = -2;
     st->rowst = (kiss_fft_state *)(st+1); /*just beyond kiss_fft2d_state struct */
     st->colst = (kiss_fft_state *)( (char*)(st->rowst) + size1 );
+    st->tmpbuf = (kiss_fft_cpx *)( (char*)(st->rowst) + size1 + size2 );
     init_state (st->rowst, ncols, inverse_fft);
     init_state (st->colst, nrows, inverse_fft);
     return st;
 }
 
-void kiss_fft2d(const void * cfg,kiss_fft_cpx *f)
+void kiss_fft2d(const void * cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
 {
-    /* 
-     f is stored row-wise
-     * */
-    fprintf(stderr,"not yet implemented\n");
-    exit(1);
-}
+    /* input buffer fin is stored row-wise */
+    kiss_fft2d_state *st = ( kiss_fft2d_state *)cfg;
+    int row,col;
+    int nrows,ncols;
+    nrows = st->colst->nfft;
+    ncols = st->rowst->nfft;
 
-void kiss_fft2d_io(const void * cfg,const kiss_fft_cpx * fin,kiss_fft_cpx * fout)
-{
-    /*just use the in-place version sinc the multi-dim needs two passes anyway*/
-    kiss_fft2d_state *st = (kiss_fft2d_state *)cfg;
-    memcpy(fout,fin,sizeof(kiss_fft_cpx) * st->rowst->nfft * st->colst->nfft );
-    kiss_fft2d(cfg,fout);
+    /*fft each column*/
+    for (col=0;col<ncols;++col) {
+        for (row=0;row< nrows ;++row) 
+            st->tmpbuf[row] = fin[row*ncols + col];
+        kiss_fft(st->colst,st->tmpbuf);
+        for (row=0;row< nrows ;++row) {
+            fout[row*ncols + col] = st->tmpbuf[row];
+        }
+    }
+
+    /*fft each row */
+    for (row=0;row< nrows ;++row) 
+        kiss_fft(st->rowst , fout + row*ncols );
 }
 
 /* original form of processing function, first release of KISS FFT was in-place.  This maintains API. */
@@ -454,7 +465,7 @@ void kiss_fft(const void * cfg,kiss_fft_cpx *f)
 {
     const kiss_fft_state * st = cfg;
     if (st->nfft < 0) {
-        kiss_fft2d(cfg,f);
+        kiss_fft2d(cfg,f,f);
     }else{
         memcpy(st->tmpbuf,f,sizeof(kiss_fft_cpx)*st->nfft);
         fft_work( f, st->tmpbuf, 1, st->factors,st );
@@ -466,7 +477,7 @@ void kiss_fft_io(const void * cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
 {
     const kiss_fft_state * st = cfg;
     if (st->nfft < 0) {
-        kiss_fft2d_io(cfg,fin,fout);
+        kiss_fft2d(cfg,fin,fout);
     }else{
         fft_work( fout, fin, 1, st->factors,st );
     }
diff --git a/kiss_fft.h b/kiss_fft.h
index ebeb246..e20a94f 100644
--- a/kiss_fft.h
+++ b/kiss_fft.h
@@ -42,6 +42,12 @@ void kiss_fft( const void* cfg_from_alloc , kiss_fft_cpx *f ); /* call for each
 /* two buffer version */
 void kiss_fft_io(const void * cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);
 
+/* allocate a 2-dimensional FFT
+   kiss_fft() is used as in the 1d case, but the data should be stored rowwise,
+   in other words, an array made up of row[0], then row[1], etc
+ */
+void * kiss_fft2d_alloc(int nrows,int ncols,int inverse_fft);
+
 /* when done with the cfg for a given fft size and direction, simply free it*/
 #define kiss_fft_free free
 
diff --git a/tools/fftutil.c b/tools/fftutil.c
index 2dd3420..edaa8b6 100644
--- a/tools/fftutil.c
+++ b/tools/fftutil.c
@@ -26,26 +26,30 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 
 
 
-void fft_file(FILE * fin,FILE * fout,int nfft,int isinverse,int useascii,int times)
+void fft_file(FILE * fin,FILE * fout,int nfft,int nrows,int isinverse,int useascii,int times)
 {
     int i;
     void *st;
     kiss_fft_cpx * buf;
     kiss_fft_cpx * bufout;
+            
 
-    buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
-    bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
-    st = kiss_fft_alloc( nfft ,isinverse );
+    buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft *nrows );
+    bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft *nrows);
+    if (nrows!=1)
+        st = kiss_fft2d_alloc( nrows,nfft ,isinverse );
+    else
+        st = kiss_fft_alloc( nfft ,isinverse );
 
-    while ( fread( buf , sizeof(kiss_fft_cpx) * nfft ,1, fin ) > 0 ) {
+    while ( fread( buf , sizeof(kiss_fft_cpx) * nfft * nrows ,1, fin ) > 0 ) {
         for (i=0;i<times;++i)
             kiss_fft_io( st , buf ,bufout);
         if (useascii) {
             int i;
-            for (i=0;i<nfft;++i) 
+            for (i=0;i<nfft*nrows;++i) 
                 fprintf(fout, "(%g,%g) ", (double)bufout[i].r,(double)bufout[i].i);
-        }else{            
-            fwrite( bufout , sizeof(kiss_fft_cpx) , nfft , fout );
+        }else{
+            fwrite( bufout , sizeof(kiss_fft_cpx) , nfft*nrows , fout );
         }
     }
     free(st);
@@ -61,13 +65,15 @@ int main(int argc,char ** argv)
     FILE *fout=stdout;
     int useascii=0;
     int times=1;
+    int nrows=1;
 
     while (1) {
-        int c=getopt(argc,argv,"n:iax:");
+        int c=getopt(argc,argv,"n:iax:c:");
         if (c==-1) break;
         switch (c) {
             case 'a':useascii=1;break;
             case 'n':nfft = atoi(optarg);break;
+            case 'c':nrows = atoi(optarg);break;
             case 'i':isinverse=1;break;
             case 'x':times=atoi(optarg);break;
         }
@@ -85,7 +91,7 @@ int main(int argc,char ** argv)
         ++optind;
     }
 
-    fft_file(fin,fout,nfft,isinverse,useascii,times);
+    fft_file(fin,fout,nfft,nrows,isinverse,useascii,times);
 
     if (fout!=stdout) fclose(fout);
     if (fin!=stdin) fclose(fin);