/* * 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 "ifft_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 ir8tx ( 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 ); } } } }