summaryrefslogtreecommitdiff
path: root/src/signalProcessing/fft
diff options
context:
space:
mode:
Diffstat (limited to 'src/signalProcessing/fft')
-rw-r--r--src/signalProcessing/fft/Makefile.am20
-rw-r--r--src/signalProcessing/fft/Makefile.in37
-rw-r--r--src/signalProcessing/fft/dfft2.c4
-rw-r--r--src/signalProcessing/fft/dfftbi.c39
-rw-r--r--src/signalProcessing/fft/fft842.c266
-rw-r--r--src/signalProcessing/fft/fft_internal.h16
-rw-r--r--src/signalProcessing/fft/r2tx.c42
-rw-r--r--src/signalProcessing/fft/r4tx.c80
-rw-r--r--src/signalProcessing/fft/r8tx.c438
-rw-r--r--src/signalProcessing/fft/testDoubleFft.c279
-rw-r--r--src/signalProcessing/fft/testFloatFft.c36
-rw-r--r--src/signalProcessing/fft/zfftma.c143
12 files changed, 1050 insertions, 350 deletions
diff --git a/src/signalProcessing/fft/Makefile.am b/src/signalProcessing/fft/Makefile.am
index 657a39dc..a67465f9 100644
--- a/src/signalProcessing/fft/Makefile.am
+++ b/src/signalProcessing/fft/Makefile.am
@@ -15,6 +15,7 @@
libFft_la_CFLAGS = -I . \
-I $(top_builddir)/type \
-I $(top_builddir)/signalProcessing/includes \
+ -I $(top_builddir)/operations/includes \
-I $(top_builddir)/auxiliaryFunctions/includes
instdir = $(top_builddir)/lib
@@ -25,7 +26,8 @@ libFft_la_SOURCES = $(HEAD) $(SRC)
HEAD = ../includes/fft.h
-SRC = dfft2.c \
+SRC = zfftma.c \
+ dfft2.c \
dfftbi.c \
dfftmx.c \
fft842.c \
@@ -38,12 +40,20 @@ SRC = dfft2.c \
# Checking Part
####
-check_INCLUDES = -I $(top_builddir)/signalProcessing/includes \
- -I $(top_builddir)/type
+check_INCLUDES = -I . \
+ -I $(top_builddir)/type \
+ -I $(top_builddir)/signalProcessing/includes \
+ -I $(top_builddir)/operations/includes \
+ -I $(top_builddir)/auxiliaryFunctions/includes
+
check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
$(top_builddir)/type/libFloatComplex.la \
- $(top_builddir)/elementaryFunctions/fft/libFft.la \
+ $(top_builddir)/lib/lapack/libscilapack.la \
+ $(top_builddir)/lib/blas/libsciblas.la \
+ $(top_builddir)/signalProcessing/fft/libFft.la \
+ $(top_builddir)/operations/addition/libAddition.la \
+ $(top_builddir)/operations/subtraction/libSubtraction.la \
@LIBMATH@
check_PROGRAMS = testFloatFft testDoubleFft
@@ -59,4 +69,4 @@ testFloatFft_LDADD = $(check_LDADD)
testDoubleFft_SOURCES = testDoubleFft.c
testDoubleFft_CFLAGS = $(check_INCLUDES)
-testDoubleFft_LDADD = $(check_LDADD) \ No newline at end of file
+testDoubleFft_LDADD = $(check_LDADD)
diff --git a/src/signalProcessing/fft/Makefile.in b/src/signalProcessing/fft/Makefile.in
index 781a5253..672e6358 100644
--- a/src/signalProcessing/fft/Makefile.in
+++ b/src/signalProcessing/fft/Makefile.in
@@ -54,9 +54,9 @@ pkglibLTLIBRARIES_INSTALL = $(INSTALL)
LTLIBRARIES = $(pkglib_LTLIBRARIES)
libFft_la_LIBADD =
am__objects_1 =
-am__objects_2 = libFft_la-dfft2.lo libFft_la-dfftbi.lo \
- libFft_la-dfftmx.lo libFft_la-fft842.lo libFft_la-r2tx.lo \
- libFft_la-r4tx.lo libFft_la-r8tx.lo
+am__objects_2 = libFft_la-zfftma.lo libFft_la-dfft2.lo \
+ libFft_la-dfftbi.lo libFft_la-dfftmx.lo libFft_la-fft842.lo \
+ libFft_la-r2tx.lo libFft_la-r4tx.lo libFft_la-r8tx.lo
am_libFft_la_OBJECTS = $(am__objects_1) $(am__objects_2)
libFft_la_OBJECTS = $(am_libFft_la_OBJECTS)
libFft_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
@@ -66,7 +66,11 @@ am_testDoubleFft_OBJECTS = testDoubleFft-testDoubleFft.$(OBJEXT)
testDoubleFft_OBJECTS = $(am_testDoubleFft_OBJECTS)
am__DEPENDENCIES_1 = $(top_builddir)/type/libDoubleComplex.la \
$(top_builddir)/type/libFloatComplex.la \
- $(top_builddir)/elementaryFunctions/fft/libFft.la
+ $(top_builddir)/lib/lapack/libscilapack.la \
+ $(top_builddir)/lib/blas/libsciblas.la \
+ $(top_builddir)/signalProcessing/fft/libFft.la \
+ $(top_builddir)/operations/addition/libAddition.la \
+ $(top_builddir)/operations/subtraction/libSubtraction.la
testDoubleFft_DEPENDENCIES = $(am__DEPENDENCIES_1)
testDoubleFft_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(CCLD) $(testDoubleFft_CFLAGS) \
@@ -209,13 +213,15 @@ top_srcdir = @top_srcdir@
libFft_la_CFLAGS = -I . \
-I $(top_builddir)/type \
-I $(top_builddir)/signalProcessing/includes \
+ -I $(top_builddir)/operations/includes \
-I $(top_builddir)/auxiliaryFunctions/includes
instdir = $(top_builddir)/lib
pkglib_LTLIBRARIES = libFft.la
libFft_la_SOURCES = $(HEAD) $(SRC)
HEAD = ../includes/fft.h
-SRC = dfft2.c \
+SRC = zfftma.c \
+ dfft2.c \
dfftbi.c \
dfftmx.c \
fft842.c \
@@ -227,12 +233,19 @@ SRC = dfft2.c \
####
# Checking Part
####
-check_INCLUDES = -I $(top_builddir)/signalProcessing/includes \
- -I $(top_builddir)/type
+check_INCLUDES = -I . \
+ -I $(top_builddir)/type \
+ -I $(top_builddir)/signalProcessing/includes \
+ -I $(top_builddir)/operations/includes \
+ -I $(top_builddir)/auxiliaryFunctions/includes
check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
$(top_builddir)/type/libFloatComplex.la \
- $(top_builddir)/elementaryFunctions/fft/libFft.la \
+ $(top_builddir)/lib/lapack/libscilapack.la \
+ $(top_builddir)/lib/blas/libsciblas.la \
+ $(top_builddir)/signalProcessing/fft/libFft.la \
+ $(top_builddir)/operations/addition/libAddition.la \
+ $(top_builddir)/operations/subtraction/libSubtraction.la \
@LIBMATH@
@@ -334,6 +347,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libFft_la-r2tx.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libFft_la-r4tx.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libFft_la-r8tx.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libFft_la-zfftma.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testDoubleFft-testDoubleFft.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testFloatFft-testFloatFft.Po@am__quote@
@@ -358,6 +372,13 @@ distclean-compile:
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(LTCOMPILE) -c -o $@ $<
+libFft_la-zfftma.lo: zfftma.c
+@am__fastdepCC_TRUE@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libFft_la_CFLAGS) $(CFLAGS) -MT libFft_la-zfftma.lo -MD -MP -MF $(DEPDIR)/libFft_la-zfftma.Tpo -c -o libFft_la-zfftma.lo `test -f 'zfftma.c' || echo '$(srcdir)/'`zfftma.c
+@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/libFft_la-zfftma.Tpo $(DEPDIR)/libFft_la-zfftma.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='zfftma.c' object='libFft_la-zfftma.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libFft_la_CFLAGS) $(CFLAGS) -c -o libFft_la-zfftma.lo `test -f 'zfftma.c' || echo '$(srcdir)/'`zfftma.c
+
libFft_la-dfft2.lo: dfft2.c
@am__fastdepCC_TRUE@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libFft_la_CFLAGS) $(CFLAGS) -MT libFft_la-dfft2.lo -MD -MP -MF $(DEPDIR)/libFft_la-dfft2.Tpo -c -o libFft_la-dfft2.lo `test -f 'dfft2.c' || echo '$(srcdir)/'`dfft2.c
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/libFft_la-dfft2.Tpo $(DEPDIR)/libFft_la-dfft2.Plo
diff --git a/src/signalProcessing/fft/dfft2.c b/src/signalProcessing/fft/dfft2.c
index 20011c04..7a5ac88a 100644
--- a/src/signalProcessing/fft/dfft2.c
+++ b/src/signalProcessing/fft/dfft2.c
@@ -33,7 +33,9 @@ void dfft2 ( double* a , double* b , int nseg , int n , int nspn , int isn , int
iw[3] = lw ;
iw[4] = 10 ;
- dfftbi ( a , b , nseg , n , nspn , isn , ierr ,iIw[0], iIw[1], iIw[2], iIw[3], iIw[4], iw, iIw);
+ dfftbi ( a , b , nseg , n , nspn ,
+ isn , ierr , iIw[0], iIw[1] , iIw[2],
+ iIw[3], iIw[4], iw, iIw );
diff --git a/src/signalProcessing/fft/dfftbi.c b/src/signalProcessing/fft/dfftbi.c
index de672e82..8d3f7846 100644
--- a/src/signalProcessing/fft/dfftbi.c
+++ b/src/signalProcessing/fft/dfftbi.c
@@ -11,10 +11,13 @@
*/
#include <stdlib.h>
+#include <stdio.h>
#include "max.h"
#include "fft_internal.h"
-void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , int ierr, int lout , int lnow , int lused , int lmax , int lbook , double* rstak , int* istak )
+void dfftbi ( double* a , double* b , int nseg , int n , int nspn ,
+ int isn , int ierr , int lout , int lnow , int lused ,
+ int lmax , int lbook , double* rstak , int* istak )
{
int nfac[15] ;
@@ -31,7 +34,7 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
int nspan ;
int nitems ;
int ntot ;
- int maxp ;
+ int maxp = 0;
int maxf ;
int itype;
int istkgt ;
@@ -40,13 +43,15 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
int nf = abs ( n ) ;
ierr = 0 ;
-
+ printf ( "debut de dfftbi \n" );
/*determine the factors of n */
if ( nf == 1)
+ {
+ printf ( "argh return 1 \n" );
return ;
-
+ }
k = nf ;
nspan = abs ( nf*nspn ) ;
@@ -55,16 +60,18 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
if ( isn*ntot == 0 )
{
ierr = 1 ;
+ printf ( "argh return 2 \n" );
return ;
}
- while ( (k & 15) == 0 )
+ do
{
m++;
- nfac[m-1] = j ;
+ printf ("m %d ,k %d ,k2 %d\n" , m , k ,(int) (k/16)*16 );
+ nfac[m-1] = 4 ;
k = k >> 4 ;
- }
+ }while ( (k- (int)(k/16)*16 ) == 0 );
do
{
@@ -81,7 +88,7 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
}while ( jj <= k);
-
+printf ( "ploa\n" );
if ( k <= 4)
{
kt = m;
@@ -102,7 +109,7 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
kt = m ;
maxp = max ( (kt+1)*2 , k-1);
j=2;
-
+printf ( "plob\n" );
do
{
if ( k%j == 0 )
@@ -122,9 +129,11 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
maxp = m + kt + 1 ;
if ( m + kt > 15)
+ {
ierr = 2 ;
+ printf ( "argh return 5 \n" );
return ;
-
+ }
if ( kt != 0 )
{
j = kt ;
@@ -144,16 +153,18 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
for ( kkk = 1 ; kkk < m ; kkk++ )
maxf = max ( maxf , nfac[kkk-1]);
- nitems = maxf * 4 ;
+ nitems = maxf * 4 ;
itype = 4 ;
istkgt = ( lnow*isize[1] -1)/isize[itype-1] + 2;
i = ( (istkgt - 1 + nitems) * isize[itype-1] -1) / isize[1] + 3 ;
-
+ printf ("i %d ,\n lmax %d\n istkgt %d\n lnow %d \n", i , lmax , istkgt , lnow ) ;
+
if ( i > lmax )
{
ierr = -i ;
+ printf ( "argh return 3 \n" );
return ;
}
@@ -178,6 +189,7 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int isn , in
if ( i > lmax )
{
ierr = -i ;
+ printf ( "argh return 4 \n" );
return ;
}
@@ -197,6 +209,7 @@ c k=2*k-1
c *********************************************
*/
+ printf ( "dfftmx me voilà tayoooooooo \n" );
dfftmx( a , b , ntot , nf , nspan , isn , m , kt , &rstak[j-1] , &rstak[jj-1] , &rstak[j2-1] , &rstak[j3-1] , &istak[k-1] , nfac);
k =2 ;
@@ -206,6 +219,7 @@ c *********************************************
if (!( lbook <= lnow && lnow <= lused && lused <= lmax ))
{
ierr = 3 ;
+ printf ( "argh return 6 \n" );
return ;
}
@@ -220,5 +234,6 @@ c *********************************************
lnow = istak[lnow-1] ;
in-- ;
}
+ printf ( "fin de dfftbi \n" );
return ;
}
diff --git a/src/signalProcessing/fft/fft842.c b/src/signalProcessing/fft/fft842.c
index ba287a54..0b9f1dbc 100644
--- a/src/signalProcessing/fft/fft842.c
+++ b/src/signalProcessing/fft/fft842.c
@@ -11,138 +11,180 @@
*/
#include "fft_internal.h"
+#include <stdio.h>
+
-void fft842 ( int _iDirect , int _iDimen , double* _pdblReal , double* _pdblImag , int _err )
-{
- int i = 0 ;
- int ipass = 1 ;
+/* get binary log of integer argument; exact if n a power of 2 */
+static int fastlog2( int n)
+{
+ int log = -1;
+ while(n) {
+ log++;
+ n >>= 1;
+ }
+ return(log);
+}
- int ij ;
- int ji ;
- int j1, j2, j3, j4, j5, j6, j7, j8, j9, j10, j11, j12, j13, j14;
- int iPow2ofDimen ;
- int iTempDimen = 1 ;
- int iTemp2Pow2ofDimen ;
- int iPow8ofDimen ;
- int iLengt ;
- int CRES = 0 ;
+/*
+ int in; FORWARD or INVERSE
+ int n; length of vector
+ DPCOMPLEX *b; input vector
+*/
+void fft842 (doubleComplex* b, int size , int in)
+{
+ double fn;
+ doubleComplex temp ;
- int l[15] ;
+ int L[16],L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15;
+ int j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14;
+ int i, j, ij, ji, ij1, ji1;
+ int n2pow, n8pow, nthpo, ipass, nxtlt, lengt;
+ n2pow = fastlog2( size );
+ nthpo = size ;
+ fn = 1.0 / (double)nthpo;
- double dblTemp ;
+ if(in==FORWARD)
+ /* take conjugate */
+ for(i=0;i< size ;i++) {
+ b[i] = DoubleComplex ( zreals( b[i]) , - zimags (b[i]));
- do
- {
- iPow2ofDimen = i ;
- if ( i != 0)
- iTempDimen *= 2 ;
- i++;
- }while ( _iDimen != iTempDimen || i < 15 );
+ /* b[i].im *= -1.0;*/
+ }
- if ( _iDimen != iTempDimen)
+ if(in==INVERSE)
+ /*scramble inputs*/
+ for(i=0,j=size/2;j<size;i++,j++)
{
- _err = 1 ;
- return ;
+ temp = DoubleComplex ( zreals ( b[j] ) , zimags( b[j] ));
+ b[j] = DoubleComplex ( zreals ( b[i] ) , zimags( b[i] ));
+ b[i] = DoubleComplex ( zreals ( temp ) , zimags( temp ));
+
+ /*
+ r = b[j].re; fi = b[j].im;
+ b[j].re = b[i].re; b[j].im = b[i].im;
+ b[i].re = r; b[i].im = fi;
+ */
}
- _err = 0 ;
+ n8pow = n2pow/3;
+
+ if(n8pow)
+ {
+ /* radix 8 iterations */
+ for(ipass=1;ipass<=n8pow;ipass++)
+ {
+ nxtlt = 0x1 << (n2pow - 3*ipass);
+ lengt = 8*nxtlt;
+ printf ( "on appelle r%dtx \n" , 8);
+ r8tx(nxtlt,nthpo,lengt,
+ b,b+nxtlt,b+2*nxtlt,
+ b+3*nxtlt,b+4*nxtlt,b+5*nxtlt,
+ b+6*nxtlt,b+7*nxtlt);
+ }
+ }
+
+ if(n2pow%3 == 1)
+ {
+ /* radix 2 iteration needed */
+ printf ( "on appelle r%dtx \n" , 2);
+ r2tx(nthpo,b,b+1);
+ }
+
+
+ if(n2pow%3 == 2)
+ {
+ /* radix 4 iteration needed */
+ printf ( "on appelle r%dtx \n" , 4);
+ r4tx(nthpo,b,b+1,b+2,b+3);
+ }
+
+
+
+ for(j=1;j<=15;j++)
+ {
+ L[j] = 1;
+ if(j-n2pow <= 0) L[j] = 0x1 << (n2pow + 1 - j);
+ }
+ L15=L[1];L14=L[2];L13=L[3];L12=L[4];L11=L[5];L10=L[6];L9=L[7];
+ L8=L[8];L7=L[9];L6=L[10];L5=L[11];L4=L[12];L3=L[13];L2=L[14];L1=L[15];
+
+ ij = 1;
+
+
+ for(j1=1;j1<=L1;j1++)
+ for(j2=j1;j2<=L2;j2+=L1)
+ for(j3=j2;j3<=L3;j3+=L2)
+ for(j4=j3;j4<=L4;j4+=L3)
+ for(j5=j4;j5<=L5;j5+=L4)
+ for(j6=j5;j6<=L6;j6+=L5)
+ for(j7=j6;j7<=L7;j7+=L6)
+ for(j8=j7;j8<=L8;j8+=L7)
+ for(j9=j8;j9<=L9;j9+=L8)
+ for(j10=j9;j10<=L10;j10+=L9)
+ for(j11=j10;j11<=L11;j11+=L10)
+ for(j12=j11;j12<=L12;j12+=L11)
+ for(j13=j12;j13<=L13;j13+=L12)
+ for(j14=j13;j14<=L14;j14+=L13)
+ for(ji=j14;ji<=L15;ji+=L14)
+ {
+ ij1 = ij-1;
+ ji1 = ji-1;
+
+ if(ij-ji<0)
+ {
+ temp = b[ij1];
+ b[ij1] = b[ji1];
+ b[ji1] = temp;
+
+ /*
+ r = b[ij1].re;
+ b[ij1].re = b[ji1].re;
+ b[ji1].re = r;
+ fi = b[ij1].im;
+ b[ij1].im = b[ji1].im;
+ b[ji1].im = fi;
+ */
+ }
+ ij++;
+ }
+
+ if(in==FORWARD) /* take conjugates & unscramble outputs */
+ for(i=0,j=size/2;j<size;i++,j++)
+ {
+ temp = DoubleComplex ( zreals ( b[j] ) ,- zimags( b[j] ));
+ b[j] = DoubleComplex ( zreals ( b[i] ) ,- zimags( b[i] ));
+ b[i] = DoubleComplex ( zreals ( temp ) , zimags( temp ));
+
+
+ /* r = b[j].re; fi = b[j].im;
+ b[j].re = b[i].re; b[j].im = -b[i].im;
+ b[i].re = r; b[i].im = -fi;
+ */
+ }
+
+
+ if(in==INVERSE) /* scale outputs */
+ for(i=0;i<nthpo;i++)
+ {
- if ( _iDirect == 0 )
- for ( i = 0 ; i < _iDimen ; i++)
- {
- _pdblImag[i] *= -1 ;
- }
+ b[i] = DoubleComplex ( zreals ( b[i] )*fn , zimags( b[i] )*fn);
- iPow8ofDimen = iPow2ofDimen / 3 ;
- if ( iPow8ofDimen != 0 )
- {
- /* development of the algo in 8-base if ... */
- for ( ipass = 1 ; i < iPow8ofDimen ; i++ )
- {
- iTemp2Pow2ofDimen = 1 << ( iPow2ofDimen - 3*ipass ) ;
- iLengt = iTemp2Pow2ofDimen * 8 ;
- r8tx( iTemp2Pow2ofDimen, _iDimen, iLengt, _pdblReal, _pdblImag );
- }
}
- CRES = iPow2ofDimen - 3*iPow8ofDimen - 1 ;
- if ( CRES >= 0)
- {
- if ( CRES == 0 )
- {
- r2tx ( _iDimen , _pdblReal , _pdblImag) ;
-
- }
- else
- {
- r4tx ( _iDimen , _pdblReal , _pdblImag) ;
- }
- }
- for ( i = 14 ; i >= 0 ; i-- )
- {
- l[i] = 1 ;
- CRES = i - iPow2ofDimen ;
- if ( CRES <= 0)
- l[i] = 1 << ( iPow2ofDimen + 1 -i ) ;
+ for ( i = 0 ; i < size /2 ; i++)
+ {
+ temp = b[i] ;
+ b[i] = b[i+(size/2)];
+ b[i+(size/2)]= temp ;
- }
- ij = 0 ;
-
-
- for ( j1=0 ; j1 < l[1-1] ; j1++ )
- for ( j2=j1 ; j2 < l[2-1] ; j2 += l[1-1] )
- for ( j3=j2 ; j3 < l[3-1] ; j3 += l[2-1] )
- for ( j4=j3 ; j4 < l[4-1] ; j4 += l[3-1] )
- for ( j5=j4 ; j5 < l[5-1] ; j5 += l[4-1] )
- for ( j6=j5 ; j6 < l[6-1] ; j6 += l[5-1] )
- for ( j7=j6 ; j7 < l[7-1] ; j7 += l[6-1] )
- for ( j8=j7 ; j8 < l[8-1] ; j8 += l[7-1] )
- for ( j9=j8 ; j9 < l[9-1] ; j9 += l[8-1] )
- for ( j10=j9 ; j10 < l[10-1] ; j10 += l[9-1] )
- for ( j11=j10 ; j11 < l[11-1] ; j11 += l[10-1] )
- for ( j12=j11 ; j12 < l[12-1] ; j12 += l[11-1] )
- for ( j13=j12 ; j13 < l[13-1] ; j13 += l[12-1] )
- for ( j14=j13 ; j14 < l[14-1] ; j14 += l[13-1] )
- for ( ji=j14 ; ji < l[15-1] ; ji += l[14-1] )
- {
- CRES = ij - ji ;
- if ( CRES < 0 )
- {
- dblTemp = _pdblReal[ij];
- _pdblReal[ij] = _pdblReal[ji];
- _pdblReal[ji] = dblTemp;
- dblTemp = _pdblImag[ij];
- _pdblImag[ij] = _pdblImag[ji];
- _pdblImag[ji] = dblTemp;
-
- }
- ij ++ ;
- }
-
-
-
- /*130*/
- if ( _iDirect == 0 )
- {
- for ( i = 0 ; i < _iDimen ; i++)
- {
- _pdblImag[i] *= -1 ;
- }
- }
- else
- {
- for ( i = 0 ; i < _iDimen ; i++)
- {
- _pdblReal[i] /= _iDimen ;
- _pdblImag[i] /= _iDimen ;
- }
- }
+ }
+
}
diff --git a/src/signalProcessing/fft/fft_internal.h b/src/signalProcessing/fft/fft_internal.h
index 1b374872..f230a003 100644
--- a/src/signalProcessing/fft/fft_internal.h
+++ b/src/signalProcessing/fft/fft_internal.h
@@ -13,6 +13,11 @@
#ifndef __FFT_INTERNAL_H__
#define __FFT_INTERNAL_H__
+#include "addition.h"
+#include "subtraction.h"
+
+#define FORWARD 0
+#define INVERSE 1
void dfft2 ( double* a , double* b , int nseg , int n , int nspn ,
int isn , int ierr, double* iw , int lw );
@@ -22,12 +27,13 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn ,
int isn , int ierr, int lout , int lnow ,
int lused ,int lmax , int lbook , double* rstak , int* istakk );
-void fft842 ( int _iDirect , int _iDimen , double* _pdblReal , double* _pdblImag , int _err );
-
-void r2tx ( int _iDimen , double* _pdblReal, double* _pdblImag );
-void r4tx ( int _iDimen , double* _pdblReal, double* _pdblImag) ;
-void r8tx ( int _iTempDimen , int _iDimen , int _iLengt , double* _pdblReal, double* _pdblImag );
+void fft842 (doubleComplex* b, int size , int in);
+void r2tx(int nthpo, doubleComplex* c0, doubleComplex* c1);
+void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, doubleComplex* c3);
+void r8tx ( int nxtlt,int nthpo,int lengt,
+ doubleComplex* cc0,doubleComplex* cc1,doubleComplex* cc2,doubleComplex* cc3,
+ doubleComplex* cc4,doubleComplex* cc5,doubleComplex* cc6,doubleComplex* cc7);
int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan,
int _iIsn, int _iM, int _iKt, double* _pdblWt, double* _pdblCk,
diff --git a/src/signalProcessing/fft/r2tx.c b/src/signalProcessing/fft/r2tx.c
index fd44c19c..c1231dc9 100644
--- a/src/signalProcessing/fft/r2tx.c
+++ b/src/signalProcessing/fft/r2tx.c
@@ -11,25 +11,37 @@
*/
#include "fft_internal.h"
-
-void r2tx ( int _iDimen , double* _pdblReal, double* _pdblImag )
+#include <stdio.h>
+/*
+** radix 2 iteration subroutine
+*/
+void r2tx(int nthpo, doubleComplex* c0, doubleComplex* c1)
{
- double dblTemp ;
+ int k,kk;
+ /* double *cr0, *ci0, *cr1, *ci1, r1, fi1;*/
+ doubleComplex temp ;
+/* cr0 = &(c0[0].re);
+ ci0 = &(c0[0].im);
+ cr1 = &(c1[0].re);
+ ci1 = &(c1[0].im);*/
- int i = 0 ;
+ for(k=0;k<nthpo;k+=2)
+ {
+ kk = k*2;
- for ( i = 0 ; i < _iDimen ; i += 2 )
- {
- /*for real part */
- dblTemp = _pdblReal[i] + _pdblReal[i+1] ;
- _pdblReal[i+1] = _pdblReal[i] - _pdblReal[i+1] ;
- _pdblReal[i+0] = dblTemp ;
- /*for imaginary one */
- dblTemp = _pdblImag[i] + _pdblImag[i+1] ;
- _pdblImag[i+1] = _pdblImag[i] - _pdblImag[i+1] ;
- _pdblImag[i] = dblTemp ;
+ temp = zadds ( c0[kk] , c1[kk] );
+ c1[kk] = zdiffs( c0[kk] , c1[kk] );
+ c0[kk] = DoubleComplex ( zreals ( temp ) , zimags( temp ));
+/*
+ r1 = cr0[kk] + cr1[kk];
+ cr1[kk] = cr0[kk] - cr1[kk];
+ cr0[kk] = r1;
- }
+ fi1 = ci0[kk] + ci1[kk];
+ ci1[kk] = ci0[kk] - ci1[kk];
+ ci0[kk] = fi1;
+*/
+ }
}
diff --git a/src/signalProcessing/fft/r4tx.c b/src/signalProcessing/fft/r4tx.c
index ed3dd896..c75b5d4e 100644
--- a/src/signalProcessing/fft/r4tx.c
+++ b/src/signalProcessing/fft/r4tx.c
@@ -12,37 +12,53 @@
#include "fft_internal.h"
-void r4tx ( int _iDimen , double* _pdblReal, double* _pdblImag)
+
+
+/*
+** radix 4 iteration subroutine
+*/
+void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, doubleComplex* c3)
{
- double dblTemp1 ;
- double dblTemp2 ;
- double dblTemp3 ;
- double dblTemp4 ;
-
- int i = 0 ;
-
- for ( i = 0 ; i < _iDimen ; i += 4 )
- {
- /*for real part */
- dblTemp1 = _pdblReal[i+0] + _pdblReal[i+2] ;
- dblTemp2 = _pdblReal[i+0] - _pdblReal[i+2] ;
- dblTemp3 = _pdblReal[i+1] + _pdblReal[i+3] ;
- dblTemp4 = _pdblReal[i+1] - _pdblReal[i+3] ;
-
- _pdblReal[i+0] = dblTemp1 + dblTemp3;
- _pdblReal[i+0] = dblTemp1 - dblTemp3;
- _pdblReal[i+2] = dblTemp2 + dblTemp4;
- _pdblReal[i+3] = dblTemp2 - dblTemp4;
-
- /*for imaginary part */
- dblTemp1 = _pdblImag[i+0] + _pdblImag[i+2] ;
- dblTemp2 = _pdblImag[i+0] - _pdblImag[i+2] ;
- dblTemp3 = _pdblImag[i+1] + _pdblImag[i+3] ;
- dblTemp4 = _pdblImag[i+1] - _pdblImag[i+3] ;
-
- _pdblImag[i+0] = dblTemp1 + dblTemp3;
- _pdblImag[i+1] = dblTemp1 - dblTemp3;
- _pdblImag[i+2] = dblTemp2 + dblTemp4;
- _pdblImag[i+3] = dblTemp2 - dblTemp4;
- }
+ int k,kk;
+ doubleComplex temp1 , temp2 , temp3 , temp4 ;
+
+ for(k=1;k<=nthpo;k+=4)
+ {
+ kk = (k-1)*2; /* real and imag parts alternate */
+
+
+ temp1 = zadds ( c0[kk] , c2[kk] ) ;
+ temp2 = zdiffs( c0[kk] , c2[kk] ) ;
+ temp3 = zadds ( c1[kk] , c3[kk] ) ;
+ temp4 = zdiffs( c1[kk] , c3[kk] ) ;
+/*
+ r1 = cr0[kk] + cr2[kk];
+ r2 = cr0[kk] - cr2[kk];
+ r3 = cr1[kk] + cr3[kk];
+ r4 = cr1[kk] - cr3[kk];
+
+
+ i1 = ci0[kk] + ci2[kk];
+ i2 = ci0[kk] - ci2[kk];
+ i3 = ci1[kk] + ci3[kk];
+ i4 = ci1[kk] - ci3[kk];
+*/
+ c0[kk] = zadds ( temp1 , temp3 );
+ c1[kk] = zdiffs( temp1 , temp3 );
+
+/*
+ cr0[kk] = r1 + r3;
+ ci0[kk] = i1 + i3;
+ cr1[kk] = r1 - r3;
+ ci1[kk] = i1 - i3;
+*/
+ c2[kk] = DoubleComplex ( zreals ( temp2 ) - zimags( temp4 ) , zimags ( temp2 ) + zreals( temp4 ) );
+ c3[kk] = DoubleComplex ( zreals ( temp2 ) + zimags( temp4 ) , zimags ( temp2 ) - zreals( temp4 ) );
+/*
+ cr2[kk] = r2 - i4;
+ ci2[kk] = i2 + r4;
+ cr3[kk] = r2 + i4;
+ ci3[kk] = i2 - r4;
+*/
+ }
}
diff --git a/src/signalProcessing/fft/r8tx.c b/src/signalProcessing/fft/r8tx.c
index 1d109de3..dd693134 100644
--- a/src/signalProcessing/fft/r8tx.c
+++ b/src/signalProcessing/fft/r8tx.c
@@ -1,3 +1,6 @@
+
+
+
/*
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2008 - INRIA - Allan SIMON
@@ -12,169 +15,284 @@
#include "fft_internal.h"
#include <math.h>
+#include <stdio.h>
+
+
+/*
+** radix 8 iteration subroutine
+*/
+void r8tx ( int nxtlt,int nthpo,int lengt,
+ doubleComplex* cc0,doubleComplex* cc1,doubleComplex* cc2,doubleComplex* cc3,
+ doubleComplex* cc4,doubleComplex* cc5,doubleComplex* cc6,doubleComplex* cc7)
-void r8tx ( int _iTempDimen , int _iDimen , int _iLengt , double* _pdblReal, double* _pdblImag )
{
- int j = 0 ;
- int i = 0 ;
- int k = 0 ;
-
-
- double* pdblReal0 = _pdblReal + ( 0*_iTempDimen );
- double* pdblReal1 = _pdblReal + ( 1*_iTempDimen );
- double* pdblReal2 = _pdblReal + ( 2*_iTempDimen );
- double* pdblReal3 = _pdblReal + ( 3*_iTempDimen );
- double* pdblReal4 = _pdblReal + ( 4*_iTempDimen );
- double* pdblReal5 = _pdblReal + ( 5*_iTempDimen );
- double* pdblReal6 = _pdblReal + ( 6*_iTempDimen );
- double* pdblReal7 = _pdblReal + ( 7*_iTempDimen );
-
- double* pdblImag0 = _pdblImag + ( 0*_iTempDimen );
- double* pdblImag1 = _pdblImag + ( 1*_iTempDimen );
- double* pdblImag2 = _pdblImag + ( 2*_iTempDimen );
- double* pdblImag3 = _pdblImag + ( 3*_iTempDimen );
- double* pdblImag4 = _pdblImag + ( 4*_iTempDimen );
- double* pdblImag5 = _pdblImag + ( 5*_iTempDimen );
- double* pdblImag6 = _pdblImag + ( 6*_iTempDimen );
- double* pdblImag7 = _pdblImag + ( 7*_iTempDimen );
-
- double arg,c1,s1,c2,s2,c3,s3,c4,s4,c5,s5,c6,s6,c7,s7 ;
- double dblAReal0, dblAReal1, dblAReal2, dblAReal3, dblAReal4, dblAReal5, dblAReal6, dblAReal7;
- double dblAImag0, dblAImag1, dblAImag2, dblAImag3, dblAImag4, dblAImag5, dblAImag6, dblAImag7 ;
-
- double dblBReal0, dblBReal1, dblBReal2, dblBReal3, dblBReal4, dblBReal5, dblBReal6, dblBReal7;
- double dblBImag0, dblBImag1, dblBImag2, dblBImag3, dblBImag4, dblBImag5, dblBImag6, dblBImag7 ;
-
- double dblTReal , dblTImag ;
-
- double dblPi2 = 8 * atan (1);
- double dblP7 = 1 / sqrt (2) ;
- double scale = dblPi2 /_iLengt ;
-
-
-
- for ( j = 0 ; j < _iTempDimen ; i++)
- {
- arg=j*scale;
-
- c1= cos(arg);
- s1= sin(arg);
-
- c2= c1*c1 - s1*s1 ;
- s2= c1*s1 + c1*s1;
-
- c3= c1*c2 - s1*s2;
- s3= c2*s1 + s2*c1;
-
- c4= c2*c2 - s2*s2 ;
- s4= c2*s2 + c2*s2;
-
- c5= c2*c3 - s2*s3;
- s5= c3*s2 + s3*c2;
-
- c6= c3*c3 - s3*s3 ;
- s6= c3*s3 + c3*s3;
-
- c7= c3*c4 - s3*s4;
- s7= c4*s3 + s4*c3;
-
- for ( k = j ; k < _iDimen ; k += _iLengt )
- {
- dblAReal0 = pdblReal0[k] + pdblReal4[k];
- dblAReal1 = pdblReal1[k] + pdblReal5[k];
- dblAReal2 = pdblReal2[k] + pdblReal6[k];
- dblAReal3 = pdblReal3[k] + pdblReal7[k];
- dblAReal4 = pdblReal0[k] - pdblReal4[k];
- dblAReal5 = pdblReal1[k] - pdblReal5[k];
- dblAReal6 = pdblReal2[k] - pdblReal6[k];
- dblAReal7 = pdblReal3[k] - pdblReal7[k];
-
- dblAImag0 = pdblImag0[k] + pdblImag4[k];
- dblAImag1 = pdblImag1[k] + pdblImag5[k];
- dblAImag2 = pdblImag2[k] + pdblImag6[k];
- dblAImag3 = pdblImag3[k] + pdblImag7[k];
- dblAImag4 = pdblImag0[k] - pdblImag4[k];
- dblAImag5 = pdblImag1[k] - pdblImag5[k];
- dblAImag6 = pdblImag2[k] - pdblImag6[k];
- dblAImag7 = pdblImag3[k] - pdblImag7[k];
-
- dblBReal0 = dblAReal0 + dblAReal2;
- dblBReal1 = dblAReal1 + dblAReal3;
- dblBReal2 = dblAReal0 - dblAReal2;
- dblBReal3 = dblAReal1 - dblAReal3;
- dblBReal4 = dblAReal4 - dblAImag6;
- dblBReal5 = dblAReal5 - dblAImag7;
- dblBReal6 = dblAReal4 + dblAImag6;
- dblBReal7 = dblAReal5 + dblAImag7;
-
- dblBImag0 = dblAImag0 + dblAImag2;
- dblBImag1 = dblAImag1 + dblAImag3;
- dblBImag2 = dblAImag0 - dblAImag2;
- dblBImag3 = dblAImag1 - dblAImag3;
- dblBImag4 = dblAImag4 + dblAReal6;
- dblBImag5 = dblAImag5 + dblAReal7;
- dblBImag6 = dblAImag4 - dblAReal6;
- dblBImag7 = dblAImag5 - dblAReal7;
-
- pdblReal0[k] = dblBReal0 + dblBReal1;
- pdblImag0[k] = dblBImag0 + dblBImag1;
-
- if ( j > 1 )
- {
- pdblReal1[k] = c4*( dblBReal0 - dblBReal1 ) - s4*( dblBImag0 - dblBImag1 );
- pdblImag1[k] = c4*( dblBImag0 - dblBImag1 ) + s4*( dblBReal0 - dblBReal1 );
- pdblReal2[k] = c2*( dblBReal2 - dblBImag3 ) - s2*( dblBImag2 + dblBReal3 );
- pdblImag2[k] = c2*( dblBImag2 + dblBReal3 ) + s2*( dblBReal2 - dblBImag3 );
- pdblReal3[k] = c6*( dblBReal2 + dblBImag3 ) - s6*( dblBImag2 - dblBReal3 );
- pdblImag3[k] = c6*( dblBImag2 - dblBReal3 ) + s6*( dblBReal2 + dblBImag3 );
-
- dblTReal = dblP7*( dblBReal5 - dblBImag5 );
- dblTImag = dblP7*( dblBReal5 + dblBImag5 );
-
- pdblReal4[k] = c1*( dblBReal4 + dblTReal ) - s1*( dblBImag4 + dblTImag );
- pdblImag4[k] = c1*( dblBImag4 + dblTImag ) + s1*( dblBReal4 + dblTReal );
- pdblReal5[k] = c5*( dblBReal4 - dblTReal ) - s5*( dblBImag4 - dblTImag );
- pdblImag5[k] = c5*( dblBImag4 - dblTImag ) + s5*( dblBReal4 - dblTReal );
-
- dblTReal = - dblP7* ( dblBReal7 + dblBImag7);
- dblTImag = dblP7* ( dblBReal7 - dblBImag7);
-
- pdblReal6[k] = c3*( dblBReal6 + dblTReal ) - s3*( dblBImag6 + dblTImag );
- pdblImag6[k] = c3*( dblBImag6 + dblTImag ) + s3*( dblBReal6 + dblTReal );
- pdblReal7[k] = c7*( dblBReal6 - dblTReal ) - s7*( dblBImag6 - dblTImag );
- pdblImag7[k] = c7*( dblBImag6 - dblTImag ) + s7*( dblBReal6 - dblTReal );
- }
- else
- {
- pdblReal1[k] = dblBReal0 - dblBReal1;
- pdblImag1[k] = dblBImag0 - dblBImag1;
- pdblReal2[k] = dblBReal2 - dblBImag3;
- pdblImag2[k] = dblBImag2 + dblBReal3;
- pdblReal3[k] = dblBReal2 + dblBImag3;
- pdblImag3[k] = dblBImag2 - dblBReal3;
-
- dblTReal = dblP7 * ( dblBReal5 - dblBImag5 );
- dblTImag = dblP7 * ( dblBReal5 + dblBImag5 );
-
- pdblReal4[k] = dblBReal4 + dblTReal;
- pdblImag4[k] = dblBImag4 + dblTImag;
- pdblReal5[k] = dblBReal4 - dblTReal;
- pdblImag5[k] = dblBImag4 - dblTImag;
-
- dblTReal = -dblP7 * ( dblBReal7 + dblBImag7 );
- dblTImag = dblP7 * ( dblBReal7 - dblBImag7 );
-
- pdblReal6[k] = dblBReal6 + dblTReal;
- pdblImag6[k] = dblBImag6 + dblTImag;
- pdblReal7[k] = dblBReal6 - dblTReal;
- pdblImag7[k] = dblBImag6 - dblTImag;
-
- }
-
- }
-
- }
+ int j,k,kk;
+ double dblP7 = 1 / sqrt (2) ;
+ double dblPi2 = 8 * atan (1);
+
+ double scale, arg;
+ double c1,c2,c3,c4,c5,c6,c7;
+ double s1,s2,s3,s4,s5,s6,s7;
+
+ doubleComplex Atemp0,Atemp1,Atemp2,Atemp3,Atemp4,Atemp5,Atemp6,Atemp7;
+ doubleComplex Btemp0,Btemp1,Btemp2,Btemp3,Btemp4,Btemp5,Btemp6,Btemp7;
+
+ doubleComplex temp ;
+
+/*
+ double *cr0,*cr1,*cr2,*cr3,*cr4,*cr5,*cr6,*cr7;
+ double *ci0,*ci1,*ci2,*ci3,*ci4,*ci5,*ci6,*ci7;
+
+ cr0 = &(cc0[0].re);
+ cr1 = &(cc1[0].re);
+ cr2 = &(cc2[0].re);
+ cr3 = &(cc3[0].re);
+ cr4 = &(cc4[0].re);
+ cr5 = &(cc5[0].re);
+ cr6 = &(cc6[0].re);
+ cr7 = &(cc7[0].re);
+
+ ci0 = &(cc0[0].im);
+ ci1 = &(cc1[0].im);
+ ci2 = &(cc2[0].im);
+ ci3 = &(cc3[0].im);
+ ci4 = &(cc4[0].im);
+ ci5 = &(cc5[0].im);
+ ci6 = &(cc6[0].im);
+ ci7 = &(cc7[0].im);
+*/
+
+ scale = dblPi2/lengt;
+
+ for(j=1;j<=nxtlt;j++)
+ {
+ arg = (j-1)*scale;
+ c1 = cos(arg);
+ s1 = sin(arg);
+ c2 = c1*c1 - s1*s1;
+ s2 = c1*s1 + c1*s1;
+ c3 = c1*c2 - s1*s2;
+ s3 = c2*s1 + s2*c1;
+ c4 = c2*c2 - s2*s2;
+ s4 = c2*s2 + c2*s2;
+ c5 = c2*c3 - s2*s3;
+ s5 = c3*s2 + s3*c2;
+ c6 = c3*c3 - s3*s3;
+ s6 = c3*s3 + c3*s3;
+ c7 = c3*c4 - s3*s4;
+ s7 = c4*s3 + s4*c3;
+
+ for(k=j;k<=nthpo;k+=lengt)
+ {
+ kk = (k-1)*2; /* index by twos; re & im alternate */
+
+ Atemp0 = zadds ( cc0[kk] , cc4[kk] ) ;
+ Atemp1 = zadds ( cc1[kk] , cc5[kk] ) ;
+ Atemp2 = zadds ( cc2[kk] , cc6[kk] ) ;
+ Atemp3 = zadds ( cc3[kk] , cc7[kk] ) ;
+
+/*
+ ar0 = cr0[kk] + cr4[kk];
+ ar1 = cr1[kk] + cr5[kk];
+ ar2 = cr2[kk] + cr6[kk];
+ ar3 = cr3[kk] + cr7[kk];
+
+
+ ai0 = ci0[kk] + ci4[kk];
+ ai1 = ci1[kk] + ci5[kk];
+ ai2 = ci2[kk] + ci6[kk];
+ ai3 = ci3[kk] + ci7[kk];
+*/
+
+ Atemp4 = zdiffs ( cc0[kk] , cc4[kk] ) ;
+ Atemp5 = zdiffs ( cc1[kk] , cc5[kk] ) ;
+ Atemp6 = zdiffs ( cc2[kk] , cc6[kk] ) ;
+ Atemp7 = zdiffs ( cc3[kk] , cc7[kk] ) ;
+/*
+ ar4 = cr0[kk] - cr4[kk];
+ ar5 = cr1[kk] - cr5[kk];
+ ar6 = cr2[kk] - cr6[kk];
+ ar7 = cr3[kk] - cr7[kk];
+
+ ai4 = ci0[kk] - ci4[kk];
+ ai5 = ci1[kk] - ci5[kk];
+ ai6 = ci2[kk] - ci6[kk];
+ ai7 = ci3[kk] - ci7[kk];
+*/
+
+
+
+ Btemp0 = zadds ( Atemp0 , Atemp2 ) ;
+ Btemp1 = zadds ( Atemp1 , Atemp3 ) ;
+ Btemp2 = zdiffs ( Atemp0 , Atemp2 ) ;
+ Btemp3 = zdiffs ( Atemp1 , Atemp3 ) ;
+
+/*
+ br0 = ar0 + ar2;
+ br1 = ar1 + ar3;
+ br2 = ar0 - ar2;
+ br3 = ar1 - ar3;
+
+ bi0 = ai0 + ai2;
+ bi1 = ai1 + ai3;
+ bi2 = ai0 - ai2;
+ bi3 = ai1 - ai3;
+*/
+
+ Btemp4 = DoubleComplex ( zreals ( Atemp4 ) - zimags( Atemp6 ) , zimags ( Atemp4 ) + zreals( Atemp6 ) );
+ Btemp5 = DoubleComplex ( zreals ( Atemp5 ) - zimags( Atemp7 ) , zimags ( Atemp5 ) + zreals( Atemp7 ) );
+ Btemp6 = DoubleComplex ( zreals ( Atemp4 ) + zimags( Atemp6 ) , zimags ( Atemp4 ) - zreals( Atemp6 ) );
+ Btemp7 = DoubleComplex ( zreals ( Atemp5 ) + zimags( Atemp7 ) , zimags ( Atemp5 ) - zreals( Atemp7 ) );
+/*
+ br4 = ar4 - ai6;
+ br5 = ar5 - ai7;
+ br6 = ar4 + ai6;
+ br7 = ar5 + ai7;
+
+ bi4 = ai4 + ar6;
+ bi5 = ai5 + ar7;
+ bi6 = ai4 - ar6;
+ bi7 = ai5 - ar7;
+*/
+ cc0[kk] = zadds ( Btemp0 , Btemp1 );
+/*
+ cr0[kk] = br0 + br1;
+ ci0[kk] = bi0 + bi1;
+*/
+
+
+
+
+ if(j>1)
+ {
+ printf ( "on sup j>1\n");
+ cc1[kk] = DoubleComplex ( c4 * (zreals(Btemp0) - zreals(Btemp1)) - s4 * (zimags(Btemp0) - zimags(Btemp1)),
+ c4 * (zimags(Btemp0) - zimags(Btemp1)) + s4 * (zreals(Btemp0) - zreals(Btemp1)));
+
+ cc2[kk] = DoubleComplex ( c2 * (zreals(Btemp2) - zimags(Btemp3)) - s2 * (zimags(Btemp2) + zreals(Btemp3)) ,
+ c2 * (zimags(Btemp2) + zreals(Btemp3)) + s2 * (zreals(Btemp2) - zimags(Btemp3)));
+
+ cc3[kk] = DoubleComplex ( c6 * (zreals(Btemp2) + zimags(Btemp3)) - s6 * (zimags(Btemp2) - zreals(Btemp3)) ,
+ c6 * (zimags(Btemp2) - zreals(Btemp3)) + s6 * (zreals(Btemp2) + zimags(Btemp3)));
+
+/*
+ cr1[kk] = c4*(br0-br1) - s4*(bi0-bi1);
+ ci1[kk] = c4*(bi0-bi1) + s4*(br0-br1);
+
+ cr2[kk] = c2*(br2-bi3) - s2*(bi2+br3);
+ ci2[kk] = c2*(bi2+br3) + s2*(br2-bi3);
+
+ cr3[kk] = c6*(br2+bi3) - s6*(bi2-br3);
+ ci3[kk] = c6*(bi2-br3) + s6*(br2+bi3);
+*/
+
+ temp = DoubleComplex ( dblP7*(zreals ( Btemp5 ) - zimags( Btemp5 )) ,
+ dblP7*(zreals ( Btemp5 ) + zimags( Btemp5 )) );
+/*
+ tr = dblP7*(br5-bi5);
+ ti = dblP7*(br5+bi5);
+*/
+
+ cc4[kk] = DoubleComplex ( c1 * (zreals (Btemp4) + zreals(temp)) - s1 * (zimags (Btemp4) + zimags(temp)) ,
+ c1 * (zimags (Btemp4) + zimags(temp)) + s1 * (zreals (Btemp4) + zreals(temp)));
+ cc5[kk] = DoubleComplex ( c5 * (zreals (Btemp4) - zreals(temp)) - s5 * (zimags (Btemp4) - zimags(temp)) ,
+ c5 * (zimags (Btemp4) - zimags(temp)) + s5 * (zreals (Btemp4) - zreals(temp)));
+
+/*
+ cr4[kk] = c1*(br4+tr) - s1*(bi4+ti);
+ ci4[kk] = c1*(bi4+ti) + s1*(br4+tr);
+ cr5[kk] = c5*(br4-tr) - s5*(bi4-ti);
+ ci5[kk] = c5*(bi4-ti) + s5*(br4-tr);
+*/
+ temp = DoubleComplex ( - dblP7*(zreals ( Btemp7 ) + zimags( Btemp7 )) ,
+ dblP7*(zreals ( Btemp7 ) - zimags( Btemp7 )) );
+/*
+ tr = -dblP7*(br7+bi7);
+ ti = dblP7*(br7-bi7);
+*/
+ cc6[kk] = DoubleComplex ( c3 * (zreals (Btemp6) + zreals(temp)) - s3 * (zimags (Btemp6) + zimags(temp)) ,
+ c3 * (zimags (Btemp6) + zimags(temp)) + s3 * (zreals (Btemp6) + zreals(temp)));
+ cc7[kk] = DoubleComplex ( c7 * (zreals (Btemp6) - zreals(temp)) - s7 * (zimags (Btemp6) - zimags(temp)) ,
+ c7 * (zimags (Btemp6) - zimags(temp)) + s7 * (zreals (Btemp6) - zreals(temp)));
+/*
+ cr6[kk] = c3*(br6+tr) - s3*(bi6+ti);
+ ci6[kk] = c3*(bi6+ti) + s3*(br6+tr);
+ cr7[kk] = c7*(br6-tr) - s7*(bi6-ti);
+ ci7[kk] = c7*(bi6-ti) + s7*(br6-tr);
+*/
+ }
+ else
+ {
+ printf ( "on oupa sup j>1\n");
+ cc1[kk] = zdiffs ( Btemp0 , Btemp1 );
+ /*
+ cr1[kk] = br0 - br1;
+ ci1[kk] = bi0 - bi1;
+ */
+
+ cc2[kk] = DoubleComplex ( zreals ( Btemp2 ) - zimags( Btemp3 ) ,
+ zimags ( Btemp2 ) + zreals( Btemp3 ) );
+
+
+ cc3[kk] = DoubleComplex ( zreals ( Btemp2 ) + zimags( Btemp3 ) ,
+ zimags ( Btemp2 ) - zreals( Btemp3 ) );
+
+ /*
+ cr2[kk] = br2 - bi3;
+ cr3[kk] = br2 + bi3;
+
+ ci2[kk] = bi2 + br3;
+ ci3[kk] = bi2 - br3;
+ */
+
+ temp = DoubleComplex ( dblP7*(zreals ( Btemp5 ) - zimags( Btemp5 )) ,
+ dblP7*(zreals ( Btemp5 ) + zimags( Btemp5 )) );
+ /*
+ tr = dblP7*(br5-bi5);
+ ti = dblP7*(br5+bi5);
+ */
+ cc4[kk] = zadds ( Btemp4 , temp );
+ cc5[kk] = zdiffs ( Btemp4 , temp );
+ /*
+ cr4[kk] = br4 + tr;
+ ci4[kk] = bi4 + ti;
+ cr5[kk] = br4 - tr;
+ ci5[kk] = bi4 - ti;
+ */
+ temp = DoubleComplex ( - dblP7*(zreals ( Btemp7 ) + zimags( Btemp7 )) ,
+ dblP7*(zreals ( Btemp7 ) - zimags( Btemp7 )) );
+ /*
+ tr = -dblP7*(br7+bi7);
+ ti = dblP7*(br7-bi7);
+ */
+ cc6[kk] = zadds ( Btemp6 , temp );
+ cc7[kk] = zdiffs ( Btemp6 , temp );
+ /*
+ cr6[kk] = br6 + tr;
+ ci6[kk] = bi6 + ti;
+ cr7[kk] = br6 - tr;
+ ci7[kk] = bi6 - ti;
+ */
+
+
+
+
+
+
+ }
+
+ }
+ }
+
+ printf ( "plop %d \t%e \t%e\n" , 0 , zreals(cc0[kk]) , zimags(cc0[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 1 , zreals(cc1[kk]) , zimags(cc1[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 2 , zreals(cc2[kk]) , zimags(cc2[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 3 , zreals(cc3[kk]) , zimags(cc3[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 4 , zreals(cc4[kk]) , zimags(cc4[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 5 , zreals(cc5[kk]) , zimags(cc5[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 6 , zreals(cc6[kk]) , zimags(cc6[kk]));
+ printf ( "plop %d \t%e \t%e\n" , 7 , zreals(cc7[kk]) , zimags(cc7[kk]));
}
diff --git a/src/signalProcessing/fft/testDoubleFft.c b/src/signalProcessing/fft/testDoubleFft.c
new file mode 100644
index 00000000..05f4c28c
--- /dev/null
+++ b/src/signalProcessing/fft/testDoubleFft.c
@@ -0,0 +1,279 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING, which
+ * you should have received as part of this distribution. The terms
+ * are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+#include <assert.h>
+#include <stdio.h>
+#include "fft.h"
+
+
+#define ROW 1
+#define COLS2 2
+#define COLS4 4
+#define COLS8 8
+#define COLS16 16
+
+#define ZREAL_IN2 { 0.00022113462910056 , 0.33032709173858166 }
+#define ZIMAG_IN2 { 0.66538110421970487 , 0.62839178834110498 }
+
+#define ZREAL_IN4 { 0.84974523587152362, 0.68573101982474327, 0.87821648130193353, 0.06837403681129217}
+#define ZIMAG_IN4 { 0.56084860628470778, 0.66235693730413914, 0.72635067673400044, 0.19851438421756029}
+
+#define ZREAL_IN8 { 0.54425731627270579, 0.23207478970289230, 0.23122371966019273, 0.21646326314657927,\
+ 0.88338878145441413, 0.65251349471509457, 0.30760907428339124, 0.93296162132173777 }
+#define ZIMAG_IN8 { 0.21460078610107303, 0.31264199689030647, 0.36163610080257058, 0.2922266637906432,\
+ 0.56642488157376647, 0.48264719732105732, 0.33217189135029912, 0.59350947011262178}
+
+#define ZREAL_IN16 {0.23122371966019273, 0.21646326314657927, 0.88338878145441413, 0.65251349471509457,\
+ 0.30760907428339124, 0.93296162132173777, 0.21460078610107303, 0.31264199689030647,\
+ 0.36163610080257058, 0.2922266637906432 , 0.56642488157376647, 0.48264719732105732,\
+ 0.33217189135029912, 0.59350947011262178, 0.50153415976092219, 0.43685875833034515}
+#define ZIMAG_IN16 {0.26931248093023896, 0.63257448654621840, 0.40519540151581168, 0.91847078315913677,\
+ 0.04373343335464597, 0.48185089323669672, 0.26395560009405017, 0.41481037065386772,\
+ 0.28064980218186975, 0.12800584640353918, 0.77831285959109664, 0.21190304495394230,\
+ 0.11213546665385365, 0.68568959552794695, 0.15312166837975383, 0.69708506017923355}
+
+
+#define ZREAL_RESULT2 { 0.33054822636768222,- 0.33010595710948110}
+#define ZIMAG_RESULT2 { 1.29377289256080985, 0.03698931587859988}
+
+#define ZREAL_RESULT4 { 2.48206677380949259, 0.43537130765616894, 0.97385666053742170, -0.49231379851698875}
+#define ZIMAG_RESULT4 { 2.14807060454040766,- 0.78285905346274376, 0.42632796149700880, 0.45185491256415844}
+
+
+#define ZREAL_RESULT8 { 4.00049206055700779,-0.43357241280891956, 0.79836636409163475,-0.91119240848798977,\
+ -0.06753427721560001,-0.18576209864995416, 0.97926024347543716, 0.17400105922003017}
+#define ZIMAG_RESULT8 { 3.15585898794233799, 0.62132445165622818, 0.35205427557229996, 0.28289917172258683,\
+ -0.20619166828691959,-1.17220193335521805,-0.17761892452836037,-1.13931807191437073 }
+
+#define ZREAL_RESULT16 { 7.31841186061501503, 0.57213963313411265,-0.54757095809921363,-0.48628670926159856,\
+ -1.24745626002550125,-0.60260425121772254,-0.09566750389725764, 1.12013387649474438,\
+ -0.52123307064175606,-0.4866536676629296 , 1.98659065302356819,-0.8626986211125984 ,\
+ -0.61915938556194305,-0.27813937201980266,-1.53103677171080510,-0.01918993749322817}
+#define ZIMAG_RESULT16 { 6.47680679336190224, 0.33111151130330035,-0.19343861330849654, 0.12474172265893407,\
+ -1.0452539175748825 , 1.29632487527975693, 1.87557979276701658,-1.82623636350346352,\
+ -1.86397336795926094,-1.03154071610913434,-0.48573205481665604, 0.44539904220706855,\
+ -0.74425477534532547,-0.54299368721281471, 0.37996440777257234, 1.11249504536330601}
+
+
+static void zfftmaTest2 (void )
+{
+ int i = 0 ;
+
+ double tRealIn [] = ZREAL_IN2;
+ double tImagIn [] = ZIMAG_IN2 ;
+
+
+
+ double tRealResult [] = ZREAL_RESULT2 ;
+ double tImagResult [] = ZIMAG_RESULT2 ;
+
+
+
+ doubleComplex* out = (doubleComplex*) malloc ( sizeof(doubleComplex) * (unsigned int) (ROW*COLS2));
+ doubleComplex* in = DoubleComplexMatrix ( tRealIn , tImagIn , ROW*COLS2 );
+ doubleComplex* Result = DoubleComplexMatrix ( tRealResult , tImagResult ,ROW*COLS2) ;
+
+
+ zfftma ( in , ROW , COLS2 , out ) ;
+
+ /* if we don't add that test assert failed if result = 0 'cause then we have |(out - 0)|/|out| = 1*/
+
+
+ for ( i = 0 ; i < (ROW*COLS2 ) ; i++ )
+ {
+ printf ( "\t\t %d out : %e\t %e\t * i result : %e\t %e\t * i assert : : %e\t %e\t * i \n" ,
+ i ,zreals(out[i]) , zimags(out[i]), zreals (Result[i]) , zimags (Result[i]),
+ fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) ,
+ fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])));
+
+ if ( zreals(out[i]) < 1e-14 && zreals (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-12 );
+
+
+ if ( zimags(out[i]) < 1e-14 && zimags (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-12 ) ;
+
+ }
+
+
+}
+
+static void zfftmaTest4 (void )
+{
+ int i = 0 ;
+
+ double tRealIn [] = ZREAL_IN4;
+ double tImagIn [] = ZIMAG_IN4 ;
+
+
+
+ double tRealResult [] = ZREAL_RESULT4 ;
+ double tImagResult [] = ZIMAG_RESULT4 ;
+
+
+
+ doubleComplex* out = (doubleComplex*) malloc ( sizeof(doubleComplex) * (unsigned int) (ROW*COLS4));
+ doubleComplex* in = DoubleComplexMatrix ( tRealIn , tImagIn , ROW*COLS4 );
+ doubleComplex* Result = DoubleComplexMatrix ( tRealResult , tImagResult ,ROW*COLS4) ;
+
+
+ zfftma ( in , ROW , COLS4 , out ) ;
+
+ /* if we don't add that test assert failed if result = 0 'cause then we have |(out - 0)|/|out| = 1*/
+ for ( i = 0 ; i < (ROW*COLS4 ) ; i++ )
+ {
+ printf ( "\t\t %d out : %e\t %e\t * i result : %e\t %e\t * i assert : : %e\t %e\t * i \n" ,
+ i ,zreals(out[i]) , zimags(out[i]), zreals (Result[i]) , zimags (Result[i]),
+ fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) ,
+ fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])));
+
+ if ( zreals(out[i]) < 1e-14 && zreals (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-12 );
+
+
+ if ( zimags(out[i]) < 1e-14 && zimags (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-12 ) ;
+
+ }
+
+
+}
+
+
+static void zfftmaTest8 (void )
+{
+ int i = 0 ;
+
+ double tRealIn [] = ZREAL_IN8;
+ double tImagIn [] = ZIMAG_IN8 ;
+
+
+
+ double tRealResult [] = ZREAL_RESULT8 ;
+ double tImagResult [] = ZIMAG_RESULT8 ;
+
+
+
+ doubleComplex* out = (doubleComplex*) malloc ( sizeof(doubleComplex) * (unsigned int) (ROW*COLS8));
+ doubleComplex* in = DoubleComplexMatrix ( tRealIn , tImagIn , ROW*COLS8 );
+ doubleComplex* Result = DoubleComplexMatrix ( tRealResult , tImagResult ,ROW*COLS8) ;
+
+
+ zfftma ( in , ROW , COLS8 , out ) ;
+
+ /* if we don't add that test assert failed if result = 0 'cause then we have |(out - 0)|/|out| = 1*/
+ for ( i = 0 ; i < (ROW*COLS8 ) ; i++ )
+ {
+ printf ( "\t\t %d out : %e\t %e\t * i result : %e\t %e\t * i assert : : %e\t %e\t * i \n" ,
+ i ,zreals(out[i]) , zimags(out[i]), zreals (Result[i]) , zimags (Result[i]),
+ fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) ,
+ fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])));
+
+ if ( zreals(out[i]) < 1e-14 && zreals (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-12 );
+
+
+ if ( zimags(out[i]) < 1e-14 && zimags (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-12 ) ;
+
+ }
+
+
+}
+
+static void zfftmaTest16 (void )
+{
+ int i = 0 ;
+
+ double tRealIn [] = ZREAL_IN16;
+ double tImagIn [] = ZIMAG_IN16 ;
+
+
+
+ double tRealResult [] = ZREAL_RESULT16 ;
+ double tImagResult [] = ZIMAG_RESULT16 ;
+
+
+
+ doubleComplex* out = (doubleComplex*) malloc ( sizeof(doubleComplex) * (unsigned int) (ROW*COLS16));
+ doubleComplex* in = DoubleComplexMatrix ( tRealIn , tImagIn , ROW*COLS16 );
+ doubleComplex* Result = DoubleComplexMatrix ( tRealResult , tImagResult ,ROW*COLS16) ;
+
+
+ zfftma ( in , ROW , COLS16 , out ) ;
+
+ /* if we don't add that test assert failed if result = 0 'cause then we have |(out - 0)|/|out| = 1*/
+ for ( i = 0 ; i < (ROW*COLS16 ) ; i++ )
+ {
+
+ printf ( "\t\t %d out : %e\t %e\t * i result : %e\t %e\t * i assert : : %e\t %e\t * i \n" ,
+ i ,zreals(out[i]) , zimags(out[i]), zreals (Result[i]) , zimags (Result[i]),
+ fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) ,
+ fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])));
+
+/*
+ if ( zreals(out[i]) < 1e-14 && zreals (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-12 );
+
+
+ if ( zimags(out[i]) < 1e-14 && zimags (Result[i]) < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-12 ) ;
+ */
+
+
+ }
+
+
+}
+
+static int testFft(void) {
+
+ printf("\n>>>> FFT Tests\n");
+ printf("\t>>>> Matrix Double Realt Tests\n");
+ /*dfftmaTest();*/
+
+ printf("\n\n\n");
+ printf("\t>>>> Vector 2 Double Complex Tests\n");
+ zfftmaTest2();
+ printf("\t>>>> Vector 4 Double Complex Tests\n");
+ zfftmaTest4();
+ printf("\t>>>> Vector 8 Double Complex Tests\n");
+ zfftmaTest8();
+ printf("\t>>>> Vector 16 Double Complex Tests\n");
+ zfftmaTest16();
+
+
+ return 0;
+}
+
+
+
+int main(void) {
+ assert(testFft() == 0);
+ return 0;
+}
diff --git a/src/signalProcessing/fft/testFloatFft.c b/src/signalProcessing/fft/testFloatFft.c
new file mode 100644
index 00000000..ee71e73d
--- /dev/null
+++ b/src/signalProcessing/fft/testFloatFft.c
@@ -0,0 +1,36 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING, which
+ * you should have received as part of this distribution. The terms
+ * are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+#include <assert.h>
+#include <stdio.h>
+#include "fft.h"
+
+
+
+
+
+
+static int testFft(void) {
+
+ printf("\n>>>> FFT Tests\n");
+ printf("\t>>>> Matrix Double Realt Tests\n");
+
+
+ return 0;
+}
+
+
+
+int main(void) {
+ assert(testFft() == 0);
+ return 0;
+}
diff --git a/src/signalProcessing/fft/zfftma.c b/src/signalProcessing/fft/zfftma.c
new file mode 100644
index 00000000..2b0a7510
--- /dev/null
+++ b/src/signalProcessing/fft/zfftma.c
@@ -0,0 +1,143 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2008 - INRIA - Allan SIMON
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING, which
+ * you should have received as part of this distribution. The terms
+ * are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+
+#include "fft.h"
+#include <stdio.h>
+
+void zfftma ( doubleComplex* in , int rows, int cols, doubleComplex* out)
+{
+
+int size = rows*cols ;
+int sizeTemp = 0;
+int sizeWorkSpace = 2* size;
+
+int rowsTemp = 0 ;
+int colsTemp = 0 ;
+
+int ierr = 0 ;
+int isn = -1;
+int i = 0;
+
+double* realIn = (double*) malloc ( sizeof (double) * (unsigned int) size );
+double* imagIn = (double*) malloc ( sizeof (double) * (unsigned int) size );
+
+
+doubleComplex* inTemp = (doubleComplex*) malloc ( sizeof (doubleComplex) * (unsigned int) size );
+
+double* workSpace = (double*) malloc ( sizeof (double) * (unsigned int) sizeWorkSpace );
+
+
+
+zimaga ( in , size , imagIn) ;
+zreala ( in , size , realIn) ;
+
+if ( rows == 1 || cols == 1 )
+{
+ printf ( "it'a vector \n" ) ;
+
+ sizeTemp = (int) pow ( 2 , (int ) (log( size + 0.5 ) /log ( 2 ))) ;
+ printf ("pow %e , temp %d \n" , pow ( 2 , (int )(log( size +0.5 ) /log ( 2 ))), sizeTemp);
+
+ if ( size == sizeTemp )
+ {
+ if ( size <= pow ( 2 , 15 ))
+ {
+ printf ( "we call fft842 \n" ) ;
+ fft842 ( in , size , 0 );
+ }
+ else
+ {
+ printf ( "we call dfft2 \n" ) ;
+ dfft2 ( realIn , imagIn , 1 , size , 1 , isn , ierr , workSpace , sizeWorkSpace );
+ }
+
+
+ }
+ else
+ {
+ printf ( "we call dfft2 2\n" ) ;
+ dfft2 ( realIn , imagIn , 1 , size , 1 , isn , ierr , workSpace , sizeWorkSpace );
+ }
+
+}
+else
+ {
+ printf ( "it'a matrix \n" ) ;
+ rowsTemp = (int) pow ( 2 , log( rows + 0.5) /log ( 2 )) ;
+ colsTemp = (int) pow ( 2 , log( cols + 0.5) /log ( 2 )) ;
+
+ if ( cols == colsTemp)
+ {
+ if ( cols <= pow ( 2 , 15 ))
+ {
+ for ( i = 0 ; i < rows ; i++ )
+ {
+ fft842 ( &in[ cols*i] , cols , 0);
+ }
+ }
+ else
+ {
+ dfft2 ( realIn, imagIn ,rows , cols , 1 , isn , ierr ,workSpace , sizeWorkSpace );
+ }
+ }
+ else
+ {
+ dfft2 ( realIn, imagIn ,rows , cols , 1 , isn , ierr ,workSpace , sizeWorkSpace );
+ }
+
+ /*second call*/
+
+ if ( 2*rows <= sizeWorkSpace )
+ {
+ if ( rowsTemp == rows )
+ {
+ if ( rows <= pow ( 2 ,15) )
+ {
+ /*compute the fft on each line of the matrix */
+ for (i = 0 ; i < cols ; i++ )
+ {
+ C2F(zcopy) ( rows, in + i, cols, inTemp , 1 );
+
+ fft842( inTemp , rows , 0);
+
+ C2F(zcopy) ( rows, inTemp , cols, in + i, 1 );
+
+ }
+ }
+ else
+ {
+ dfft2 ( realIn, imagIn, 1, rows, cols, isn, ierr, workSpace, sizeWorkSpace);
+ }
+ }
+ else
+ {
+ dfft2 ( realIn, imagIn, 1, rows, cols, isn, ierr, workSpace, sizeWorkSpace);
+ }
+ }
+ else
+ {
+ dfft2 ( realIn, imagIn, 1, rows, cols, isn, ierr, workSpace, sizeWorkSpace);
+ }
+ }
+
+
+
+
+
+for ( i = 0 ; i < size ; i++)
+{
+out[i] = DoubleComplex ( zreals(in[i]) , zimags(in[i]) );
+}
+
+
+}
generated by cgit (git 2.34.1) at 2025-03-26 15:27:21 +0000