diff options
Diffstat (limited to 'src/signalProcessing/fft/r8tx.c')
| -rw-r--r-- | src/signalProcessing/fft/r8tx.c | 438 |
1 files changed, 278 insertions, 160 deletions
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])); } |