diff options
Diffstat (limited to '2.3-1/src/c/signalProcessing/fft/r8tx.c')
| -rw-r--r-- | 2.3-1/src/c/signalProcessing/fft/r8tx.c | 169 |
1 files changed, 169 insertions, 0 deletions
diff --git a/2.3-1/src/c/signalProcessing/fft/r8tx.c b/2.3-1/src/c/signalProcessing/fft/r8tx.c new file mode 100644 index 00000000..44866ce6 --- /dev/null +++ b/2.3-1/src/c/signalProcessing/fft/r8tx.c @@ -0,0 +1,169 @@ + +/* + * 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_internal.h" +#include <math.h> + + + +/* +** radix 8 iteration subroutine +*/ + +/* this function do in one turn the same computation that do radix 2 in three turns */ + +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 j , kk; + double dblP7 = 1 / sqrt(2.0) ; + double dblPi2 = 8 * atan(1.0); + + 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 ; + + + scale = dblPi2/lengt; + + + for(j=0;j<nxtlt;j++) + { + + 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(kk=j;kk<nthpo;kk+=lengt) + { + /* (k-1)*2*/ /* index by twos; re & im alternate */ + + + /* first turn the same as calling radix 2 with the input vector */ + /* but radix2 will have do it in three turn , radix8 do it in one */ + Atemp0 = zadds ( cc0[kk] , cc4[kk] ) ; + Atemp1 = zadds ( cc1[kk] , cc5[kk] ) ; + Atemp2 = zadds ( cc2[kk] , cc6[kk] ) ; + Atemp3 = zadds ( cc3[kk] , cc7[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] ) ; + + /* second turn the same as calling radix 2 with the vector transformed by a previous call of radix2 */ + /* the same here , three turns in one */ + Btemp0 = zadds ( Atemp0 , Atemp2 ) ; + Btemp1 = zadds ( Atemp1 , Atemp3 ) ; + Btemp2 = zdiffs ( Atemp0 , Atemp2 ) ; + Btemp3 = zdiffs ( Atemp1 , Atemp3 ) ; + + 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 ) ); + + /*third turn the same as calling radix 2 with the vector transformed by two previous call of radix2 */ + cc0[kk] = zadds ( Btemp0 , Btemp1 ); + + + + /* if we are not in the first turn */ + + 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))); + + 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))); + + + temp = DoubleComplex ( dblP7*(zreals ( Btemp5 ) - zimags( Btemp5 )) , + dblP7*(zreals ( Btemp5 ) + zimags( Btemp5 )) ); + + 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))); + + + temp = DoubleComplex ( - dblP7*(zreals ( Btemp7 ) + zimags( Btemp7 )) , + dblP7*(zreals ( Btemp7 ) - zimags( Btemp7 )) ); + + 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))); + + + + + } + else + { + cc1[kk] = zdiffs ( Btemp0 , Btemp1 ); + + cc2[kk] = DoubleComplex ( zreals ( Btemp2 ) - zimags( Btemp3 ) , + zimags ( Btemp2 ) + zreals( Btemp3 ) ); + + + cc3[kk] = DoubleComplex ( zreals ( Btemp2 ) + zimags( Btemp3 ) , + zimags ( Btemp2 ) - zreals( Btemp3 ) ); + + + temp = DoubleComplex ( dblP7*(zreals ( Btemp5 ) - zimags( Btemp5 )) , + dblP7*(zreals ( Btemp5 ) + zimags( Btemp5 )) ); + + cc4[kk] = zadds ( Btemp4 , temp ); + cc5[kk] = zdiffs ( Btemp4 , temp ); + + temp = DoubleComplex ( - dblP7*(zreals ( Btemp7 ) + zimags( Btemp7 )) , + dblP7*(zreals ( Btemp7 ) - zimags( Btemp7 )) ); + + cc6[kk] = zadds ( Btemp6 , temp ); + cc7[kk] = zdiffs ( Btemp6 , temp ); + + + } + + } + } + + +} |