summaryrefslogtreecommitdiff
path: root/src/signalProcessing/fft
diff options
context:
space:
mode:
authorsimon2008-09-26 11:39:49 +0000
committersimon2008-09-26 11:39:49 +0000
commit13381564d0570ccc23df6d042a613f4f025edc67 (patch)
treec93d527ff447cc84ccbed8fd4336be5ce24a8018 /src/signalProcessing/fft
parent9efc2ee1a7e5373e0e5504eb2cf23ebdd3eb2233 (diff)
downloadscilab2c-13381564d0570ccc23df6d042a613f4f025edc67.tar.gz
scilab2c-13381564d0570ccc23df6d042a613f4f025edc67.tar.bz2
scilab2c-13381564d0570ccc23df6d042a613f4f025edc67.zip
corrected all tests for fft842 everything ok for that function,now time to see the general function
Diffstat (limited to 'src/signalProcessing/fft')
-rw-r--r--src/signalProcessing/fft/fft842.c6
-rw-r--r--src/signalProcessing/fft/r2tx.c6
-rw-r--r--src/signalProcessing/fft/r4tx.c6
-rw-r--r--src/signalProcessing/fft/r8tx.c160
-rw-r--r--src/signalProcessing/fft/testDoubleFft.c95
5 files changed, 121 insertions, 152 deletions
diff --git a/src/signalProcessing/fft/fft842.c b/src/signalProcessing/fft/fft842.c
index 5aef5752..12aa8002 100644
--- a/src/signalProcessing/fft/fft842.c
+++ b/src/signalProcessing/fft/fft842.c
@@ -78,7 +78,7 @@ void fft842 (doubleComplex* b, int size , int in)
{
nxtlt = 0x1 << (n2pow - 3*ipass);
lengt = 8*nxtlt;
- printf ( "on appelle r%dtx \n" , 8);
+
r8tx(nxtlt,nthpo,lengt,
@@ -93,7 +93,6 @@ void fft842 (doubleComplex* b, int size , int in)
if(n2pow%3 == 1)
{
/* radix 2 iteration needed */
- printf ( "on appelle r%dtx \n" , 2);
r2tx(nthpo,b,b+1);
@@ -103,7 +102,7 @@ void fft842 (doubleComplex* b, int size , int in)
if(n2pow%3 == 2)
{
/* radix 4 iteration needed */
- printf ( "on appelle r%dtx \n" , 4);
+
r4tx(nthpo,b,b+1,b+2,b+3);
}
@@ -135,6 +134,7 @@ void fft842 (doubleComplex* b, int size , int in)
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;
diff --git a/src/signalProcessing/fft/r2tx.c b/src/signalProcessing/fft/r2tx.c
index 8c945061..74039da5 100644
--- a/src/signalProcessing/fft/r2tx.c
+++ b/src/signalProcessing/fft/r2tx.c
@@ -17,7 +17,7 @@
*/
void r2tx(int nthpo, doubleComplex* c0, doubleComplex* c1)
{
- int k,kk;
+ int kk;
/* double *cr0, *ci0, *cr1, *ci1, r1, fi1;*/
doubleComplex temp ;
@@ -26,9 +26,9 @@ void r2tx(int nthpo, doubleComplex* c0, doubleComplex* c1)
cr1 = &(c1[0].re);
ci1 = &(c1[0].im);*/
- for(k=0;k<nthpo;k+=2)
+ for(kk=0;kk<nthpo;kk+=2)
{
- kk = k;
+
temp = zadds ( c0[kk] , c1[kk] );
c1[kk] = zdiffs( c0[kk] , c1[kk] );
diff --git a/src/signalProcessing/fft/r4tx.c b/src/signalProcessing/fft/r4tx.c
index e923d338..f7b6400d 100644
--- a/src/signalProcessing/fft/r4tx.c
+++ b/src/signalProcessing/fft/r4tx.c
@@ -19,12 +19,12 @@
*/
void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, doubleComplex* c3)
{
- int k,kk;
+ int kk;
doubleComplex temp1 , temp2 , temp3 , temp4 ;
- for(k=1;k<=nthpo;k+=4)
+ for(kk=0;kk<nthpo;kk+=4)
{
- kk = k -1; /* real and imag parts alternate */
+ /* real and imag parts alternate */
temp1 = zadds ( c0[kk] , c2[kk] ) ;
diff --git a/src/signalProcessing/fft/r8tx.c b/src/signalProcessing/fft/r8tx.c
index 042b39a8..8c44488e 100644
--- a/src/signalProcessing/fft/r8tx.c
+++ b/src/signalProcessing/fft/r8tx.c
@@ -27,7 +27,7 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
{
- int j,k,kk;
+ int j , kk;
double dblP7 = 1 / sqrt (2) ;
double dblPi2 = 8 * atan (1);
@@ -40,36 +40,14 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
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++)
+ for(j=0;j<nxtlt;j++)
{
- arg = (j-1)*scale;
+ arg = j*scale;
c1 = cos(arg);
s1 = sin(arg);
c2 = c1*c1 - s1*s1;
@@ -85,9 +63,9 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
c7 = c3*c4 - s3*s4;
s7 = c4*s3 + s4*c3;
- for(k=j;k<=nthpo;k+=lengt)
+ for(kk=j;kk<nthpo;kk+=lengt)
{
- kk = k - 1;/* (k-1)*2*/ /* index by twos; re & im alternate */
+ ;/* (k-1)*2*/ /* index by twos; re & im alternate */
Atemp0 = zadds ( cc0[kk] , cc4[kk] ) ;
@@ -95,78 +73,34 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
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)
+
+ if(j>0)
{
cc1[kk] = DoubleComplex ( (c4 * (zreals(Btemp0) - zreals(Btemp1))) - (s4 * (zimags(Btemp0) - zimags(Btemp1))),
c4 * (zimags(Btemp0) - zimags(Btemp1)) + s4 * (zreals(Btemp0) - zreals(Btemp1)));
@@ -177,61 +111,31 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
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
{
cc1[kk] = zdiffs ( Btemp0 , Btemp1 );
- /*
- cr1[kk] = br0 - br1;
- ci1[kk] = bi0 - bi1;
- */
cc2[kk] = DoubleComplex ( zreals ( Btemp2 ) - zimags( Btemp3 ) ,
zimags ( Btemp2 ) + zreals( Btemp3 ) );
@@ -240,46 +144,18 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
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;
- */
-
-
-
-
}
diff --git a/src/signalProcessing/fft/testDoubleFft.c b/src/signalProcessing/fft/testDoubleFft.c
index a55c3f20..374071ae 100644
--- a/src/signalProcessing/fft/testDoubleFft.c
+++ b/src/signalProcessing/fft/testDoubleFft.c
@@ -20,6 +20,7 @@
#define COLS4 4
#define COLS8 8
#define COLS16 16
+#define COLS32 32
#define ZREAL_IN2 { 0.00022113462910056 , 0.33032709173858166 }
#define ZIMAG_IN2 { 0.66538110421970487 , 0.62839178834110498 }
@@ -43,6 +44,24 @@
0.11213546665385365, 0.68568959552794695, 0.15312166837975383, 0.69708506017923355}
+#define ZREAL_IN32 {0.21132486546412110,0.75604385416954756,0.00022113462910056,0.33032709173858166,\
+ 0.66538110421970487,0.62839178834110498,0.84974523587152362,0.68573101982474327,\
+ 0.87821648130193353,0.06837403681129217,0.56084860628470778,0.66235693730413914,\
+ 0.72635067673400044,0.19851438421756029,0.54425731627270579,0.23207478970289230,\
+ 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_IN32 {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,\
+ 0.84155184263363481,0.40620247554033995,0.40948254754766822,0.87841258011758327,\
+ 0.11383596854284406,0.19983377400785685,0.56186607433483005,0.58961773291230202,\
+ 0.68539796629920602,0.89062247332185507,0.50422128057107329,0.34936154074966908,\
+ 0.38737787725403905,0.92228986788541079,0.94881842611357570,0.34353372454643250}
+
#define ZREAL_RESULT2 { 0.33054822636768222,- 0.33010595710948110}
#define ZIMAG_RESULT2 { 1.29377289256080985, 0.03698931587859988}
@@ -64,6 +83,25 @@
-1.86397336795926094,-1.03154071610913434,-0.48573205481665604, 0.44539904220706855,\
-0.74425477534532547,-0.54299368721281471, 0.37996440777257234, 1.11249504536330601}
+#define ZREAL_RESULT32 { 15.3165711835026741,-1.79021577127059173,-1.66659611407065089, 0.17525916470909797,\
+ -1.16958628014871602, 0.58684741669397522, 0.03947542161511042, 0.99740008842981942,\
+ -0.46323241293430328, 2.122539701124051 ,-1.52963914564883940, 0.87990417229605744,\
+ 0.58569127383151542,-0.18198535589432135, 0.26043384746900655,-1.11204765363415392,\
+ 0.35329844802618027,-1.47568616310628631,-2.03487116744967844,-3.19495610958970166,\
+ 0.81026376203844086,-0.46366666776372734,-1.63150209835186510, 0.77334707088593369,\
+ -0.35098156332969666,-1.63498270669406387, 0.67411467120679691, 2.80538085483913147,\
+ -1.62281507315555107,-0.5600265995962992 ,-0.48984739061140237, 1.75450689143393301}
+
+
+
+#define ZIMAG_RESULT32 { 15.509232945740223 ,-0.47962381296807621, 0.21213951866464975, 0.88442937061831350,\
+ 1.52924554070524898,-0.6313403060045536 ,-2.25908603874729419,-1.3836292677373856 ,\
+ -1.18231281638145447,-2.14767090006699668, 1.5452016553381984 ,-3.15355126536920993,\
+ 2.10199273301496747, 0.67530605269461363,-2.03603600735261558, 1.2309547869577584 ,\
+ -1.9912955537438393 ,-1.9668221895811833 , 0.29575245179739662, 1.34815224953105273,\
+ 1.508921339902356 ,-0.40084285801706099, 2.96716476331614754, 1.08125713762201059,\
+ -1.39964522421360016,-1.30777696073860294,-1.13169784714423916,-2.00872755010475013,\
+ 0.53915777133569487,-2.45178696294021004 , 1.56509394479014063, 1.5558426888499468}
static void zfftmaTest2 (void )
{
@@ -262,6 +300,60 @@ static void zfftmaTest16 (void )
}
+static void zfftmaTest32 (void )
+{
+ int i = 0 ;
+
+ double tRealIn [] = ZREAL_IN32;
+ double tImagIn [] = ZIMAG_IN32 ;
+
+
+
+ double tRealResult [] = ZREAL_RESULT32 ;
+ double tImagResult [] = ZIMAG_RESULT32 ;
+
+
+
+ doubleComplex* out = (doubleComplex*) malloc ( sizeof(doubleComplex) * (unsigned int) (ROW*COLS32));
+ doubleComplex* in = DoubleComplexMatrix ( tRealIn , tImagIn , ROW*COLS32 );
+
+
+
+ zfftma ( in , ROW , COLS32 , 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*COLS32 ) ; 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]),
+ tRealResult[i] ,
+ tImagResult[i],
+ fabs( zreals(out[i]) - tRealResult[i] ) / fabs (zreals (out[i])) ,
+ fabs( zimags(out[i]) - tImagResult[i] ) / fabs (zimags (out[i])));
+
+
+ if ( zreals(out[i]) < 1e-14 && tRealResult[i] < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zreals(out[i]) - tRealResult[i] ) / fabs (zreals (out[i])) < 1e-12 );
+
+
+ if ( zimags(out[i]) < 1e-14 && tImagResult[i] < 1e-18 )
+ assert ( 1 ) ;
+ else
+ assert ( fabs( zimags(out[i]) - tImagResult[i] ) / fabs (zimags (out[i])) < 1e-12 ) ;
+
+ }
+
+}
+
static int testFft(void) {
printf("\n>>>> FFT Tests\n");
@@ -278,7 +370,8 @@ static int testFft(void) {
printf("\t>>>> Vector 8 Double Complex Tests\n");
zfftmaTest8();
-
+ printf("\t>>>> Vector 32 Double Complex Tests\n");
+ zfftmaTest32();
return 0;
}