From 7bcc74849a84b87464ca6cc48cdbde8206a3a602 Mon Sep 17 00:00:00 2001 From: torset Date: Fri, 19 Dec 2008 15:26:21 +0000 Subject: debug fft and add new tests --- src/signalProcessing/fft/Makefile.am | 9 +- src/signalProcessing/fft/Makefile.in | 37 ++- src/signalProcessing/fft/dfftbi.c | 4 +- src/signalProcessing/fft/dfftmx.c | 517 +++++++++++++--------------------- src/signalProcessing/fft/testMatFft.c | 219 ++++++++++++++ src/signalProcessing/fft/zfftma.c | 57 ++-- 6 files changed, 483 insertions(+), 360 deletions(-) create mode 100644 src/signalProcessing/fft/testMatFft.c (limited to 'src') diff --git a/src/signalProcessing/fft/Makefile.am b/src/signalProcessing/fft/Makefile.am index c4d5faf4..46d75349 100644 --- a/src/signalProcessing/fft/Makefile.am +++ b/src/signalProcessing/fft/Makefile.am @@ -58,9 +58,9 @@ check_LDADD = $(top_builddir)/type/libDoubleComplex.la \ $(top_builddir)/operations/subtraction/libSubtraction.la \ @LIBMATH@ -check_PROGRAMS = testFloatFft testDoubleFft +check_PROGRAMS = testFloatFft testDoubleFft testMatFft -TESTS = testFloatFft testDoubleFft +TESTS = testFloatFft testDoubleFft testMatFft # # -*- Fftine Tests -*- @@ -72,3 +72,8 @@ testFloatFft_LDADD = $(check_LDADD) testDoubleFft_SOURCES = testDoubleFft.c testDoubleFft_CFLAGS = $(check_INCLUDES) testDoubleFft_LDADD = $(check_LDADD) + + +testMatFft_SOURCES = testMatFft.c +testMatFft_CFLAGS = $(check_INCLUDES) +testMatFft_LDADD = $(check_LDADD) diff --git a/src/signalProcessing/fft/Makefile.in b/src/signalProcessing/fft/Makefile.in index bd2a10dd..4be18b9d 100644 --- a/src/signalProcessing/fft/Makefile.in +++ b/src/signalProcessing/fft/Makefile.in @@ -32,8 +32,10 @@ PRE_UNINSTALL = : POST_UNINSTALL = : build_triplet = @build@ host_triplet = @host@ -check_PROGRAMS = testFloatFft$(EXEEXT) testDoubleFft$(EXEEXT) -TESTS = testFloatFft$(EXEEXT) testDoubleFft$(EXEEXT) +check_PROGRAMS = testFloatFft$(EXEEXT) testDoubleFft$(EXEEXT) \ + testMatFft$(EXEEXT) +TESTS = testFloatFft$(EXEEXT) testDoubleFft$(EXEEXT) \ + testMatFft$(EXEEXT) subdir = signalProcessing/fft DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 @@ -82,6 +84,12 @@ testFloatFft_DEPENDENCIES = $(am__DEPENDENCIES_1) testFloatFft_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \ $(LIBTOOLFLAGS) --mode=link $(CCLD) $(testFloatFft_CFLAGS) \ $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@ +am_testMatFft_OBJECTS = testMatFft-testMatFft.$(OBJEXT) +testMatFft_OBJECTS = $(am_testMatFft_OBJECTS) +testMatFft_DEPENDENCIES = $(am__DEPENDENCIES_1) +testMatFft_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \ + $(LIBTOOLFLAGS) --mode=link $(CCLD) $(testMatFft_CFLAGS) \ + $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@ DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)/includes depcomp = $(SHELL) $(top_srcdir)/config/depcomp am__depfiles_maybe = depfiles @@ -95,9 +103,9 @@ LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \ --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \ $(LDFLAGS) -o $@ SOURCES = $(libFft_la_SOURCES) $(testDoubleFft_SOURCES) \ - $(testFloatFft_SOURCES) + $(testFloatFft_SOURCES) $(testMatFft_SOURCES) DIST_SOURCES = $(libFft_la_SOURCES) $(testDoubleFft_SOURCES) \ - $(testFloatFft_SOURCES) + $(testFloatFft_SOURCES) $(testMatFft_SOURCES) ETAGS = etags CTAGS = ctags DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) @@ -262,6 +270,9 @@ testFloatFft_LDADD = $(check_LDADD) testDoubleFft_SOURCES = testDoubleFft.c testDoubleFft_CFLAGS = $(check_INCLUDES) testDoubleFft_LDADD = $(check_LDADD) +testMatFft_SOURCES = testMatFft.c +testMatFft_CFLAGS = $(check_INCLUDES) +testMatFft_LDADD = $(check_LDADD) all: all-am .SUFFIXES: @@ -337,6 +348,9 @@ testDoubleFft$(EXEEXT): $(testDoubleFft_OBJECTS) $(testDoubleFft_DEPENDENCIES) testFloatFft$(EXEEXT): $(testFloatFft_OBJECTS) $(testFloatFft_DEPENDENCIES) @rm -f testFloatFft$(EXEEXT) $(testFloatFft_LINK) $(testFloatFft_OBJECTS) $(testFloatFft_LDADD) $(LIBS) +testMatFft$(EXEEXT): $(testMatFft_OBJECTS) $(testMatFft_DEPENDENCIES) + @rm -f testMatFft$(EXEEXT) + $(testMatFft_LINK) $(testMatFft_OBJECTS) $(testMatFft_LDADD) $(LIBS) mostlyclean-compile: -rm -f *.$(OBJEXT) @@ -355,6 +369,7 @@ distclean-compile: @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@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatFft-testMatFft.Po@am__quote@ .c.o: @am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $< @@ -468,6 +483,20 @@ testFloatFft-testFloatFft.obj: testFloatFft.c @AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testFloatFft_CFLAGS) $(CFLAGS) -c -o testFloatFft-testFloatFft.obj `if test -f 'testFloatFft.c'; then $(CYGPATH_W) 'testFloatFft.c'; else $(CYGPATH_W) '$(srcdir)/testFloatFft.c'; fi` +testMatFft-testMatFft.o: testMatFft.c +@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatFft_CFLAGS) $(CFLAGS) -MT testMatFft-testMatFft.o -MD -MP -MF $(DEPDIR)/testMatFft-testMatFft.Tpo -c -o testMatFft-testMatFft.o `test -f 'testMatFft.c' || echo '$(srcdir)/'`testMatFft.c +@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/testMatFft-testMatFft.Tpo $(DEPDIR)/testMatFft-testMatFft.Po +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='testMatFft.c' object='testMatFft-testMatFft.o' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatFft_CFLAGS) $(CFLAGS) -c -o testMatFft-testMatFft.o `test -f 'testMatFft.c' || echo '$(srcdir)/'`testMatFft.c + +testMatFft-testMatFft.obj: testMatFft.c +@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatFft_CFLAGS) $(CFLAGS) -MT testMatFft-testMatFft.obj -MD -MP -MF $(DEPDIR)/testMatFft-testMatFft.Tpo -c -o testMatFft-testMatFft.obj `if test -f 'testMatFft.c'; then $(CYGPATH_W) 'testMatFft.c'; else $(CYGPATH_W) '$(srcdir)/testMatFft.c'; fi` +@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/testMatFft-testMatFft.Tpo $(DEPDIR)/testMatFft-testMatFft.Po +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='testMatFft.c' object='testMatFft-testMatFft.obj' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatFft_CFLAGS) $(CFLAGS) -c -o testMatFft-testMatFft.obj `if test -f 'testMatFft.c'; then $(CYGPATH_W) 'testMatFft.c'; else $(CYGPATH_W) '$(srcdir)/testMatFft.c'; fi` + mostlyclean-libtool: -rm -f *.lo diff --git a/src/signalProcessing/fft/dfftbi.c b/src/signalProcessing/fft/dfftbi.c index e98b3bdd..ad0fb0f1 100644 --- a/src/signalProcessing/fft/dfftbi.c +++ b/src/signalProcessing/fft/dfftbi.c @@ -89,7 +89,7 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , int maxf ; int itype; int istkgt ; - int isize[] = {1,1,1,2,2} ; + int nf = abs ( n ) ; @@ -263,7 +263,7 @@ void dfftbi ( double* a , double* b , int nseg , int n , int nspn , /*istkgt = ( lnow*isize[1] -1)/isize[1] + 2;*/ istkgt = lnow + 1 ; /*i = ( (istkgt - 1 + nitems) * isize[1] -1) / isize[1] + 3 ;*/ - i = ( ( lnow + nitems) * isize[1] -1) / isize[1] + 3 ; + i = lnow + nitems + 2 ; istak = (int*) realloc ( istak ,sizeof (int) * (unsigned int) i); rstak = (double*) malloc ( sizeof (double) * (unsigned int) i); diff --git a/src/signalProcessing/fft/dfftmx.c b/src/signalProcessing/fft/dfftmx.c index 60129b1c..74cfcb4a 100644 --- a/src/signalProcessing/fft/dfftmx.c +++ b/src/signalProcessing/fft/dfftmx.c @@ -97,7 +97,6 @@ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, int _iIsn, int _iM, int _iKt, double* _pdblWt, double* _pdblCk, double* _pdblBt, double* _pdblSk, int* _piNp, int* _piNfac) { - int retVal = 0 ; a = _pdblA ; b = _pdblB ; @@ -119,41 +118,23 @@ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, nt = inc*ntot ; ks = inc*nspan; rad = atan ( 1 ); - c72 = cos (rad/0.6250); s72 = sin (rad/0.6250); s120= sqrt(0.750); - - - - preliminaryWork () ; - while ( retVal == 0 ) - { - - retVal = factorTransform ( ) ; - } - + while ( retVal == 0 ) retVal = factorTransform ( ) ; np[0] = ks ; - if ( kt != 0) - { - permute_stage1 ( ) ; - } - - if ( 2*kt + 1 < m ) - { - - permute_stage2 ( ) ; - } - + if ( kt != 0) permute_stage1 ( ) ; + + if ( 2*kt + 1 < m ) permute_stage2 ( ) ; _pdblA = a ; _pdblB = b ; @@ -170,14 +151,9 @@ Sous-Fonctions /* this function only set the value of variable */ void preliminaryWork (void) { - - int lim ; - - - - s72 = -s72 ; - s120= -s120; - rad = -rad ; + s72 = -s72 ; + s120= -s120; + rad = -rad ; kspan = ks ; nn = nt -inc ; @@ -185,19 +161,12 @@ void preliminaryWork (void) /* sin , cos values are re-initialized each lim steps */ - lim = 32 ; - klim = lim*jc ; + klim = 32*jc ; i = 0; jf = 0 ; - maxf = m -kt ; - maxf = nfac[maxf-1] ; - - - - if ( kt > 0 ) - maxf = max ( nfac[kt-1] , maxf ); - + if ( kt > 0 ) maxf = max ( nfac[kt-1] , nfac[m-kt-1] ); + else maxf = nfac[m-kt-1] ; } @@ -214,8 +183,6 @@ int factorTransform (void) int retVal = 42; - - dr = 8 * (double)jc/(double)kspan ; cd = 2 * sin(0.5*dr*rad)*sin(0.5*dr*rad); sd = sin(dr*rad) ; @@ -227,38 +194,32 @@ int factorTransform (void) switch ( nfac[i-1] ) { - case 2 : + case 2 :printf("passge en 2\n"); /*transform for factor of 2 (including rotation factor)*/ - retVal = pre_fOf2Trans( ) ; - if ( retVal == 0 ) - { - - factorOf2Transform ( ) ; + retVal = pre_fOf2Trans() ; + if ( retVal == 0 ) factorOf2Transform () ; - } break ; - case 4 : + case 4 :printf("passge en 4\n"); /*transform for factor of 4 */ kspnn = kspan ; kspan = kspan >> 2 ; /*kspan /= 4 */ - - retVal = factorOf4Transform ( ) ; + retVal = factorOf4Transform () ; break ; - case 3 : - + case 3 :printf("passge en 3\n"); +printf("valeurs init : %d and %d and %d\n",nfac[i-1],kspan,nn); k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; - - + factorOf3Transform ( ) ; break ; - case 5 : + case 5 :printf("passge en 5\n"); k = nfac[i-1] ; kspnn = kspan ; @@ -267,19 +228,14 @@ switch ( nfac[i-1] ) factorOf5Transform ( ) ; break ; - default : + default :printf("passge en autre\n"); k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; - - if ( nfac[i-1] != jf) - { - - preFOtherTransform ( ) ; - } - + if ( nfac[i-1] != jf) preFOtherTransform ( ) ; + factorOfOtherTransform ( ) ; break ; } @@ -288,24 +244,12 @@ switch ( nfac[i-1] ) if ( retVal == 42 ) { - - - if ( i != m) - { - - retVal = mulByRotationFactor ( ) ; - } - else - retVal = 1 ; + if ( i != m) retVal = mulByRotationFactor ( ) ; + else retVal = 1 ; } - - if ( retVal == 1 ) - return 1 ; /*goto permute */ - else - return 0 ; /*goto factor_transform => once again*/ - - + if ( retVal == 1 ) return 1 ; /*goto permute */ + else return 0 ; /*goto factor_transform => once again*/ } @@ -315,59 +259,45 @@ void permute_stage1 (void) int retVal = 1 ; - pre_sqFactor2NormlOrder ( ) ; + pre_sqFactor2NormlOrder () ; if ( n == ntot ) /*permutation for single-variate transform (optional code)*/ while ( retVal == 1) { - single_sqFactor2NormlOrder ( ) ; + single_sqFactor2NormlOrder () ; retVal = post_sqFactor2NormlOrder () ; } else /*permutation for multivariate transform*/ - while ( retVal == 1) - { - retVal = multi_sqFactor2NormlOrder ( ); - } - + while ( retVal == 1) retVal = multi_sqFactor2NormlOrder (); } void permute_stage2 (void) { - - - kspnn = np[kt] ; - - - /*permutation for square-free facotrs of n */ - - nonSqFactor2NormOrder ( ) ; - - /*determine the permutation cycles of length greater than 1*/ - - detPermutCycles ( ); - - - - j = k3 + 1; - nt -= kspnn ; - i = nt - inc + 1 ; - while ( nt >= 0 ) - { - - reorderMatrix ( ) ; - - j = k3 + 1 ; - nt -= kspnn ; - i = nt - inc + 1 ; + kspnn = np[kt] ; + + /*permutation for square-free facotrs of n */ + nonSqFactor2NormOrder () ; + + /*determine the permutation cycles of length greater than 1*/ + detPermutCycles (); + + j = k3 + 1; + nt -= kspnn ; + i = nt - inc + 1 ; + while ( nt >= 0 ) + { + reorderMatrix ( ) ; + + j = k3 + 1 ; + nt -= kspnn ; + i = nt - inc + 1 ; } - - } -/** ************************************** +/***************************************** Sous-Sous-Fonctions ******************************************/ @@ -377,42 +307,30 @@ Sous-Sous-Fonctions int pre_fOf2Trans (void) { - int ktemp = 0 ; + kspan /= 2; + k1 = kspan + 2 ; + /*50*/ + do{ + do{ + k2 = kk + kspan ; + ak = a[k2-1] ; + bk = b[k2-1] ; - kspan /= 2; - k1 = kspan + 2 ; -/*50*/ - do - { + a[k2-1] = a[kk-1] - ak; + b[k2-1] = b[kk-1] - bk; - k2 = kk + kspan ; - ak = a[k2-1] ; - bk = b[k2-1] ; + a[kk-1] = a[kk-1] + ak; + b[kk-1] = b[kk-1] + bk; - a[k2-1] = a[kk-1] - ak; - b[k2-1] = b[kk-1] - bk; + kk = k2 + kspan ; + }while (kk <= nn); - a[kk-1] = a[kk-1] + ak; - b[kk-1] = b[kk-1] + bk; - - kk = k2 + kspan ; - ktemp = kk ; - - if ( kk > nn ) - { - kk -= nn ; - } - }while (ktemp <= nn ||( kk <= jc && ktemp <= nn)); - - - - if ( kk > kspan ) - return 1 ; /*goto350*/ - else - - return 0 ;/*goto60*/ + kk -= nn ; + }while (kk <= jc); + if ( kk > kspan ) return 1 ; /*goto350*/ + else return 0 ; /*goto60*/ } @@ -421,79 +339,66 @@ int pre_fOf2Trans (void) int factorOf2Transform (void) { -int ktemp = 0 ; - - - - -do /*60*/ - { - c1 = 1 - cd ; - s1 = sd ; - mm = min( k1/2 , klim); - - do/* do 80 */ - { - do - { - k2 = kk + kspan; - - ak = a[kk-1] - a[k2-1]; - bk = b[kk-1] - b[k2-1]; - - a[kk-1] = a[kk-1] + a[k2-1]; - b[kk-1] = b[kk-1] + b[k2-1]; - - a[k2-1] = c1*ak - s1*bk; - b[k2-1] = s1*ak + c1*bk; - - - kk = k2 + kspan; - ktemp = kk ; - if (kk >= nt) - { - k2 = kk - nt; - c1 = -c1; - kk = k1 - k2; - - } - - }while ( ktemp < nt || (kk > k2 && ( ktemp >= nt)) ); - - kk += jc; - - if ( kk <= mm ) /* 70 */ - { - ak = c1 - ( cd*c1+sd*s1) ; - s1 += (sd*c1-cd*s1) ; - /*c the following three statements compensate for truncation - c error. if rounded arithmetic is used, substitute - c c1=ak*/ - c1 = 0.5/(ak*ak+s1*s1) + 0.5 ; - s1 *= c1 ; - c1 *= ak ; - } - else - { - if ( kk < k2 ) /*90*/ - { - - s1 = dr*rad*((double)(kk-1)/(double)jc); - c1 = cos(s1) ; - s1 = sin(s1) ; - mm = min(k1/2,mm+klim); - } - } - - }while ( kk <= mm || ( kk > mm && kk < k2 )); - - k1 += (inc+inc) ; - kk = (k1-kspan)/2 + jc; - - } while ( kk <= jc*2 ); - - - return 0 ; /*goto40*/ + do /*60*/ {/*while ( kk <= jc*2 )*/ + c1 = 1 - cd ; + s1 = sd ; + mm = min( k1/2 , klim); + + do/* do 80 */ {/*while ( kk <= mm || ( kk > mm && kk < k2 ))*/ + do {/*while(kk > k2) */ + do { /*while ( kk < nt )*/ + k2 = kk + kspan; + + ak = a[kk-1] - a[k2-1]; + bk = b[kk-1] - b[k2-1]; + + a[kk-1] = a[kk-1] + a[k2-1]; + b[kk-1] = b[kk-1] + b[k2-1]; + + a[k2-1] = c1*ak - s1*bk; + b[k2-1] = s1*ak + c1*bk; + + kk = k2 + kspan; + }while ( kk < nt ); + + k2 = kk - nt; + c1 = -c1; + kk = k1 - k2; + + + }while (kk > k2); + + kk += jc; + + if ( kk <= mm ) /* 70 */ + { + ak = c1 - ( cd*c1+sd*s1) ; + s1 += (sd*c1-cd*s1) ; + /*c the following three statements compensate for truncation + c error. if rounded arithmetic is used, substitute + c c1=ak*/ + c1 = 0.5/(ak*ak+s1*s1) + 0.5 ; + s1 *= c1 ; + c1 *= ak ; + } + else { + if ( kk < k2 ) /*90*/ { + s1 = dr*rad*((double)(kk-1)/(double)jc); + c1 = cos(s1) ; + s1 = sin(s1) ; + mm = min(k1/2,mm+klim); + } + } + + } while ( kk <= mm || ( kk > mm && kk < k2 )); + + k1 += (inc+inc) ; + kk = (k1-kspan)/2 + jc; + + } while ( kk <= jc*2 ); + + + return 0 ; /*goto40*/ } @@ -501,7 +406,6 @@ do /*60*/ int factorOf4Transform (void) { - int return_value = 0 ; /*120*/ @@ -537,8 +441,6 @@ int factorOf4Transform (void) void f4t_150 (void) { - int sign = 1 ; - do{ k1 = kk + kspan ; k2 = k1 + kspan ; @@ -562,17 +464,11 @@ void f4t_150 (void) b[kk-1] = bkp + bjp ; bjp = bkp - bjp ; - if ( isn < 0 ) - sign = 1 ; - else - sign = -1 ; - + akp = akm + bjm ; + akm = akm - bjm ; - akp = akm +(sign * bjm ); - akm = akm -(sign * bjm ); - - bkp = bkm +(sign * ajm) ; - bkm = bkm -(sign * ajm) ; + bkp = bkm - ajm ; + bkm = bkm + ajm ; if ( s1 == 0 )/*190*/ { @@ -597,14 +493,14 @@ void f4t_150 (void) a[k2-1] = bjp*c2 + ajp*s2 ; a[k3-1] = bkm*c3 + akm*s3 ; } - }while ( kk < nt ) ; + kk=k3+kspan; + }while ( kk <= nt ) ; } int f4t_170 (void) { - kk += ( jc - nt ) ; if ( kk <= mm ) @@ -620,12 +516,13 @@ int f4t_170 (void) c1 = 0.5/(c2*c2+s1*s1) + 0.5 ; s1 *= c1 ; - c2 *= c1 ; + c1 *= c2 ; /*140*/ c2 = c1*c1 - s1*s1 ; s2 = c1*s1*2 ; + c3 = c2*c1 - s2*s1 ; s3 = c2*s1 + s2*c1 ; @@ -637,7 +534,7 @@ int f4t_170 (void) if ( kk <= kspan ) { s1 = dr*rad * (kk-1)/jc ; - c2 = cos (s1) ; + c1 = cos (s1) ; s1 = sin (s1) ; mm = min ( kspan , mm + klim ); @@ -645,6 +542,7 @@ int f4t_170 (void) c2 = c1*c1 - s1*s1 ; s2 = c1*s1*2 ; + c3 = c2*c1 - s2*s1 ; s3 = c2*s1 + s2*c1 ; return 0 ; @@ -659,114 +557,99 @@ int f4t_170 (void) void factorOf3Transform (void) { -int ktemp = 0 ; -do - { - - k1 = kk + kspan ; - k2 = k1 + kspan ; + do{ + do{ + k1 = kk + kspan ; + k2 = k1 + kspan ; - ak = a[kk-1] ; - bk = b[kk-1] ; + ak = a[kk-1] ; + bk = b[kk-1] ; - aj = a[k1-1] + a[k2-1] ; - bj = b[k1-1] + b[k2-1] ; + aj = a[k1-1] + a[k2-1] ; + bj = b[k1-1] + b[k2-1] ; - a[kk-1] = ak + aj ; - b[kk-1] = bk + bj ; + a[kk-1] = ak + aj ; + b[kk-1] = bk + bj ; - ak = -0.5*aj + ak ; - bk = -0.5*bj + bk ; - - aj = (a[k1-1] - a[k2-1])*s120 ; - bj = (b[k1-1] - b[k2-1])*s120 ; - - a[k1-1] = ak - bj ; - b[k1-1] = bk + aj ; - a[k2-1] = ak + bj ; - b[k2-1] = bk - aj ; - - kk = k2 + kspan ; - ktemp = kk ; + ak = -0.5*aj + ak ; + bk = -0.5*bj + bk ; + aj = (a[k1-1] - a[k2-1])*s120 ; + bj = (b[k1-1] - b[k2-1])*s120 ; + a[k1-1] = ak - bj ; + b[k1-1] = bk + aj ; + a[k2-1] = ak + bj ; + b[k2-1] = bk - aj ; - if ( kk >= nn ) - kk -= nn ; - - }while ( ktemp < nn || (kk <= kspan && ( ktemp >= nn)) ); + kk = k2 + kspan ; + } while (kk < nn); + + kk -= nn ; + }while (kk <= kspan); } void factorOf5Transform (void) { - int ktemp ; - - c2 = c72*c72 - s72 *s72 ; - s2 = 2 * c72*s72; - - do - { - k1 = kk + kspan ; - k2 = k1 + kspan ; - k3 = k2 + kspan ; - k4 = k3 + kspan ; + c2 = c72*c72 - s72 *s72 ; + s2 = 2 * c72*s72; + do{ + do{ + k1 = kk + kspan ; + k2 = k1 + kspan ; + k3 = k2 + kspan ; + k4 = k3 + kspan ; - akp = a[k1-1] + a[k4-1] ; - akm = a[k1-1] - a[k4-1] ; - bkp = b[k1-1] + b[k4-1] ; - bkm = b[k1-1] - b[k4-1] ; + akp = a[k1-1] + a[k4-1] ; + akm = a[k1-1] - a[k4-1] ; - ajp = a[k2-1] + a[k3-1] ; - ajm = a[k2-1] - a[k3-1] ; + bkp = b[k1-1] + b[k4-1] ; + bkm = b[k1-1] - b[k4-1] ; - bjp = b[k2-1] + b[k3-1] ; - bjm = b[k2-1] - b[k3-1] ; + ajp = a[k2-1] + a[k3-1] ; + ajm = a[k2-1] - a[k3-1] ; - aa = a[kk-1] ; - bb = b[kk-1] ; + bjp = b[k2-1] + b[k3-1] ; + bjm = b[k2-1] - b[k3-1] ; - a[kk-1] = aa + akp + ajp; - b[kk-1] = bb + bkp + bjp; + aa = a[kk-1] ; + bb = b[kk-1] ; - ak = akp*c72 + ajp*c2 + aa ; - bk = bkp*c72 + bjp*c2 + bb ; + a[kk-1] = aa + akp + ajp; + b[kk-1] = bb + bkp + bjp; - aj = akm*s72 + ajm*s2 ; - bj = bkm*s72 + bjm*s2 ; + ak = akp*c72 + ajp*c2 + aa ; + bk = bkp*c72 + bjp*c2 + bb ; - a[k1-1] = ak - bj ; - a[k4-1] = ak + bj ; - b[k1-1] = bk + aj ; - b[k4-1] = bk - aj ; - - ak = akp*c2 + ajp*c72 + aa ; - bk = bkp*c2 + bjp*c72 + bb ; - - aj = akm*s2 - ajm*s72 ; + aj = akm*s72 + ajm*s2 ; + bj = bkm*s72 + bjm*s2 ; - bj = bkm*s2 - bjm*s72 ; + a[k1-1] = ak - bj ; + a[k4-1] = ak + bj ; + b[k1-1] = bk + aj ; + b[k4-1] = bk - aj ; - a[k2-1] = ak - bj ; - a[k3-1] = ak + bj ; - b[k2-1] = bk + aj ; - b[k3-1] = bk - aj ; + ak = akp*c2 + ajp*c72 + aa ; + bk = bkp*c2 + bjp*c72 + bb ; - kk = k4 + kspan; - ktemp = kk ; + aj = akm*s2 - ajm*s72 ; + bj = bkm*s2 - bjm*s72 ; - if ( kk >= nn ) - kk -= nn ; - - - }while (ktemp < nn || ( kk <= kspan && ktemp >= nn)); - + a[k2-1] = ak - bj ; + a[k3-1] = ak + bj ; + b[k2-1] = bk + aj ; + b[k3-1] = bk - aj ; + kk = k4 + kspan; + }while (kk < nn); + kk -= nn ; + }while (kk <= kspan); } /* this function is the general case of non factor of 2 factor , the factorof3transform and factorof5trandform are just @@ -775,7 +658,7 @@ special case of this one */ void preFOtherTransform (void) { - +printf("0.k=%d \n",k); jf = k ; s1 = (rad*8)/k ; diff --git a/src/signalProcessing/fft/testMatFft.c b/src/signalProcessing/fft/testMatFft.c new file mode 100644 index 00000000..728074b5 --- /dev/null +++ b/src/signalProcessing/fft/testMatFft.c @@ -0,0 +1,219 @@ + +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008-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 +#include +#include +#include +#include "fft.h" + +#define test1 {0.2113249000000000099586,0.3303270999999999846253,0.8497451999999999783242,0.0683740000000000042180,\ + 0.7560438999999999909463,0.6653810999999999475918,0.6857309999999999794440,0.5608486000000000304411,\ + 0.0002211000000000000075,0.6283917999999999448590,0.8782164999999999555058,0.6623569000000000261963} + +#define test2 {0.2113249000000000099586,0.6857309999999999794440,\ + 0.3303270999999999846253,0.5608486000000000304411,\ + 0.8497451999999999783242,0.0002211000000000000075,\ + 0.0683740000000000042180,0.6283917999999999448590,\ + 0.7560438999999999909463,0.8782164999999999555058,\ + 0.6653810999999999475918,0.6623569000000000261963} + +#define test3 {0.2113249000000000099586,0.7560438999999999909463,0.0002211000000000000075,\ + 0.3303270999999999846253,0.6653810999999999475918,0.6283917999999999448590,\ + 0.8497451999999999783242,0.6857309999999999794440,0.8782164999999999555058,\ + 0.0683740000000000042180,0.5608486000000000304411,0.6623569000000000261963} + +#define test4 {0.2113249000000000099586,0.0683740000000000042180,0.6857309999999999794440,0.6283917999999999448590,\ + 0.3303270999999999846253,0.7560438999999999909463,0.5608486000000000304411,0.8782164999999999555058,\ + 0.8497451999999999783242,0.6653810999999999475918,0.0002211000000000000075,0.6623569000000000261963} + +#define test6 {0.2113249000000000099586,0.8497451999999999783242,0.7560438999999999909463,0.6857309999999999794440,0.0002211000000000000075,0.8782164999999999555058,\ + 0.3303270999999999846253,0.0683740000000000042180,0.6653810999999999475918,0.5608486000000000304411,0.6283917999999999448590,0.6623569000000000261963} + +#define test9 {1,2,3,4,5,6,7,8,9} + +#define RRESULT1 {6.2969621000000000066166,-0.1908142591446193858662,0.0676336000000001824617,\ +-1.4461027999999997994962,-0.7577484999999999359943,0.2136987591446193046618,\ +0.4656031000000000474870,0.2136987591446193046618,-0.7577484999999999359943,\ +-1.4461027999999997994962,0.0676336000000001824617,-0.1908142591446193858662} + +#define IRESULT1 {0,-0.0438615735192828815059,-1.0572926567728160307524,0.3325204999999997745874,\ +0.3011445725036896448046,-1.3036713264807171519521,0,1.3036713264807171519521,\ +-0.3011445725036896448046,-0.3325204999999997745874,1.0572926567728160307524,\ +0.0438615735192828815059} + +#define RRESULT2 {6.2969621000000000066166,-0.5345696999999998144659,\ +0.0676336000000000714394,-0.3918127000000001247670,\ +-0.7577485000000003800835,-1.2843509999999995763176,\ +0.4656030999999999364647,1.0404605000000000103455,\ +-0.7577485000000003800835,-1.2843509999999995763176,\ +0.0676336000000000714394,-0.3918127000000001247670} + +#define IRESULT2 {0,0,\ +-1.05729265677281625280,0.6392575178288871784460,\ +0.3011445725036897003157,-1.043770536118125757952,\ +0,0,\ +-0.3011445725036897003157,1.043770536118125757952,\ +1.05729265677281625280,-0.6392575178288871784460} + +#define RRESULT3 {6.2969621000000000066166,-0.9588242499999999335358,-0.9588242499999999335358,\ +-1.4461027999999997994962,-0.3545214899631756733101,-0.1146366100368243767793,\ +0.4656031000000000474870,0.7607519500000001766438,0.7607519500000001766438,\ +-1.4461027999999997994962,-0.1146366100368243767793,-0.3545214899631756733101} + +#define IRESULT3 {0,0.4319893196725672623870,-0.4319893196725672623870,\ +0.3325204999999997745874,1.0479284784196345459151,-0.5945896784196345041806,\ +0,0.5437395057261036335206,-0.5437395057261036335206,\ +-0.3325204999999997745874,0.5945896784196345041806,-1.0479284784196345459151} + +#define RRESULT4 {6.2969621000000000066166,0.1445965000000000166835,-1.0205663000000000373291,\ + 0.1445965000000000166835,\ +-0.7577485000000003800835,-0.6754837907274804376101,0.8107183000000001138474,\ + -0.8923310092725194042629,\ +-0.7577485000000003800835,-0.8923310092725194042629,0.8107183000000001138474,\ + -0.6754837907274804376101} + +#define IRESULT4 {0,-0.6791661999999999421718,0,0.6791661999999999421718,\ +0.3011445725036897003157,-1.4357905268456061520510,-0.2297677979542606663266,\ + -0.4349033268456063971996,\ +-0.3011445725036897003157,0.4349033268456063971996,0.2297677979542606663266,\ + 1.4357905268456061520510} + +#define RRESULT6 {6.2969621000000000066166,-0.5006002499999999688640,-0.4661336499999999549715,\ + -1.113582300000000024909,-0.4661336499999999549715,-0.5006002499999999688640,\ +0.46560309999999915931,0.5234847499999999431708,-0.2239812499999997985611,\ + -1.7786232999999995740836,-0.2239812499999997985611,0.5234847499999999431708} + +#define IRESULT6 {0,0.1475342691353691004785,-1.2256565689200087376776,0,1.2256565689200087376776,\ + -0.1475342691353691004785,\ +0,1.112275483826065114457,-0.1327806603564972709464,0,0.1327806603564972709464,\ + -1.112275483826065114457} + +#define RRESULT9 {45,-4.5,-4.5,-13.5,0,0,-13.5,0,0} + +#define IRESULT9 {0,2.598076211353316011810,- 2.598076211353316011810,7.7942286340599471472501,0,0,- 7.7942286340599471472501,0,0} + +static int testFft(void){ + int i; + + double inR1[]=test1; + double inR2[]=test2; + double inR3[]=test3; + double inR4[]=test4; + double inR6[]=test6; + double inR9[]=test9; + + double resR1[]=RRESULT1; + double resI1[]=IRESULT1; + double resR2[]=RRESULT2; + double resI2[]=IRESULT2; + double resR3[]=RRESULT3; + double resI3[]=IRESULT3; + double resR4[]=RRESULT4; + double resI4[]=IRESULT4; + double resR6[]=RRESULT6; + double resI6[]=IRESULT6; + double resR9[]=RRESULT9; + double resI9[]=IRESULT9; + + doubleComplex *in1, *in2, *in3, *in4, *in6, *in9, *out1, *out2, *out3, *out4, *out6, *out9; + + in1=malloc((uint)12*sizeof(doubleComplex)); + in2=malloc((uint)12*sizeof(doubleComplex)); + in3=malloc((uint)12*sizeof(doubleComplex)); + in4=malloc((uint)12*sizeof(doubleComplex)); + in6=malloc((uint)12*sizeof(doubleComplex)); + in9=malloc((uint)9*sizeof(doubleComplex)); + out1=malloc((uint)12*sizeof(doubleComplex)); + out2=malloc((uint)12*sizeof(doubleComplex)); + out3=malloc((uint)12*sizeof(doubleComplex)); + out4=malloc((uint)12*sizeof(doubleComplex)); + out6=malloc((uint)12*sizeof(doubleComplex)); + out9=malloc((uint)9*sizeof(doubleComplex)); + + for (i=0;i<12;i++){ + in1[i]=DoubleComplex(inR1[i],0); + in2[i]=DoubleComplex(inR2[i],0); + in3[i]=DoubleComplex(inR3[i],0); + in4[i]=DoubleComplex(inR4[i],0); + in6[i]=DoubleComplex(inR6[i],0); + } + for (i=0;i<9;i++){ + in9[i]=DoubleComplex(inR9[i],0); + + } + + + zfftma(in1, 1, 12, out1); + zfftma(in2, 2, 6, out2); + zfftma(in3, 3, 4, out3); + zfftma(in4, 4, 3, out4); + zfftma(in6, 6, 2, out6); + zfftma(in9, 3, 3, out9); + + /* !!!!!!!!!!!!!!!!!!!!!!! + for the imaginary part, the assert is out + res instead of out - res + cause I export the transposate of the result matrix and the transposate change the sign + of the imaginary part. + And instead of change all the define, I only change the sign of the assert.*/ + for (i=0;i<12;i++){ + if (zreals(out1[i])>1e-16) assert( (fabs(zreals(out1[i])-resR1[i]) / fabs(zreals(out1[i]))) < 1e-14 ); + else assert(1); + if (zimags(out1[i])>1e-16) assert( (fabs(zimags(out1[i])+resI1[i]) / fabs(zimags(out1[i]))) < 1e-14 ); + else assert(1); + } + + for (i=0;i<12;i++){ + if (zreals(out2[i])>1e-16) assert( (fabs(zreals(out2[i])-resR2[i]) / fabs(zreals(out2[i]))) < 1e-14 ); + else assert(1); + if (zimags(out2[i])>1e-16) assert( (fabs(zimags(out2[i])+resI2[i]) / fabs(zimags(out2[i]))) < 1e-14 ); + else assert(1); + } + + for (i=0;i<12;i++){ + if (zreals(out3[i])>1e-16) assert( (fabs(zreals(out3[i])-resR3[i]) / fabs(zreals(out3[i]))) < 1e-14 ); + else assert(1); + if (zimags(out3[i])>1e-16) assert( (fabs(zimags(out3[i])+resI3[i]) / fabs(zimags(out3[i]))) < 1e-14 ); + else assert(1); + } + + for (i=0;i<12;i++){ + if (zreals(out4[i])>1e-16) assert( (fabs(zreals(out4[i])-resR4[i]) / fabs(zreals(out4[i]))) < 1e-14 ); + else assert(1); + if (zimags(out4[i])>1e-16) assert( (fabs(zimags(out4[i])+resI4[i]) / fabs(zimags(out4[i]))) < 1e-14 ); + else assert(1); + } + + for (i=0;i<12;i++){ + if (zreals(out6[i])>1e-16) assert( (fabs(zreals(out6[i])-resR6[i]) / fabs(zreals(out6[i]))) < 1e-14 ); + else assert(1); + if (zimags(out6[i])>1e-16) assert( (fabs(zimags(out6[i])+resI6[i]) / fabs(zimags(out6[i]))) < 1e-14 ); + else assert(1); + } + + for (i=0;i<9;i++){ + if (zreals(out9[i])>1e-16) assert( (fabs(zreals(out9[i])-resR9[i]) / fabs(zreals(out9[i]))) < 1e-16 ); + else assert(1); + + if (zimags(out9[i])>1e-15) assert( (fabs(zimags(out9[i])-resI9[i]) / fabs(zimags(out9[i]))) < 1e-15 ); + else assert(1); + } + return 0; +} + + +int main(void) { + printf(">>> Fft Matrices Double Tests <<<"); + assert(testFft() == 0); + return 0; +} + diff --git a/src/signalProcessing/fft/zfftma.c b/src/signalProcessing/fft/zfftma.c index a2d552e1..66573455 100644 --- a/src/signalProcessing/fft/zfftma.c +++ b/src/signalProcessing/fft/zfftma.c @@ -11,7 +11,7 @@ */ #define FFT842 1 #define DFFT2 0 - +#include #include #include #include "fft.h" @@ -25,9 +25,7 @@ void zfftma ( doubleComplex* in , int rows, int cols, doubleComplex* out) int size = rows*cols ; int sizeTemp = 0; - /* - int sizeWorkSpace = size *2; - */ + int rowsTemp = 0 ; int colsTemp = 0 ; @@ -42,11 +40,6 @@ void zfftma ( doubleComplex* in , int rows, int cols, doubleComplex* out) 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) ; @@ -62,69 +55,63 @@ void zfftma ( doubleComplex* in , int rows, int cols, doubleComplex* out) } else { - dfft2 ( realIn , imagIn , 1 , size , 1 , isn , ierr /*, workSpace , sizeWorkSpace*/ ); + dfft2 ( realIn , imagIn , 1 , size , 1 , isn , ierr ); } } else { - dfft2 ( realIn , imagIn , 1 , size , 1 , isn , ierr /*, workSpace , sizeWorkSpace */); + dfft2 ( realIn , imagIn , 1 , size , 1 , isn , ierr ); } } else { - rowsTemp = (int) pow ( 2 , log( rows + 0.5) /log ( 2 )) ; - colsTemp = (int) pow ( 2 , log( cols + 0.5) /log ( 2 )) ; - if ( cols == colsTemp) + rowsTemp = (int) pow ( 2 ,(int ) log( rows + 0.5) /log ( 2 )) ; + colsTemp = (int) pow ( 2 ,(int ) log( cols + 0.5) /log ( 2 )) ; + if ( rows == rowsTemp) { - if ( cols <= pow ( 2 , 15 )) + if ( rows <= pow ( 2 , 15 )) { - for ( i = 0 ; i < rows ; i++ ) + for ( i = 0 ; i < cols ; i++ ) { - fft842 ( &in[ cols*i] , cols , 0); + fft842 ( &in[ rows*i] , rows , 0); choosenAlgo = FFT842 ; } } else { - dfft2 ( realIn, imagIn ,rows , cols , 1 , isn , ierr/* ,workSpace , sizeWorkSpace */); + dfft2 ( realIn, imagIn ,cols , rows , 1 , isn , ierr); } } else { - dfft2 ( realIn, imagIn ,rows , cols , 1 , isn , ierr/* ,workSpace , sizeWorkSpace*/ ); + dfft2 ( realIn, imagIn ,cols , rows , 1 , isn , ierr); } /*second call*/ - if ( 2*rows <= 0 /* sizeWorkSpace*/ ) - { - if ( rowsTemp == rows ) + if ( colsTemp == cols ) { - if ( rows <= pow ( 2 ,15) ) + if ( cols <= 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 ); + for (i = 0 ; i < rows ; i++ ) + { + C2F(zcopy) ( cols, in + i, rows, inTemp , 1 ); - fft842( inTemp , rows , 0); + fft842( inTemp , cols , 0); choosenAlgo = FFT842 ; - C2F(zcopy) ( rows, inTemp , cols, in + i, 1 ); + C2F(zcopy) ( cols, inTemp , 1, in + i, rows ); } } else { - dfft2 ( realIn, imagIn, 1, rows, cols, isn, ierr/*, workSpace, sizeWorkSpace*/); + dfft2 ( realIn, imagIn, 1, cols, rows, isn, ierr); } } else { - dfft2 ( realIn, imagIn, 1, rows, cols, isn, ierr/*, workSpace, sizeWorkSpace*/); + dfft2 ( realIn, imagIn, 1, cols, rows, isn, ierr); } - } - else - { - dfft2 ( realIn, imagIn, 1, rows, cols, isn, ierr/*, workSpace, sizeWorkSpace*/); - } + } -- cgit