00001 /* 00002 * THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE 00003 * CHECKED IT OUT USING THE COMMAND CHECKOUT. 00004 * 00005 * $Id: fft99_8c-source.html 2161 2006-05-19 16:55:03Z paulf $ 00006 * 00007 * Revision history: 00008 * $Log$ 00008 * Revision 1.1 2006/05/19 16:55:01 paulf 00008 * first inclusion 00008 * 00009 * Revision 1.1 2001/03/31 00:46:48 lombard 00010 * Initial revision 00011 * 00012 * 00013 * 00014 */ 00015 00016 /* fft99.f -- translated by f2c (version 19991025).*/ 00017 /* 00018 * Hand-editted after f2c to make it look a little more like C. 00019 * But there is no hiding the fact that this was once Fortran. 00020 * Read the comments with caution! Pete Lombard, January 2001 00021 */ 00022 #include <math.h> 00023 #include "fft99.h" 00024 00025 /* Internal function prototypes */ 00026 static void fact(long n, long *ifax); 00027 static void cftrig(long n, double *trigs); 00028 static void vpassm(double *a, double *b, double *c, double *d, 00029 double *trigs, long inc1, long inc2, long inc3, long inc4, 00030 long lot, long n, long ifac, long la); 00031 static void fft99a(double *a, double *work, double *trigs, long inc, long jump, 00032 long n, long lot); 00033 static void fft99b(double *work, double *a, double *trigs, long inc, long jump, 00034 long n, long lot); 00035 static void fax(long *ifax, long n); 00036 static void fftrig(double *trigs, long n); 00037 00038 00039 /* Function definitions */ 00040 void cfft99(double *a, double *work, double *trigs, long *ifax, long inc, 00041 long jump, long n, long lot, long isign) 00042 { 00043 /* Local variables */ 00044 static long nfax, i, j, k, l, m, ibase, jbase; 00045 static double himag, hreal; 00046 static long ilast, i1, i2, la, nh, nn; 00047 static long igo, ink, jnk, jum, jst; 00048 00049 00050 /* PURPOSE PERFORMS MULTIPLE FAST FOURIER TRANSFORMS. THIS PACKAGE 00051 * WILL PERFORM A NUMBER OF SIMULTANEOUS COMPLEX PERIODIC 00052 * FOURIER TRANSFORMS OR CORRESPONDING INVERSE TRANSFORMS. 00053 * THAT IS, GIVEN A SET OF COMPLEX GRIDPOINT VECTORS, THE 00054 * PACKAGE RETURNS A SET OF COMPLEX FOURIER 00055 * COEFFICIENT VECTORS, OR VICE VERSA. THE LENGTH OF THE 00056 * TRANSFORMS MUST BE A NUMBER GREATER THAN 1 THAT HAS 00057 * NO PRIME FACTORS OTHER THAN 2, 3, AND 5. 00058 * 00059 * THE PACKAGE CFFT99 CONTAINS SEVERAL USER-LEVEL ROUTINES: 00060 * 00061 * CFTFAX 00062 * AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE 00063 * BEFORE A SEQUENCE OF CALLS TO CFFT99 00064 * (PROVIDED THAT N IS NOT CHANGED). 00065 * 00066 * CFFT99 00067 * THE ACTUAL TRANSFORM ROUTINE ROUTINE, CABABLE OF 00068 * PERFORMING BOTH THE TRANSFORM AND ITS INVERSE. 00069 * HOWEVER, AS THE TRANSFORMS ARE NOT NORMALIZED, 00070 * THE APPLICATION OF A TRANSFORM FOLLOWED BY ITS 00071 * INVERSE WILL YIELD THE ORIGINAL VALUES MULTIPLIED 00072 * BY N. 00073 * 00074 * 00075 * 00076 * USAGE LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 0, 00077 * Q .GE. 0, AND R .GE. 0. THEN A TYPICAL SEQUENCE OF 00078 * CALLS TO TRANSFORM A GIVEN SET OF COMPLEX VECTORS OF 00079 * LENGTH N TO A SET OF (UNSCALED) COMPLEX FOURIER 00080 * COEFFICIENT VECTORS OF LENGTH N IS 00081 * 00082 * long ifax[13]; 00083 * double trigs[2*n]; 00084 * double a[...], work[...] 00085 * 00086 * cftfax (n, ifax, trigs); 00087 * cfft99 (a,work,trigs,ifax,inc,jump,n,lot,isign); 00088 * 00089 * THE OUTPUT VECTORS OVERWRITE THE INPUT VECTORS, AND 00090 * THESE ARE STORED IN A. WITH APPROPRIATE CHOICES FOR 00091 * THE OTHER ARGUMENTS, THESE VECTORS MAY BE CONSIDERED 00092 * EITHER THE ROWS OR THE COLUMNS OF THE ARRAY A. 00093 * SEE THE INDIVIDUAL WRITE-UPS FOR CFTFAX AND 00094 * CFFT99 BELOW, FOR A DETAILED DESCRIPTION OF THE 00095 * ARGUMENTS. 00096 * 00097 * HISTORY THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN 00098 * NOVEMBER, 1978. IT WAS MODIFIED, DOCUMENTED, AND TESTED 00099 * FOR NCAR BY RUSS REW IN SEPTEMBER, 1980. IT WAS 00100 * FURTHER MODIFIED FOR THE FULLY COMPLEX CASE BY DAVE 00101 * FULKER IN NOVEMBER, 1980. 00102 * 00103 * ----------------------------------------------------------------------- 00104 * 00105 * cftfax (n,ifax,trigs) 00106 * 00107 * PURPOSE A SET-UP ROUTINE FOR CFFT99. IT NEED ONLY BE 00108 * CALLED ONCE BEFORE A SEQUENCE OF CALLS TO CFFT99, 00109 * PROVIDED THAT N IS NOT CHANGED. 00110 * 00111 * ARGUMENT long ifax[13]; double trigs[2*n]; 00112 * DIMENSIONS 00113 * 00114 * ARGUMENTS 00115 * 00116 * ON INPUT N 00117 * AN EVEN NUMBER GREATER THAN 1 THAT HAS NO PRIME FACTOR 00118 * GREATER THAN 5. N IS THE LENGTH OF THE TRANSFORMS (SEE 00119 * THE DOCUMENTATION FOR CFFT99 FOR THE DEFINITION OF 00120 * THE TRANSFORMS). 00121 * 00122 * IFAX 00123 * AN INTEGER ARRAY. THE NUMBER OF ELEMENTS ACTUALLY USED 00124 * WILL DEPEND ON THE FACTORIZATION OF N. DIMENSIONING 00125 * IFAX FOR 13 SUFFICES FOR ALL N LESS THAN 1 MILLION. 00126 * 00127 * TRIGS 00128 * AN ARRAY OF DIMENSION 2*N 00129 * 00130 * ON OUTPUT IFAX 00131 * CONTAINS THE FACTORIZATION OF N. IFAX[0] IS THE 00132 * NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED 00133 * IN IFAX[1],IFAX[2],... IF N HAS ANY PRIME FACTORS 00134 * GREATER THAN 5, IFAX[0] IS SET TO -99. 00135 * 00136 * TRIGS 00137 * AN ARRAY OF TRIGONOMETRIC FUNCTION VALUES SUBSEQUENTLY 00138 * USED BY THE CFT ROUTINES. 00139 * 00140 * ----------------------------------------------------------------------- 00141 * 00142 * cfft99 (a,work,trigs,ifax,inc,jump,n,lot,isign) 00143 * 00144 * PURPOSE PERFORM A NUMBER OF SIMULTANEOUS (UNNORMALIZED) COMPLEX 00145 * PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE 00146 * TRANSFORMS. GIVEN A SET OF COMPLEX GRIDPOINT 00147 * VECTORS, THE PACKAGE RETURNS A SET OF 00148 * COMPLEX FOURIER COEFFICIENT VECTORS, OR VICE 00149 * VERSA. THE LENGTH OF THE TRANSFORMS MUST BE A 00150 * NUMBER HAVING NO PRIME FACTORS OTHER THAN 00151 * 2, 3, AND 5. THIS ROUTINE IS 00152 * OPTIMIZED FOR USE ON THE CRAY-1. 00153 * 00154 * ARGUMENT double a[2*n*inc+(lot-1)*jump], work[2*n*lot] 00155 * DIMENSIONS double trigs[2*n], long ifax[13] 00156 * 00157 * ARGUMENTS 00158 * 00159 * ON INPUT A 00160 * A COMPLEX ARRAY OF LENGTH N*INC+(LOT-1)*JUMP CONTAINING 00161 * THE INPUT GRIDPOINT OR COEFFICIENT VECTORS. THIS ARRAY 00162 * OVERWRITTEN BY THE RESULTS. 00163 * 00164 * N.B. ALTHOUGH THE ARRAY A IS USUALLY CONSIDERED TO BE OF 00165 * TYPE COMPLEX IN THE CALLING PROGRAM, IT IS TREATED AS 00166 * REAL WITHIN THE TRANSFORM PACKAGE. THIS REQUIRES THAT 00167 * SUCH TYPE CONFLICTS ARE PERMITTED IN THE USER^S 00168 * ENVIRONMENT, AND THAT THE STORAGE OF COMPLEX NUMBERS 00169 * MATCHES THE ASSUMPTIONS OF THIS ROUTINE. THIS ROUTINE 00170 * ASSUMES THAT THE REAL AND IMAGINARY PORTIONS OF A 00171 * COMPLEX NUMBER OCCUPY ADJACENT ELEMENTS OF MEMORY. IF 00172 * THESE CONDITIONS ARE NOT MET, THE USER MUST TREAT THE 00173 * ARRAY A AS REAL (AND OF TWICE THE ABOVE LENGTH), AND 00174 * WRITE THE CALLING PROGRAM TO TREAT THE REAL AND 00175 * IMAGINARY PORTIONS EXPLICITLY. 00176 * 00177 * WORK 00178 * A COMPLEX WORK ARRAY OF LENGTH N*LOT OR A REAL ARRAY 00179 * OF LENGTH 2*N*LOT. SEE N.B. ABOVE. 00180 * 00181 * TRIGS 00182 * AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST. 00183 * 00184 * IFAX 00185 * AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST. 00186 * 00187 * 00188 * N.B. IN THE FOLLOWING ARGUMENTS, INCREMENTS ARE MEASURED 00189 * IN WORD PAIRS, BECAUSE EACH COMPLEX ELEMENT IS ASSUMED 00190 * TO OCCUPY AN ADJACENT PAIR OF WORDS IN MEMORY. 00191 * 00192 * INC 00193 * THE INCREMENT (IN WORD PAIRS) BETWEEN SUCCESSIVE ELEMENT 00194 * OF EACH (COMPLEX) GRIDPOINT OR COEFFICIENT VECTOR 00195 * (E.G. INC=1 FOR CONSECUTIVELY STORED DATA). 00196 * 00197 * JUMP 00198 * THE INCREMENT (IN WORD PAIRS) BETWEEN THE FIRST ELEMENTS 00199 * OF SUCCESSIVE DATA OR COEFFICIENT VECTORS. ON THE CRAY- 00200 * TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8 00201 * (TO AVOID MEMORY BANK CONFLICTS). FOR CLARIFICATION OF 00202 * INC AND JUMP, SEE THE EXAMPLES BELOW. 00203 * 00204 * N 00205 * THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF 00206 * TRANSFORMS, BELOW). 00207 * 00208 * LOT 00209 * THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY. 00210 * 00211 * ISIGN 00212 * = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER 00213 * COEFFICIENTS. 00214 * = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO 00215 * GRIDPOINT VALUES. 00216 * 00217 * ON OUTPUT A 00218 * IF ISIGN = -1, AND LOT GRIDPOINT VECTORS ARE SUPPLIED, 00219 * EACH CONTAINING THE COMPLEX SEQUENCE: 00220 * 00221 * G(0),G(1), ... ,G(N-1) (N COMPLEX VALUES) 00222 * 00223 * THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH 00224 * CONTAINING THE CORRESPONDING N COEFFICIENT VALUES: 00225 * 00226 * C(0),C(1), ... ,C(N-1) (N COMPLEX VALUES) 00227 * 00228 * DEFINED BY: 00229 * C(K) = SUM(J=0,...,N-1)( G(J)*EXP(-2*I*J*K*PI/N) ) 00230 * WHERE I = SQRT(-1) 00231 * 00232 * 00233 * IF ISIGN = +1, AND LOT COEFFICIENT VECTORS ARE SUPPLIED, 00234 * EACH CONTAINING THE COMPLEX SEQUENCE: 00235 * 00236 * C(0),C(1), ... ,C(N-1) (N COMPLEX VALUES) 00237 * 00238 * THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH 00239 * CONTAINING THE CORRESPONDING N GRIDPOINT VALUES: 00240 * 00241 * G(0),G(1), ... ,G(N-1) (N COMPLEX VALUES) 00242 * 00243 * DEFINED BY: 00244 * G(J) = SUM(K=0,...,N-1)( G(K)*EXP(+2*I*J*K*PI/N) ) 00245 * WHERE I = SQRT(-1) 00246 * 00247 * 00248 * A CALL WITH ISIGN=-1 FOLLOWED BY A CALL WITH ISIGN=+1 00249 * (OR VICE VERSA) RETURNS THE ORIGINAL DATA, MULTIPLIED 00250 * BY THE FACTOR N. 00251 * 00252 * 00253 * EXAMPLE GIVEN A 64 BY 9 GRID OF COMPLEX VALUES, STORED IN 00254 * A 66 BY 9 COMPLEX ARRAY, A, COMPUTE THE TWO DIMENSIONAL 00255 * FOURIER TRANSFORM OF THE GRID. FROM TRANSFORM THEORY, 00256 * IT IS KNOWN THAT A TWO DIMENSIONAL TRANSFORM CAN BE 00257 * OBTAINED BY FIRST TRANSFORMING THE GRID ALONG ONE 00258 * DIRECTION, THEN TRANSFORMING THESE RESULTS ALONG THE 00259 * ORTHOGONAL DIRECTION. 00260 * 00261 * COMPLEX A(66,9), WORK(64,9) 00262 * REAL TRIGS1(128), TRIGS2(18) 00263 * INTEGER IFAX1(13), IFAX2(13) 00264 * 00265 * SET UP THE IFAX AND TRIGS ARRAYS FOR EACH DIRECTION: 00266 * 00267 * CALL CFTFAX(64, IFAX1, TRIGS1) 00268 * CALL CFTFAX( 9, IFAX2, TRIGS2) 00269 * 00270 * IN THIS CASE, THE COMPLEX VALUES OF THE GRID ARE 00271 * STORED IN MEMORY AS FOLLOWS (USING U AND V TO 00272 * DENOTE THE REAL AND IMAGINARY COMPONENTS, AND 00273 * ASSUMING CONVENTIONAL FORTRAN STORAGE): 00274 * 00275 * U(1,1), V(1,1), U(2,1), V(2,1), ... U(64,1), V(64,1), 4 NULLS, 00276 * 00277 * U(1,2), V(1,2), U(2,2), V(2,2), ... U(64,2), V(64,2), 4 NULLS, 00278 * 00279 * . . . . . . . . 00280 * . . . . . . . . 00281 * . . . . . . . . 00282 * 00283 * U(1,9), V(1,9), U(2,9), V(2,9), ... U(64,9), V(64,9), 4 NULLS. 00284 * 00285 * WE CHOOSE (ARBITRARILY) TO TRANSORM FIRST ALONG THE 00286 * DIRECTION OF THE FIRST SUBSCRIPT. THUS WE DEFINE 00287 * THE LENGTH OF THE TRANSFORMS, N, TO BE 64, THE 00288 * NUMBER OF TRANSFORMS, LOT, TO BE 9, THE INCREMENT 00289 * BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 1, 00290 * AND THE INCREMENT BETWEEN THE STARTING POINTS 00291 * FOR EACH TRANSFORM, JUMP, TO BE 66 (THE FIRST 00292 * DIMENSION OF A). 00293 * 00294 * CALL CFFT99( A, WORK, TRIGS1, IFAX1, 1, 66, 64, 9, -1) 00295 * 00296 * TO TRANSFORM ALONG THE DIRECTION OF THE SECOND SUBSCRIPT 00297 * THE ROLES OF THE INCREMENTS ARE REVERSED. THUS WE DEFIN 00298 * THE LENGTH OF THE TRANSFORMS, N, TO BE 9, THE 00299 * NUMBER OF TRANSFORMS, LOT, TO BE 64, THE INCREMENT 00300 * BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 66, 00301 * AND THE INCREMENT BETWEEN THE STARTING POINTS 00302 * FOR EACH TRANSFORM, JUMP, TO BE 1 00303 * 00304 * CALL CFFT99( A, WORK, TRIGS2, IFAX2, 66, 1, 9, 64, -1) 00305 * 00306 * THESE TWO SEQUENTIAL STEPS RESULTS IN THE TWO-DIMENSIONA 00307 * FOURIER COEFFICIENT ARRAY OVERWRITING THE INPUT 00308 * GRIDPOINT ARRAY, A. THE SAME TWO STEPS APPLIED AGAIN 00309 * WITH ISIGN = +1 WOULD RESULT IN THE RECONSTRUCTION OF 00310 * THE GRIDPOINT ARRAY (MULTIPLIED BY A FACTOR OF 64*9). 00311 * 00312 * 00313 * ----------------------------------------------------------------------- 00314 */ 00315 00316 /* function ^cfft99^ - MULTIPLE FAST COMPLEX FOURIER TRANSFORM 00317 00318 * A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA 00319 * WORK IS AN AREA OF SIZE N*LOT 00320 * TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES 00321 * IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N 00322 * INC IS THE INCREMENT WITHIN EACH DATA #VECTOR# 00323 * (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) 00324 * JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR 00325 * N IS THE LENGTH OF THE DATA VECTORS 00326 * LOT IS THE NUMBER OF DATA VECTORS 00327 * ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT 00328 * = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL 00329 */ 00330 00331 /* VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN PARALLEL. */ 00332 00333 /* Parameter adjustments */ 00334 --work; 00335 --a; 00336 00337 /* Function Body */ 00338 nn = n + n; 00339 ink = inc + inc; 00340 jum = jump + jump; 00341 nfax = ifax[0]; 00342 jnk = 2; 00343 jst = 2; 00344 if (isign >= 0) { 00345 goto L30; 00346 } 00347 00348 /* THE INNERMOST TEMPERTON ROUTINES HAVE NO FACILITY FOR THE */ 00349 /* FORWARD (ISIGN = -1) TRANSFORM. THEREFORE, THE INPUT MUST BE */ 00350 /* REARRANGED AS FOLLOWS: */ 00351 00352 /* THE ORDER OF EACH INPUT VECTOR, */ 00353 00354 /* G(0), G(1), G(2), ... , G(N-2), G(N-1) */ 00355 00356 /* IS REVERSED (EXCLUDING G(0)) TO YIELD */ 00357 00358 /* G(0), G(N-1), G(N-2), ... , G(2), G(1). */ 00359 00360 /* WITHIN THE TRANSFORM, THE CORRESPONDING EXPONENTIAL MULTIPLIER */ 00361 /* IS THEN PRECISELY THE CONJUGATE OF THAT FOR THE NORMAL */ 00362 /* ORDERING. THUS THE FORWARD (ISIGN = -1) TRANSFORM IS */ 00363 /* ACCOMPLISHED */ 00364 00365 /* FOR NFAX ODD, THE INPUT MUST BE TRANSFERRED TO THE WORK ARRAY, */ 00366 /* AND THE REARRANGEMENT CAN BE DONE DURING THE MOVE. */ 00367 00368 jnk = -2; 00369 jst = nn - 2; 00370 if (nfax % 2 == 1) { 00371 goto L40; 00372 } 00373 00374 /* FOR NFAX EVEN, THE REARRANGEMENT MUST BE APPLIED DIRECTLY TO */ 00375 /* THE INPUT ARRAY. THIS CAN BE DONE BY SWAPPING ELEMENTS. */ 00376 00377 ibase = 1; 00378 ilast = (n - 1) * ink; 00379 nh = n / 2; 00380 for (l = 1; l <= lot; ++l) { 00381 i1 = ibase + ink; 00382 i2 = ibase + ilast; 00383 /* DIR$ IVDEP */ 00384 for (m = 1; m <= nh; ++m) { 00385 /* SWAP REAL AND IMAGINARY PORTIONS */ 00386 hreal = a[i1]; 00387 himag = a[i1 + 1]; 00388 a[i1] = a[i2]; 00389 a[i1 + 1] = a[i2 + 1]; 00390 a[i2] = hreal; 00391 a[i2 + 1] = himag; 00392 i1 += ink; 00393 i2 -= ink; 00394 /* L10: */ 00395 } 00396 ibase += jum; 00397 /* L20: */ 00398 } 00399 goto L100; 00400 00401 L30: 00402 if (nfax % 2 == 0) { 00403 goto L100; 00404 } 00405 00406 L40: 00407 00408 /* DURING THE TRANSFORM PROCESS, NFAX STEPS ARE TAKEN, AND THE */ 00409 /* RESULTS ARE STORED ALTERNATELY IN WORK AND IN A. IF NFAX IS */ 00410 /* ODD, THE INPUT DATA ARE FIRST MOVED TO WORK SO THAT THE FINAL */ 00411 /* RESULT (AFTER NFAX STEPS) IS STORED IN ARRAY A. */ 00412 00413 ibase = 1; 00414 jbase = 1; 00415 for (l = 1; l <= lot; ++l) { 00416 /* MOVE REAL AND IMAGINARY PORTIONS OF ELEMENT ZERO */ 00417 work[jbase] = a[ibase]; 00418 work[jbase + 1] = a[ibase + 1]; 00419 i = ibase + ink; 00420 j = jbase + jst; 00421 /* DIR$ IVDEP */ 00422 for (m = 2; m <= n; ++m) { 00423 /* MOVE REAL AND IMAGINARY PORTIONS OF OTHER ELEMENTS (POSSIBLY IN */ 00424 /* REVERSE ORDER, DEPENDING ON JST AND JNK) */ 00425 work[j] = a[i]; 00426 work[j + 1] = a[i + 1]; 00427 i += ink; 00428 j += jnk; 00429 /* L50: */ 00430 } 00431 ibase += jum; 00432 jbase += nn; 00433 /* L60: */ 00434 } 00435 00436 L100: 00437 00438 /* PERFORM THE TRANSFORM PASSES, ONE PASS FOR EACH FACTOR. DURING */ 00439 /* EACH PASS THE DATA ARE MOVED FROM A TO WORK OR FROM WORK TO A. */ 00440 00441 /* FOR NFAX EVEN, THE FIRST PASS MOVES FROM A TO WORK */ 00442 igo = 110; 00443 /* FOR NFAX ODD, THE FIRST PASS MOVES FROM WORK TO A */ 00444 if (nfax % 2 == 1) { 00445 igo = 120; 00446 } 00447 la = 1; 00448 for (k = 1; k <= nfax; ++k) { 00449 if (igo == 120) { 00450 goto L120; 00451 } 00452 /* L110: */ 00453 vpassm(&a[1], &a[2], &work[1], &work[2], trigs, ink, 2, jum, nn, lot, n, 00454 ifax[k], la); 00455 igo = 120; 00456 goto L130; 00457 L120: 00458 vpassm(&work[1], &work[2], &a[1], &a[2], trigs, 2, ink, nn, 00459 jum, lot, n, ifax[k], la); 00460 igo = 110; 00461 L130: 00462 la *= ifax[k]; 00463 /* L140: */ 00464 } 00465 00466 /* AT THIS POINT THE FINAL TRANSFORM RESULT IS STORED IN A. */ 00467 00468 return; 00469 } /* cfft99 */ 00470 00471 void cftfax(long n, long *ifax, double *trigs) 00472 { 00473 /* Local variables */ 00474 static long k; 00475 00476 /* THIS ROUTINE WAS MODIFIED FROM TEMPERTON^S ORIGINAL */ 00477 /* BY DAVE FULKER. IT NO LONGER PRODUCES FACTORS IN ASCENDING */ 00478 /* ORDER, AND THERE ARE NONE OF THE ORIGINAL #MODE# OPTIONS. */ 00479 00480 /* ON INPUT N */ 00481 /* THE LENGTH OF EACH COMPLEX TRANSFORM TO BE PERFORMED */ 00482 00483 /* N MUST BE GREATER THAN 1 AND CONTAIN NO PRIME */ 00484 /* FACTORS GREATER THAN 5. */ 00485 00486 /* ON OUTPUT IFAX */ 00487 /* IFAX(1) */ 00488 /* THE NUMBER OF FACTORS CHOSEN OR -99 IN CASE OF ERROR */ 00489 /* IFAX(2) THRU IFAX( IFAX(1)+1 ) */ 00490 /* THE FACTORS OF N IN THE FOLLOWIN ORDER: APPEARING */ 00491 /* FIRST ARE AS MANY FACTORS OF 4 AS CAN BE OBTAINED. */ 00492 /* SUBSEQUENT FACTORS ARE PRIMES, AND APPEAR IN */ 00493 /* ASCENDING ORDER, EXCEPT FOR MULTIPLE FACTORS. */ 00494 00495 /* TRIGS */ 00496 /* 2N SIN AND COS VALUES FOR USE BY THE TRANSFORM ROUTINE */ 00497 00498 /* Function Body */ 00499 fact(n, ifax); 00500 k = ifax[0]; 00501 if (k < 1 || ifax[k] > 5) { 00502 ifax[0] = -99; 00503 } 00504 if (ifax[0] <= 0) 00505 return; 00506 00507 cftrig(n, trigs); 00508 return; 00509 } /* cftfax */ 00510 00511 static void fact(long n, long *ifax) 00512 { 00513 /* Local variables */ 00514 static long k, l, nn, inc, max; 00515 00516 /* FACTORIZATION ROUTINE THAT FIRST EXTRACTS ALL FACTORS OF 4 */ 00517 /* Function Body */ 00518 if (n > 1) { 00519 goto L10; 00520 } 00521 ifax[0] = 0; 00522 if (n < 1) { 00523 ifax[0] = -99; 00524 } 00525 return; 00526 L10: 00527 nn = n; 00528 k = 0; 00529 /* TEST FOR FACTORS OF 4 */ 00530 L20: 00531 if (nn % 4 != 0) { 00532 goto L30; 00533 } 00534 ++k; 00535 ifax[k] = 4; 00536 nn /= 4; 00537 if (nn == 1) { 00538 goto L80; 00539 } 00540 goto L20; 00541 /* TEST FOR EXTRA FACTOR OF 2 */ 00542 L30: 00543 if (nn % 2 != 0) { 00544 goto L40; 00545 } 00546 ++k; 00547 ifax[k] = 2; 00548 nn /= 2; 00549 if (nn == 1) { 00550 goto L80; 00551 } 00552 /* TEST FOR FACTORS OF 3 */ 00553 L40: 00554 if (nn % 3 != 0) { 00555 goto L50; 00556 } 00557 ++k; 00558 ifax[k] = 3; 00559 nn /= 3; 00560 if (nn == 1) { 00561 goto L80; 00562 } 00563 goto L40; 00564 /* NOW FIND REMAINING FACTORS */ 00565 L50: 00566 l = 5; 00567 max = sqrt((double) nn); 00568 inc = 2; 00569 /* INC ALTERNATELY TAKES ON VALUES 2 AND 4 */ 00570 L60: 00571 if (nn % l != 0) { 00572 goto L70; 00573 } 00574 ++k; 00575 ifax[k] = l; 00576 nn /= l; 00577 if (nn == 1) { 00578 goto L80; 00579 } 00580 goto L60; 00581 L70: 00582 if (l > max) { 00583 goto L75; 00584 } 00585 l += inc; 00586 inc = 6 - inc; 00587 goto L60; 00588 L75: 00589 ++k; 00590 ifax[k] = nn; 00591 L80: 00592 ifax[0] = k; 00593 /* IFAX(0) NOW CONTAINS NUMBER OF FACTORS */ 00594 return; 00595 } /* fact */ 00596 00597 static void cftrig(long n, double *trigs) 00598 { 00599 /* Local variables */ 00600 static long i, n2; 00601 static double angle, del; 00602 00603 /* Function Body */ 00604 del = (PI + PI) / (double) (n); 00605 n2 = n + n; 00606 for (i = 0; i < n2; i += 2) { 00607 angle = (double) i * 0.5 * del; 00608 trigs[i] = cos(angle); 00609 trigs[i + 1] = sin(angle); 00610 } 00611 return; 00612 } /* cftrig */ 00613 00614 static void vpassm(double *a, double *b, double *c, double *d, 00615 double *trigs, long inc1, long inc2, long inc3, long inc4, 00616 long lot, long n, long ifac, long la) 00617 { 00618 /* Initialized data */ 00619 00620 static double sin36 = 0.587785252292473; 00621 static double cos36 = 0.809016994374947; 00622 static double sin72 = 0.951056516295154; 00623 static double cos72 = 0.309016994374947; 00624 static double sin60 = 0.866025403784437; 00625 00626 /* Local variables */ 00627 static long iink, jink, jump, i, j, k, l, m, ibase, jbase; 00628 static double c1, c2, c3, c4, s1, s2, s3, s4; 00629 static long ia, ja, ib, jb, kb, ic, jc, kc, id, jd, kd, ie, je, ke, 00630 la1, ijk, igo; 00631 00632 00633 /* SUBROUTINE ^VPASSM^ - MULTIPLE VERSION OF ^VPASSA^ */ 00634 /* PERFORMS ONE PASS THROUGH DATA */ 00635 /* AS PART OF MULTIPLE COMPLEX (INVERSE) FFT ROUTINE */ 00636 /* A IS FIRST REAL INPUT VECTOR */ 00637 /* B IS FIRST IMAGINARY INPUT VECTOR */ 00638 /* C IS FIRST REAL OUTPUT VECTOR */ 00639 /* D IS FIRST IMAGINARY OUTPUT VECTOR */ 00640 /* TRIGS IS PRECALCULATED TABLE OF SINES \ COSINES */ 00641 /* INC1 IS ADDRESSING INCREMENT FOR A AND B */ 00642 /* INC2 IS ADDRESSING INCREMENT FOR C AND D */ 00643 /* INC3 IS ADDRESSING INCREMENT BETWEEN A^S \ B^S */ 00644 /* INC4 IS ADDRESSING INCREMENT BETWEEN C^S \ D^S */ 00645 /* LOT IS THE NUMBER OF VECTORS */ 00646 /* N IS LENGTH OF VECTORS */ 00647 /* IFAC IS CURRENT FACTOR OF N */ 00648 /* LA IS PRODUCT OF PREVIOUS FACTORS */ 00649 00650 /* Parameter adjustments */ 00651 --d; 00652 --c; 00653 --b; 00654 --a; 00655 00656 /* Function Body */ 00657 00658 m = n / ifac; 00659 iink = m * inc1; 00660 jink = la * inc2; 00661 jump = (ifac - 1) * jink; 00662 ibase = 0; 00663 jbase = 0; 00664 igo = ifac - 1; 00665 if (igo > 4) { 00666 return; 00667 } 00668 switch ((int)igo) { 00669 case 1: goto L10; 00670 case 2: goto L50; 00671 case 3: goto L90; 00672 case 4: goto L130; 00673 } 00674 00675 /* CODING FOR FACTOR 2 */ 00676 00677 L10: 00678 ia = 1; 00679 ja = 1; 00680 ib = ia + iink; 00681 jb = ja + jink; 00682 for (l = 1; l <= la; ++l) { 00683 i = ibase; 00684 j = jbase; 00685 /* DIR$ IVDEP */ 00686 for (ijk = 1; ijk <= lot; ++ijk) { 00687 c[ja + j] = a[ia + i] + a[ib + i]; 00688 d[ja + j] = b[ia + i] + b[ib + i]; 00689 c[jb + j] = a[ia + i] - a[ib + i]; 00690 d[jb + j] = b[ia + i] - b[ib + i]; 00691 i += inc3; 00692 j += inc4; 00693 /* L15: */ 00694 } 00695 ibase += inc1; 00696 jbase += inc2; 00697 /* L20: */ 00698 } 00699 if (la == m) { 00700 return; 00701 } 00702 la1 = la + 1; 00703 jbase += jump; 00704 for (k = la1; la < 0 ? k >= m : k <= m; k += la) { 00705 kb = k + k - 2; 00706 c1 = trigs[kb]; 00707 s1 = trigs[kb + 1]; 00708 for (l = 1; l <= la; ++l) { 00709 i = ibase; 00710 j = jbase; 00711 /* DIR$ IVDEP */ 00712 for (ijk = 1; ijk <= lot; ++ijk) { 00713 c[ja + j] = a[ia + i] + a[ib + i]; 00714 d[ja + j] = b[ia + i] + b[ib + i]; 00715 c[jb + j] = c1 * (a[ia + i] - a[ib + i]) - s1 * 00716 (b[ia +i] - b[ib + i]); 00717 d[jb + j] = s1 * (a[ia + i] - a[ib + i]) + c1 * 00718 (b[ia +i] - b[ib + i]); 00719 i += inc3; 00720 j += inc4; 00721 /* L25: */ 00722 } 00723 ibase += inc1; 00724 jbase += inc2; 00725 /* L30: */ 00726 } 00727 jbase += jump; 00728 /* L40: */ 00729 } 00730 return; 00731 00732 /* CODING FOR FACTOR 3 */ 00733 00734 L50: 00735 ia = 1; 00736 ja = 1; 00737 ib = ia + iink; 00738 jb = ja + jink; 00739 ic = ib + iink; 00740 jc = jb + jink; 00741 for (l = 1; l <= la; ++l) { 00742 i = ibase; 00743 j = jbase; 00744 /* DIR$ IVDEP */ 00745 for (ijk = 1; ijk <= lot; ++ijk) { 00746 c[ja + j] = a[ia + i] + (a[ib + i] + a[ic + i]); 00747 d[ja + j] = b[ia + i] + (b[ib + i] + b[ic + i]); 00748 c[jb + j] = a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 - 00749 sin60 * (b[ib + i] - b[ic + i]); 00750 c[jc + j] = a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 + 00751 sin60 * (b[ib + i] - b[ic + i]); 00752 d[jb + j] = b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 + 00753 sin60 * (a[ib + i] - a[ic + i]); 00754 d[jc + j] = b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 - 00755 sin60 * (a[ib + i] - a[ic + i]); 00756 i += inc3; 00757 j += inc4; 00758 /* L55: */ 00759 } 00760 ibase += inc1; 00761 jbase += inc2; 00762 /* L60: */ 00763 } 00764 if (la == m) { 00765 return; 00766 } 00767 la1 = la + 1; 00768 jbase += jump; 00769 for (k = la1; la < 0 ? k >= m : k <= m; k += la) { 00770 kb = k + k - 2; 00771 kc = kb + kb; 00772 c1 = trigs[kb]; 00773 s1 = trigs[kb + 1]; 00774 c2 = trigs[kc]; 00775 s2 = trigs[kc + 1]; 00776 for (l = 1; l <= la; ++l) { 00777 i = ibase; 00778 j = jbase; 00779 /* DIR$ IVDEP */ 00780 for (ijk = 1; ijk <= lot; ++ijk) { 00781 c[ja + j] = a[ia + i] + (a[ib + i] + a[ic + i]); 00782 d[ja + j] = b[ia + i] + (b[ib + i] + b[ic + i]); 00783 c[jb + j] = c1 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 - 00784 sin60 * (b[ib + i] - b[ic + i])) - 00785 s1 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 + 00786 sin60 * (a[ib + i] - a[ic + i])); 00787 d[jb + j] = s1 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 - 00788 sin60 * (b[ib + i] - b[ic + i])) + 00789 c1 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 + 00790 sin60 * (a[ib + i] - a[ic + i])); 00791 c[jc + j] = c2 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 + 00792 sin60 * (b[ib + i] - b[ic + i])) - 00793 s2 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 - 00794 sin60 * (a[ib + i] - a[ic + i])); 00795 d[jc + j] = s2 * (a[ia + i] - (a[ib + i] + a[ic + i]) * 0.5 + 00796 sin60 * (b[ib + i] - b[ic + i])) + 00797 c2 * (b[ia + i] - (b[ib + i] + b[ic + i]) * 0.5 - 00798 sin60 * (a[ib + i] - a[ic + i])); 00799 i += inc3; 00800 j += inc4; 00801 /* L65: */ 00802 } 00803 ibase += inc1; 00804 jbase += inc2; 00805 /* L70: */ 00806 } 00807 jbase += jump; 00808 /* L80: */ 00809 } 00810 return; 00811 00812 /* CODING FOR FACTOR 4 */ 00813 00814 L90: 00815 ia = 1; 00816 ja = 1; 00817 ib = ia + iink; 00818 jb = ja + jink; 00819 ic = ib + iink; 00820 jc = jb + jink; 00821 id = ic + iink; 00822 jd = jc + jink; 00823 for (l = 1; l <= la; ++l) { 00824 i = ibase; 00825 j = jbase; 00826 /* DIR$ IVDEP */ 00827 for (ijk = 1; ijk <= lot; ++ijk) { 00828 c[ja + j] = a[ia + i] + a[ic + i] + (a[ib + i] + a[id + i]); 00829 c[jc + j] = a[ia + i] + a[ic + i] - (a[ib + i] + a[id + i]); 00830 d[ja + j] = b[ia + i] + b[ic + i] + (b[ib + i] + b[id + i]); 00831 d[jc + j] = b[ia + i] + b[ic + i] - (b[ib + i] + b[id + i]); 00832 c[jb + j] = a[ia + i] - a[ic + i] - (b[ib + i] - b[id + i]); 00833 c[jd + j] = a[ia + i] - a[ic + i] + (b[ib + i] - b[id + i]); 00834 d[jb + j] = b[ia + i] - b[ic + i] + (a[ib + i] - a[id + i]); 00835 d[jd + j] = b[ia + i] - b[ic + i] - (a[ib + i] - a[id + i]); 00836 i += inc3; 00837 j += inc4; 00838 /* L95: */ 00839 } 00840 ibase += inc1; 00841 jbase += inc2; 00842 /* L100: */ 00843 } 00844 if (la == m) { 00845 return; 00846 } 00847 la1 = la + 1; 00848 jbase += jump; 00849 for (k = la1; la < 0 ? k >= m : k <= m; k += la) { 00850 kb = k + k - 2; 00851 kc = kb + kb; 00852 kd = kc + kb; 00853 c1 = trigs[kb]; 00854 s1 = trigs[kb + 1]; 00855 c2 = trigs[kc]; 00856 s2 = trigs[kc + 1]; 00857 c3 = trigs[kd]; 00858 s3 = trigs[kd + 1]; 00859 for (l = 1; l <= la; ++l) { 00860 i = ibase; 00861 j = jbase; 00862 /* DIR$ IVDEP */ 00863 for (ijk = 1; ijk <= lot; ++ijk) { 00864 c[ja + j] = a[ia + i] + a[ic + i] + (a[ib + i] + a[id + i]); 00865 d[ja + j] = b[ia + i] + b[ic + i] + (b[ib + i] + b[id + i]); 00866 c[jc + j] = c2 * (a[ia + i] + a[ic + i] - (a[ib + i] + a[id + i])) - 00867 s2 * (b[ia + i] + b[ic + i] - (b[ib + i] + b[id + i])); 00868 d[jc + j] = s2 * (a[ia + i] + a[ic + i] - (a[ib + i] + a[id + i])) + 00869 c2 * (b[ia + i] + b[ic + i] - (b[ib + i] + b[id + i])); 00870 c[jb + j] = c1 * (a[ia + i] - a[ic + i] - (b[ib + i] - b[id + i])) - 00871 s1 * (b[ia + i] - b[ic + i] + (a[ib + i] - a[id + i])); 00872 d[jb + j] = s1 * (a[ia + i] - a[ic + i] - (b[ib + i] - b[id + i])) + 00873 c1 * (b[ia + i] - b[ic + i] + (a[ib + i] - a[id + i])); 00874 c[jd + j] = c3 * (a[ia + i] - a[ic + i] + (b[ib + i] - b[id + i])) - 00875 s3 * (b[ia + i] - b[ic + i] - (a[ib + i] - a[id + i])); 00876 d[jd + j] = s3 * (a[ia + i] - a[ic + i] + (b[ib + i] - b[id + i])) + 00877 c3 * (b[ia + i] - b[ic + i] - (a[ib + i] - a[id + i])); 00878 i += inc3; 00879 j += inc4; 00880 /* L105: */ 00881 } 00882 ibase += inc1; 00883 jbase += inc2; 00884 /* L110: */ 00885 } 00886 jbase += jump; 00887 /* L120: */ 00888 } 00889 return; 00890 00891 /* CODING FOR FACTOR 5 */ 00892 00893 L130: 00894 ia = 1; 00895 ja = 1; 00896 ib = ia + iink; 00897 jb = ja + jink; 00898 ic = ib + iink; 00899 jc = jb + jink; 00900 id = ic + iink; 00901 jd = jc + jink; 00902 ie = id + iink; 00903 je = jd + jink; 00904 for (l = 1; l <= la; ++l) { 00905 i = ibase; 00906 j = jbase; 00907 /* DIR$ IVDEP */ 00908 for (ijk = 1; ijk <= lot; ++ijk) { 00909 c[ja + j] = a[ia + i] + (a[ib + i] + a[ie + i]) + 00910 (a[ic + i] + a[id + i]); 00911 d[ja + j] = b[ia + i] + (b[ib + i] + b[ie + i]) + 00912 (b[ic + i] + b[id + i]); 00913 c[jb + j] = a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) - 00914 cos36 * (a[ic + i] + a[id + i]) - (sin72 * (b[ib + i] - b[ie + i]) + 00915 sin36 * (b[ic + i] - b[id + i])); 00916 c[je + j] = a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) - 00917 cos36 * (a[ic + i] + a[id + i]) + (sin72 * (b[ib + i] - b[ie + i]) + 00918 sin36 * (b[ic + i] - b[id + i])); 00919 d[jb + j] = b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) - 00920 cos36 * (b[ic + i] + b[id + i]) + (sin72 * (a[ib + i] - a[ie + i]) + 00921 sin36 * (a[ic + i] - a[id + i])); 00922 d[je + j] = b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) - 00923 cos36 * (b[ic + i] + b[id + i]) - (sin72 * (a[ib + i] - a[ie + i]) + 00924 sin36 * (a[ic + i] - a[id + i])); 00925 c[jc + j] = a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) + 00926 cos72 * (a[ic + i] + a[id + i]) - (sin36 * (b[ib + i] - b[ie + i]) - 00927 sin72 * (b[ic + i] - b[id + i])); 00928 c[jd + j] = a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) + 00929 cos72 * (a[ic + i] + a[id + i]) + (sin36 * (b[ib + i] - b[ie + i]) - 00930 sin72 * (b[ic + i] - b[id + i])); 00931 d[jc + j] = b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) + 00932 cos72 * (b[ic + i] + b[id + i]) + (sin36 * (a[ib + i] - a[ie + i]) - 00933 sin72 * (a[ic + i] - a[id + i])); 00934 d[jd + j] = b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) + 00935 cos72 * (b[ic + i] + b[id + i]) - (sin36 * (a[ib + i] - a[ie + i]) - 00936 sin72 * (a[ic + i] - a[id + i])); 00937 i += inc3; 00938 j += inc4; 00939 /* L135: */ 00940 } 00941 ibase += inc1; 00942 jbase += inc2; 00943 /* L140: */ 00944 } 00945 if (la == m) { 00946 return; 00947 } 00948 la1 = la + 1; 00949 jbase += jump; 00950 for (k = la1; la < 0 ? k >= m : k <= m; k += la) { 00951 kb = k + k - 2; 00952 kc = kb + kb; 00953 kd = kc + kb; 00954 ke = kd + kb; 00955 c1 = trigs[kb]; 00956 s1 = trigs[kb + 1]; 00957 c2 = trigs[kc]; 00958 s2 = trigs[kc + 1]; 00959 c3 = trigs[kd]; 00960 s3 = trigs[kd + 1]; 00961 c4 = trigs[ke]; 00962 s4 = trigs[ke + 1]; 00963 for (l = 1; l <= la; ++l) { 00964 i = ibase; 00965 j = jbase; 00966 /* DIR$ IVDEP */ 00967 for (ijk = 1; ijk <= lot; ++ijk) { 00968 c[ja + j] = a[ia + i] + (a[ib + i] + a[ie + i]) + 00969 (a[ic + i] + a[id + i]); 00970 d[ja + j] = b[ia + i] + (b[ib + i] + b[ie + i]) + 00971 (b[ic + i] + b[id + i]); 00972 c[jb + j] = 00973 c1 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) - 00974 cos36 * (a[ic + i] + a[id + i]) - 00975 (sin72 * (b[ib + i] - b[ie + i]) + 00976 sin36 * (b[ic + i] - b[id + i]))) - 00977 s1 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) - 00978 cos36 * (b[ic + i] + b[id + i]) + 00979 (sin72 * (a[ib + i] - a[ie + i]) + 00980 sin36 * (a[ic + i] - a[id + i]))); 00981 d[jb + j] = 00982 s1 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) - 00983 cos36 * (a[ic + i] + a[id + i]) - 00984 (sin72 * (b[ib + i] - b[ie + i]) + 00985 sin36 * (b[ic + i] - b[id + i]))) + 00986 c1 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) - 00987 cos36 * (b[ic + i] + b[id + i]) + 00988 (sin72 * (a[ib + i] - a[ie + i]) + 00989 sin36 * (a[ic + i] - a[id + i]))); 00990 c[je + j] = 00991 c4 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) - 00992 cos36 * (a[ic + i] + a[id + i]) + 00993 (sin72 * (b[ib + i] - b[ie + i]) + 00994 sin36 * (b[ic + i] - b[id + i]))) - 00995 s4 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) - 00996 cos36 * (b[ic + i] + b[id + i]) - 00997 (sin72 * (a[ib + i] - a[ie + i]) + 00998 sin36 * (a[ic + i] - a[id + i]))); 00999 d[je + j] = 01000 s4 * (a[ia + i] + cos72 * (a[ib + i] + a[ie + i]) - 01001 cos36 * (a[ic + i] + a[id + i]) + 01002 (sin72 * (b[ib + i] - b[ie + i]) + 01003 sin36 * (b[ic + i] - b[id + i]))) + 01004 c4 * (b[ia + i] + cos72 * (b[ib + i] + b[ie + i]) - 01005 cos36 * (b[ic + i] + b[id + i]) - 01006 (sin72 * (a[ib + i] - a[ie + i]) + 01007 sin36 * (a[ic + i] - a[id + i]))); 01008 c[jc + j] = 01009 c2 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) + 01010 cos72 * (a[ic + i] + a[id + i]) - 01011 (sin36 * (b[ib + i] - b[ie + i]) - 01012 sin72 * (b[ic + i] - b[id + i]))) - 01013 s2 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) + 01014 cos72 * (b[ic + i] + b[id + i]) + 01015 (sin36 * (a[ib + i] - a[ie + i]) - 01016 sin72 * (a[ic + i] - a[id + i]))); 01017 d[jc + j] = 01018 s2 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) + 01019 cos72 * (a[ic + i] + a[id + i]) - 01020 (sin36 * (b[ib + i] - b[ie + i]) - 01021 sin72 * (b[ic + i] - b[id + i]))) + 01022 c2 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) + 01023 cos72 * (b[ic + i] + b[id + i]) + 01024 (sin36 * (a[ib + i] - a[ie + i]) - 01025 sin72 * (a[ic + i] - a[id + i]))); 01026 c[jd + j] = 01027 c3 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) + 01028 cos72 * (a[ic + i] + a[id + i]) + 01029 (sin36 * (b[ib + i] - b[ie + i]) - 01030 sin72 * (b[ic +i] - b[id + i]))) - 01031 s3 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) + 01032 cos72 * (b[ic + i] + b[id + i]) - 01033 (sin36 * (a[ib + i] - a[ie + i]) - 01034 sin72 * (a[ic + i] - a[id + i]))); 01035 d[jd + j] = 01036 s3 * (a[ia + i] - cos36 * (a[ib + i] + a[ie + i]) + 01037 cos72 * (a[ic + i] + a[id + i]) + 01038 (sin36 * (b[ib + i] - b[ie + i]) - 01039 sin72 * (b[ic + i] - b[id + i]))) + 01040 c3 * (b[ia + i] - cos36 * (b[ib + i] + b[ie + i]) + 01041 cos72 * (b[ic + i] + b[id + i]) - 01042 (sin36 * (a[ib + i] - a[ie + i]) - 01043 sin72 * (a[ic + i] - a[id + i]))); 01044 i += inc3; 01045 j += inc4; 01046 /* L145: */ 01047 } 01048 ibase += inc1; 01049 jbase += inc2; 01050 /* L150: */ 01051 } 01052 jbase += jump; 01053 /* L160: */ 01054 } 01055 return; 01056 } /* vpassm */ 01057 01058 void fft99(double *a, double *work, double *trigs, long *ifax, long inc, 01059 long jump, long n, long lot, long isign) 01060 { 01061 /* Local variables */ 01062 static long nfax, i, j, k, l, m, ibase, jbase; 01063 static long ia, ib, la, nh, nx; 01064 static long igo, ink; 01065 01066 01067 /* PURPOSE PERFORMS MULTIPLE FAST FOURIER TRANSFORMS. THIS PACKAGE 01068 * WILL PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX 01069 * PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE 01070 * TRANSFORMS, I.E. GIVEN A SET OF REAL DATA VECTORS, THE 01071 * PACKAGE RETURNS A SET OF #HALF-COMPLEX# FOURIER 01072 * COEFFICIENT VECTORS, OR VICE VERSA. THE LENGTH OF THE 01073 * TRANSFORMS MUST BE AN EVEN NUMBER GREATER THAN 4 THAT HAS 01074 * NO OTHER FACTORS EXCEPT POSSIBLY POWERS OF 2, 3, AND 5. 01075 * 01076 * THE PACKAGE FFT99F CONTAINS SEVERAL USER-LEVEL ROUTINES: 01077 * 01078 * fftfax 01079 * AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE 01080 * BEFORE A SEQUENCE OF CALLS TO THE FFT ROUTINES 01081 * (PROVIDED THAT N IS NOT CHANGED). 01082 * 01083 * fft99 AND fft991 01084 * TWO FFT ROUTINES THAT RETURN SLIGHTLY DIFFERENT 01085 * ARRANGEMENTS OF THE DATA IN GRIDPOINT SPACE. 01086 * 01087 * USAGE LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 1, 01088 * Q .GE. 0, AND R .GE. 0. THEN A TYPICAL SEQUENCE OF 01089 * CALLS TO TRANSFORM A GIVEN SET OF REAL VECTORS OF LENGTH 01090 * N TO A SET OF #HALF-COMPLEX# FOURIER COEFFICIENT VECTORS 01091 * OF LENGTH N IS 01092 * 01093 * long ifax[13]; 01094 * double trigs[3*n/2+1],a[m*(n+2)], work[m*(n+1)]; 01095 * 01096 * fftfax (n, ifax, trigs); 01097 * fft99 (a,work,trigs,ifax,inc,jump,n,m,isign); 01098 * 01099 * SEE THE INDIVIDUAL WRITE-UPS FOR FFTFAX, FFT99, AND 01100 * FFT991 BELOW, FOR A DETAILED DESCRIPTION OF THE 01101 * ARGUMENTS. 01102 * 01103 * HISTORY THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN 01104 * NOVEMBER, 1978. IT WAS MODIFIED, DOCUMENTED, AND TESTED 01105 * FOR NCAR BY RUSS REW IN SEPTEMBER, 1980. 01106 * 01107 * ----------------------------------------------------------------------- 01108 * 01109 * fftfax (n,ifax,trigs) 01110 * 01111 * PURPOSE A SET-UP ROUTINE FOR FFT99 AND FFT991. IT NEED ONLY BE 01112 * CALLED ONCE BEFORE A SEQUENCE OF CALLS TO THE FFT 01113 * ROUTINES (PROVIDED THAT N IS NOT CHANGED). 01114 * 01115 * ARGUMENT IFAX(13),TRIGS(3*N/2+1) 01116 * DIMENSIONS 01117 * 01118 * ARGUMENTS 01119 * 01120 * ON INPUT N 01121 * AN EVEN NUMBER GREATER THAN 4 THAT HAS NO PRIME FACTOR 01122 * GREATER THAN 5. N IS THE LENGTH OF THE TRANSFORMS (SEE 01123 * THE DOCUMENTATION FOR FFT99 AND FFT991 FOR THE 01124 * DEFINITIONS OF THE TRANSFORMS). 01125 * 01126 * IFAX 01127 * AN INTEGER ARRAY. THE NUMBER OF ELEMENTS ACTUALLY USED 01128 * WILL DEPEND ON THE FACTORIZATION OF N. DIMENSIONING 01129 * IFAX FOR 13 SUFFICES FOR ALL N LESS THAN A MILLION. 01130 * 01131 * TRIGS 01132 * A FLOATING POINT ARRAY OF DIMENSION 3*N/2 IF N/2 IS 01133 * EVEN, OR 3*N/2+1 IF N/2 IS ODD. 01134 * 01135 * ON OUTPUT IFAX 01136 * CONTAINS THE FACTORIZATION OF N/2. IFAX(1) IS THE 01137 * NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED 01138 * IN IFAX(2),IFAX(3),... IF FFTFAX IS CALLED WITH N ODD, 01139 * OR IF N HAS ANY PRIME FACTORS GREATER THAN 5, IFAX(1) 01140 * IS SET TO -99. 01141 * 01142 * TRIGS 01143 * AN ARRAY OF TRIGNOMENTRIC FUNCTION VALUES SUBSEQUENTLY 01144 * USED BY THE FFT ROUTINES. 01145 * 01146 * ----------------------------------------------------------------------- 01147 * 01148 * fft991 (a,work,trigs,ifax,inc,jump,n,m,isign) 01149 * AND 01150 * fft99 (a,work,trigs,ifax,inc,jump,n,m,isign) 01151 * 01152 * PURPOSE PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX 01153 * PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE 01154 * TRANSFORMS, USING ORDINARY SPATIAL ORDER OF GRIDPOINT 01155 * VALUES (FFT991) OR EXPLICIT CYCLIC CONTINUITY IN THE 01156 * GRIDPOINT VALUES (FFT99). GIVEN A SET 01157 * OF REAL DATA VECTORS, THE PACKAGE RETURNS A SET OF 01158 * #HALF-COMPLEX# FOURIER COEFFICIENT VECTORS, OR VICE 01159 * VERSA. THE LENGTH OF THE TRANSFORMS MUST BE AN EVEN 01160 * NUMBER THAT HAS NO OTHER FACTORS EXCEPT POSSIBLY POWERS 01161 * OF 2, 3, AND 5. THESE VERSION OF FFT991 AND FFT99 ARE 01162 * OPTIMIZED FOR USE ON THE CRAY-1. 01163 * 01164 * ARGUMENT A(M*(N+2)), WORK(M*(N+1)), TRIGS(3*N/2+1), IFAX(13) 01165 * DIMENSIONS 01166 * 01167 * ARGUMENTS 01168 * 01169 * ON INPUT A 01170 * AN ARRAY OF LENGTH M*(N+2) CONTAINING THE INPUT DATA 01171 * OR COEFFICIENT VECTORS. THIS ARRAY IS OVERWRITTEN BY 01172 * THE RESULTS. 01173 * 01174 * WORK 01175 * A WORK ARRAY OF DIMENSION M*(N+1) 01176 * 01177 * TRIGS 01178 * AN ARRAY SET UP BY FFTFAX, WHICH MUST BE CALLED FIRST. 01179 * 01180 * IFAX 01181 * AN ARRAY SET UP BY FFTFAX, WHICH MUST BE CALLED FIRST. 01182 * 01183 * INC 01184 * THE INCREMENT (IN WORDS) BETWEEN SUCCESSIVE ELEMENTS OF 01185 * EACH DATA OR COEFFICIENT VECTOR (E.G. INC=1 FOR 01186 * CONSECUTIVELY STORED DATA). 01187 * 01188 * JUMP 01189 * THE INCREMENT (IN WORDS) BETWEEN THE FIRST ELEMENTS OF 01190 * SUCCESSIVE DATA OR COEFFICIENT VECTORS. ON THE CRAY-1, 01191 * TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8 01192 * (TO AVOID MEMORY BANK CONFLICTS). FOR CLARIFICATION OF 01193 * INC AND JUMP, SEE THE EXAMPLES BELOW. 01194 * 01195 * N 01196 * THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF 01197 * TRANSFORMS, BELOW). 01198 * 01199 * M 01200 * THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY. 01201 * 01202 * ISIGN 01203 * = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO 01204 * GRIDPOINT VALUES. 01205 * = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER 01206 * COEFFICIENTS. 01207 * 01208 * ON OUTPUT A 01209 * IF ISIGN = +1, AND M COEFFICIENT VECTORS ARE SUPPLIED 01210 * EACH CONTAINING THE SEQUENCE: 01211 * 01212 * A(0),B(0),A(1),B(1),...,A(N/2),B(N/2) (N+2 VALUES) 01213 * 01214 * THEN THE RESULT CONSISTS OF M DATA VECTORS EACH 01215 * CONTAINING THE CORRESPONDING N+2 GRIDPOINT VALUES: 01216 * 01217 * FOR FFT991, X(0), X(1), X(2),...,X(N-1),0,0. 01218 * FOR FFT99, X(N-1),X(0),X(1),X(2),...,X(N-1),X(0). 01219 * (EXPLICIT CYCLIC CONTINUITY) 01220 * 01221 * WHEN ISIGN = +1, THE TRANSFORM IS DEFINED BY: 01222 * X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) 01223 * WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) 01224 * AND I=SQRT (-1) 01225 * 01226 * IF ISIGN = -1, AND M DATA VECTORS ARE SUPPLIED EACH 01227 * CONTAINING A SEQUENCE OF GRIDPOINT VALUES X(J) AS 01228 * DEFINED ABOVE, THEN THE RESULT CONSISTS OF M VECTORS 01229 * EACH CONTAINING THE CORRESPONDING FOURIER COFFICIENTS 01230 * A(K), B(K), 0 .LE. K .LE N/2. 01231 * 01232 * WHEN ISIGN = -1, THE INVERSE TRANSFORM IS DEFINED BY: 01233 * C(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*EXP(-2*I*J*K*PI/N)) 01234 * WHERE C(K)=A(K)+I*B(K) AND I=SQRT(-1) 01235 * 01236 * A CALL WITH ISIGN=+1 FOLLOWED BY A CALL WITH ISIGN=-1 01237 * (OR VICE VERSA) RETURNS THE ORIGINAL DATA. 01238 * 01239 * NOTE: THE FACT THAT THE GRIDPOINT VALUES X(J) ARE REAL 01240 * IMPLIES THAT B(0)=B(N/2)=0. FOR A CALL WITH ISIGN=+1, 01241 * IT IS NOT ACTUALLY NECESSARY TO SUPPLY THESE ZEROS. 01242 * 01243 * EXAMPLES GIVEN 19 DATA VECTORS EACH OF LENGTH 64 (+2 FOR EXPLICIT 01244 * CYCLIC CONTINUITY), COMPUTE THE CORRESPONDING VECTORS OF 01245 * FOURIER COEFFICIENTS. THE DATA MAY, FOR EXAMPLE, BE 01246 * ARRANGED LIKE THIS: 01247 * 01248 * FIRST DATA A(1)= . . . A(66)= A(70) 01249 * VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS) 01250 * 01251 * SECOND DATA A(71)= . . . A(140) 01252 * VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS) 01253 * 01254 * AND SO ON. HERE INC=1, JUMP=70, N=64, M=19, ISIGN=-1, 01255 * AND FFT99 SHOULD BE USED (BECAUSE OF THE EXPLICIT CYCLIC 01256 * CONTINUITY). 01257 * 01258 * ALTERNATIVELY THE DATA MAY BE ARRANGED LIKE THIS: 01259 * 01260 * FIRST SECOND LAST 01261 * DATA DATA DATA 01262 * VECTOR VECTOR VECTOR 01263 * 01264 * A(1)= A(2)= A(19)= 01265 * 01266 * X(63) X(63) . . . X(63) 01267 * A(20)= X(0) X(0) . . . X(0) 01268 * A(39)= X(1) X(1) . . . X(1) 01269 * . . . 01270 * . . . 01271 * . . . 01272 * 01273 * IN WHICH CASE WE HAVE INC=19, JUMP=1, AND THE REMAINING 01274 * PARAMETERS ARE THE SAME AS BEFORE. IN EITHER CASE, EACH 01275 * COEFFICIENT VECTOR OVERWRITES THE CORRESPONDING INPUT 01276 * DATA VECTOR. 01277 * 01278 * ----------------------------------------------------------------------- 01279 * 01280 * SUBROUTINE ^FFT99^ - MULTIPLE FAST REAL PERIODIC TRANSFORM 01281 * CORRESPONDING TO OLD SCALAR ROUTINE FFT9 01282 * PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM 01283 * IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 01284 * (1970), 315-337) 01285 * 01286 * A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA 01287 * WORK IS AN AREA OF SIZE (N+1)*LOT 01288 * TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES 01289 * IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 01290 * INC IS THE INCREMENT WITHIN EACH DATA #VECTOR# 01291 * (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) 01292 * JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR 01293 * N IS THE LENGTH OF THE DATA VECTORS 01294 * LOT IS THE NUMBER OF DATA VECTORS 01295 * ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT 01296 * = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL 01297 * 01298 * ORDERING OF COEFFICIENTS: 01299 * A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) 01300 * WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED 01301 * 01302 * ORDERING OF DATA: 01303 * X(N-1),X(0),X(1),X(2),...,X(N),X(0) 01304 * I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED 01305 * 01306 * VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN 01307 * PARALLEL 01308 * 01309 * *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER 01310 * 01311 * DEFINITION OF TRANSFORMS: 01312 * ------------------------- 01313 * 01314 * ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) 01315 * WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) 01316 * 01317 * ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) 01318 * B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) 01319 */ 01320 01321 /* Parameter adjustments */ 01322 --work; 01323 --a; 01324 01325 /* Function Body */ 01326 nfax = ifax[0]; 01327 nx = n + 1; 01328 nh = n / 2; 01329 ink = inc + inc; 01330 if (isign == 1) { 01331 goto L30; 01332 } 01333 01334 /* IF NECESSARY, TRANSFER DATA TO WORK AREA */ 01335 igo = 50; 01336 if (nfax % 2 == 1) { 01337 goto L40; 01338 } 01339 ibase = inc + 1; 01340 jbase = 1; 01341 for (l = 1; l <= lot; ++l) { 01342 i = ibase; 01343 j = jbase; 01344 /* DIR$ IVDEP */ 01345 for (m = 1; m <= n; ++m) { 01346 work[j] = a[i]; 01347 i += inc; 01348 ++j; 01349 /* L10: */ 01350 } 01351 ibase += jump; 01352 jbase += nx; 01353 /* L20: */ 01354 } 01355 01356 igo = 60; 01357 goto L40; 01358 01359 /* PREPROCESSING (ISIGN=+1) */ 01360 /* ------------------------ */ 01361 01362 L30: 01363 fft99a(&a[1], &work[1], trigs, inc, jump, n, lot); 01364 igo = 60; 01365 01366 /* COMPLEX TRANSFORM */ 01367 /* ----------------- */ 01368 01369 L40: 01370 ia = inc + 1; 01371 la = 1; 01372 for (k = 1; k <= nfax; ++k) { 01373 if (igo == 60) { 01374 goto L60; 01375 } 01376 /* L50: */ 01377 vpassm(&a[ia], &a[ia + inc], &work[1], &work[2], trigs, ink, 2, jump, nx, 01378 lot, nh, ifax[k], la); 01379 igo = 60; 01380 goto L70; 01381 L60: 01382 vpassm(&work[1], &work[2], &a[ia], &a[ia + inc], trigs, 2, ink, nx, jump, 01383 lot, nh, ifax[k], la); 01384 igo = 50; 01385 L70: 01386 la *= ifax[k]; 01387 /* L80: */ 01388 } 01389 01390 if (isign == -1) { 01391 goto L130; 01392 } 01393 01394 /* IF NECESSARY, TRANSFER DATA FROM WORK AREA */ 01395 if (nfax % 2 == 1) { 01396 goto L110; 01397 } 01398 ibase = 1; 01399 jbase = ia; 01400 for (l = 1; l <= lot; ++l) { 01401 i = ibase; 01402 j = jbase; 01403 /* DIR$ IVDEP */ 01404 for (m = 1; m <= n; ++m) { 01405 a[j] = work[i]; 01406 ++i; 01407 j += inc; 01408 /* L90: */ 01409 } 01410 ibase += nx; 01411 jbase += jump; 01412 /* L100: */ 01413 } 01414 01415 /* FILL IN CYCLIC BOUNDARY POINTS */ 01416 L110: 01417 ia = 1; 01418 ib = n * inc + 1; 01419 /* DIR$ IVDEP */ 01420 for (l = 1; l <= lot; ++l) { 01421 a[ia] = a[ib]; 01422 a[ib + inc] = a[ia + inc]; 01423 ia += jump; 01424 ib += jump; 01425 /* L120: */ 01426 } 01427 goto L140; 01428 01429 /* POSTPROCESSING (ISIGN=-1): */ 01430 /* -------------------------- */ 01431 01432 L130: 01433 fft99b(&work[1], &a[1], trigs, inc, jump, n, lot); 01434 01435 L140: 01436 return; 01437 } /* fft99 */ 01438 01439 static void fft99a(double *a, double *work, double *trigs, long inc, long jump, 01440 long n, long lot) 01441 { 01442 /* Local variables */ 01443 static double c; 01444 static long k, l; 01445 static double s; 01446 static long ia, ib, ja, jb, iabase, ibbase, nh, jabase, jbbase, nx, 01447 ink; 01448 01449 01450 /* SUBROUTINE FFT99A - PREPROCESSING STEP FOR FFT99, ISIGN=+1 */ 01451 /* (SPECTRAL TO GRIDPOINT TRANSFORM) */ 01452 01453 /* Parameter adjustments */ 01454 --work; 01455 --a; 01456 01457 /* Function Body */ 01458 nh = n / 2; 01459 nx = n + 1; 01460 ink = inc + inc; 01461 01462 /* A(0) AND A(N/2) */ 01463 ia = 1; 01464 ib = n * inc + 1; 01465 ja = 1; 01466 jb = 2; 01467 /* DIR$ IVDEP */ 01468 for (l = 1; l <= lot; ++l) { 01469 work[ja] = a[ia] + a[ib]; 01470 work[jb] = a[ia] - a[ib]; 01471 ia += jump; 01472 ib += jump; 01473 ja += nx; 01474 jb += nx; 01475 /* L10: */ 01476 } 01477 01478 /* REMAINING WAVENUMBERS */ 01479 iabase = (inc << 1) + 1; 01480 ibbase = (n - 2) * inc + 1; 01481 jabase = 3; 01482 jbbase = n - 1; 01483 01484 for (k = 3; k <= nh; k += 2) { 01485 ia = iabase; 01486 ib = ibbase; 01487 ja = jabase; 01488 jb = jbbase; 01489 c = trigs[n + k - 1]; 01490 s = trigs[n + k]; 01491 /* DIR$ IVDEP */ 01492 for (l = 1; l <= lot; ++l) { 01493 work[ja] = a[ia] + a[ib] - (s * (a[ia] - a[ib]) + 01494 c * (a[ia + inc] + a[ib + inc])); 01495 work[jb] = a[ia] + a[ib] + (s * (a[ia] - a[ib]) + 01496 c * (a[ia + inc] + a[ib + inc])); 01497 work[ja + 1] = c * (a[ia] - a[ib]) - 01498 s * (a[ia + inc] + a[ib + inc]) + (a[ia + inc] - a[ib + inc]); 01499 work[jb + 1] = c * (a[ia] - a[ib]) - 01500 s * (a[ia + inc] + a[ib + inc]) - (a[ia + inc] - a[ib + inc]); 01501 ia += jump; 01502 ib += jump; 01503 ja += nx; 01504 jb += nx; 01505 /* L20: */ 01506 } 01507 iabase += ink; 01508 ibbase -= ink; 01509 jabase += 2; 01510 jbbase += -2; 01511 /* L30: */ 01512 } 01513 01514 if (iabase != ibbase) { 01515 goto L50; 01516 } 01517 /* WAVENUMBER N/4 (IF IT EXISTS) */ 01518 ia = iabase; 01519 ja = jabase; 01520 /* DIR$ IVDEP */ 01521 for (l = 1; l <= lot; ++l) { 01522 work[ja] = a[ia] * 2.0; 01523 work[ja + 1] = a[ia + inc] * -2.0; 01524 ia += jump; 01525 ja += nx; 01526 /* L40: */ 01527 } 01528 01529 L50: 01530 return; 01531 } /* fft99a */ 01532 01533 static void fft99b(double *work, double *a, double *trigs, long inc, long jump, 01534 long n, long lot) 01535 { 01536 /* Local variables */ 01537 static double c; 01538 static long k, l; 01539 static double s, scale; 01540 static long ia, ib, ja, jb, iabase, ibbase, nh, jabase, jbbase, nx, 01541 ink; 01542 01543 01544 /* SUBROUTINE FFT99B - POSTPROCESSING STEP FOR FFT99, ISIGN=-1 */ 01545 /* (GRIDPOINT TO SPECTRAL TRANSFORM) */ 01546 01547 /* Parameter adjustments */ 01548 --a; 01549 --work; 01550 01551 /* Function Body */ 01552 nh = n / 2; 01553 nx = n + 1; 01554 ink = inc + inc; 01555 01556 /* A(0) AND A(N/2) */ 01557 scale = 1.0 / (double) n; 01558 ia = 1; 01559 ib = 2; 01560 ja = 1; 01561 jb = n * inc + 1; 01562 /* DIR$ IVDEP */ 01563 for (l = 1; l <= lot; ++l) { 01564 a[ja] = scale * (work[ia] + work[ib]); 01565 a[jb] = scale * (work[ia] - work[ib]); 01566 a[ja + inc] = 0.0; 01567 a[jb + inc] = 0.0; 01568 ia += nx; 01569 ib += nx; 01570 ja += jump; 01571 jb += jump; 01572 /* L10: */ 01573 } 01574 01575 /* REMAINING WAVENUMBERS */ 01576 scale *= 0.5; 01577 iabase = 3; 01578 ibbase = n - 1; 01579 jabase = (inc << 1) + 1; 01580 jbbase = (n - 2) * inc + 1; 01581 01582 for (k = 3; k <= nh; k += 2) { 01583 ia = iabase; 01584 ib = ibbase; 01585 ja = jabase; 01586 jb = jbbase; 01587 c = trigs[n + k - 1]; 01588 s = trigs[n + k]; 01589 /* DIR$ IVDEP */ 01590 for (l = 1; l <= lot; ++l) { 01591 a[ja] = scale * (work[ia] + work[ib] + 01592 (c * (work[ia + 1] + work[ib + 1]) + 01593 s * (work[ia] - work[ib]))); 01594 a[jb] = scale * (work[ia] + work[ib] - 01595 (c * (work[ia + 1] + work[ib + 1]) + 01596 s * (work[ia] - work[ib]))); 01597 a[ja + inc] = scale * (c * (work[ia] - work[ib]) - 01598 s * (work[ia + 1] + work[ib + 1]) + 01599 (work[ib + 1] - work[ia + 1])); 01600 a[jb + inc] = scale * (c * (work[ia] - work[ib]) - 01601 s * (work[ia + 1] + work[ib + 1]) - 01602 (work[ib + 1] - work[ia + 1])); 01603 ia += nx; 01604 ib += nx; 01605 ja += jump; 01606 jb += jump; 01607 /* L20: */ 01608 } 01609 iabase += 2; 01610 ibbase += -2; 01611 jabase += ink; 01612 jbbase -= ink; 01613 /* L30: */ 01614 } 01615 01616 if (iabase != ibbase) { 01617 goto L50; 01618 } 01619 /* WAVENUMBER N/4 (IF IT EXISTS) */ 01620 ia = iabase; 01621 ja = jabase; 01622 scale *= 2.0; 01623 /* DIR$ IVDEP */ 01624 for (l = 1; l <= lot; ++l) { 01625 a[ja] = scale * work[ia]; 01626 a[ja + inc] = -scale * work[ia + 1]; 01627 ia += nx; 01628 ja += jump; 01629 /* L40: */ 01630 } 01631 01632 L50: 01633 return; 01634 } /* fft99b */ 01635 01636 void fft991(double *a, double *work, double *trigs, long *ifax, long inc, 01637 long jump, long n, long lot, long isign) 01638 { 01639 /* Local variables */ 01640 static long nfax, i, j, k, l, m, ibase, jbase; 01641 static long ia, ib, la, nh, nx; 01642 static long igo, ink; 01643 01644 01645 /* SUBROUTINE ^FFT991^ - MULTIPLE REAL/HALF-COMPLEX PERIODIC */ 01646 /* FAST FOURIER TRANSFORM */ 01647 01648 /* SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO */ 01649 /* THAT IN MRFFT2 */ 01650 01651 /* PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM */ 01652 /* IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 */ 01653 /* (1970), 315-337) */ 01654 01655 /* A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA */ 01656 /* WORK IS AN AREA OF SIZE (N+1)*LOT */ 01657 /* TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES */ 01658 /* IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 */ 01659 /* INC IS THE INCREMENT WITHIN EACH DATA #VECTOR# */ 01660 /* (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) */ 01661 /* JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR */ 01662 /* N IS THE LENGTH OF THE DATA VECTORS */ 01663 /* LOT IS THE NUMBER OF DATA VECTORS */ 01664 /* ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT */ 01665 /* = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL */ 01666 01667 /* ORDERING OF COEFFICIENTS: */ 01668 /* A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) */ 01669 /* WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED */ 01670 01671 /* ORDERING OF DATA: */ 01672 /* X(0),X(1),X(2),...,X(N-1) */ 01673 01674 /* VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN */ 01675 /* PARALLEL */ 01676 01677 /* *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER */ 01678 01679 /* DEFINITION OF TRANSFORMS: */ 01680 /* ------------------------- */ 01681 01682 /* ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) */ 01683 /* WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) */ 01684 01685 /* ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) */ 01686 /* B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) */ 01687 01688 /* Parameter adjustments */ 01689 --work; 01690 --a; 01691 01692 /* Function Body */ 01693 nfax = ifax[0]; 01694 nx = n + 1; 01695 nh = n / 2; 01696 ink = inc + inc; 01697 if (isign == 1) { 01698 goto L30; 01699 } 01700 01701 /* IF NECESSARY, TRANSFER DATA TO WORK AREA */ 01702 igo = 50; 01703 if (nfax % 2 == 1) { 01704 goto L40; 01705 } 01706 ibase = 1; 01707 jbase = 1; 01708 for (l = 1; l <= lot; ++l) { 01709 i = ibase; 01710 j = jbase; 01711 /* DIR$ IVDEP */ 01712 for (m = 1; m <= n; ++m) { 01713 work[j] = a[i]; 01714 i += inc; 01715 ++j; 01716 /* L10: */ 01717 } 01718 ibase += jump; 01719 jbase += nx; 01720 /* L20: */ 01721 } 01722 01723 igo = 60; 01724 goto L40; 01725 01726 /* PREPROCESSING (ISIGN=+1) */ 01727 /* ------------------------ */ 01728 01729 L30: 01730 fft99a(&a[1], &work[1], trigs, inc, jump, n, lot); 01731 igo = 60; 01732 01733 /* COMPLEX TRANSFORM */ 01734 /* ----------------- */ 01735 01736 L40: 01737 ia = 1; 01738 la = 1; 01739 for (k = 1; k <= nfax; ++k) { 01740 if (igo == 60) { 01741 goto L60; 01742 } 01743 /* L50: */ 01744 vpassm(&a[ia], &a[ia + inc], &work[1], &work[2], trigs, ink, 2, jump, nx, 01745 lot, nh, ifax[k], la); 01746 igo = 60; 01747 goto L70; 01748 L60: 01749 vpassm(&work[1], &work[2], &a[ia], &a[ia + inc], trigs, 2, ink, nx, jump, 01750 lot, nh, ifax[k], la); 01751 igo = 50; 01752 L70: 01753 la *= ifax[k]; 01754 /* L80: */ 01755 } 01756 01757 if (isign == -1) { 01758 goto L130; 01759 } 01760 01761 /* IF NECESSARY, TRANSFER DATA FROM WORK AREA */ 01762 if (nfax % 2 == 1) { 01763 goto L110; 01764 } 01765 ibase = 1; 01766 jbase = 1; 01767 for (l = 1; l <= lot; ++l) { 01768 i = ibase; 01769 j = jbase; 01770 /* DIR$ IVDEP */ 01771 for (m = 1; m <= n; ++m) { 01772 a[j] = work[i]; 01773 ++i; 01774 j += inc; 01775 /* L90: */ 01776 } 01777 ibase += nx; 01778 jbase += jump; 01779 /* L100: */ 01780 } 01781 01782 /* FILL IN ZEROS AT END */ 01783 L110: 01784 ib = n * inc + 1; 01785 /* DIR$ IVDEP */ 01786 for (l = 1; l <= lot; ++l) { 01787 a[ib] = 0.0; 01788 a[ib + inc] = 0.0; 01789 ib += jump; 01790 /* L120: */ 01791 } 01792 goto L140; 01793 01794 /* POSTPROCESSING (ISIGN=-1): */ 01795 /* -------------------------- */ 01796 01797 L130: 01798 fft99b(&work[1], &a[1], trigs, inc, jump, n, lot); 01799 01800 L140: 01801 return; 01802 } /* fft991 */ 01803 01804 void fftfax(long n, long *ifax, double *trigs) 01805 { 01806 /* Local variables */ 01807 static long i; 01808 01809 /* Function Body */ 01810 fax(ifax, n); 01811 i = ifax[0]; 01812 if (ifax[i] > 5 || n <= 4) { 01813 ifax[0] = -99; 01814 } 01815 if (ifax[0] <= 0) 01816 return; 01817 01818 fftrig(trigs, n); 01819 return; 01820 } /* fftfax */ 01821 01822 static void fax(long *ifax, long n) 01823 { 01824 /* Local variables */ 01825 static long nfax, item, i, k, l, istop, ii, nn, inc; 01826 01827 /* Function Body */ 01828 nn = n / 2; 01829 if (nn + nn == n) { 01830 goto L10; 01831 } 01832 ifax[0] = -99; 01833 return; 01834 L10: 01835 k = 0; 01836 /* TEST FOR FACTORS OF 4 */ 01837 L20: 01838 if (nn % 4 != 0) { 01839 goto L30; 01840 } 01841 ++k; 01842 ifax[k] = 4; 01843 nn /= 4; 01844 if (nn == 1) { 01845 goto L80; 01846 } 01847 goto L20; 01848 /* TEST FOR EXTRA FACTOR OF 2 */ 01849 L30: 01850 if (nn % 2 != 0) { 01851 goto L40; 01852 } 01853 ++k; 01854 ifax[k] = 2; 01855 nn /= 2; 01856 if (nn == 1) { 01857 goto L80; 01858 } 01859 /* TEST FOR FACTORS OF 3 */ 01860 L40: 01861 if (nn % 3 != 0) { 01862 goto L50; 01863 } 01864 ++k; 01865 ifax[k] = 3; 01866 nn /= 3; 01867 if (nn == 1) { 01868 goto L80; 01869 } 01870 goto L40; 01871 /* NOW FIND REMAINING FACTORS */ 01872 L50: 01873 l = 5; 01874 inc = 2; 01875 /* INC ALTERNATELY TAKES ON VALUES 2 AND 4 */ 01876 L60: 01877 if (nn % l != 0) { 01878 goto L70; 01879 } 01880 ++k; 01881 ifax[k] = l; 01882 nn /= l; 01883 if (nn == 1) { 01884 goto L80; 01885 } 01886 goto L60; 01887 L70: 01888 l += inc; 01889 inc = 6 - inc; 01890 goto L60; 01891 L80: 01892 ifax[0] = k; 01893 /* IFAX(0) CONTAINS NUMBER OF FACTORS */ 01894 nfax = ifax[0]; 01895 /* SORT FACTORS INTO ASCENDING ORDER */ 01896 if (nfax == 1) 01897 return; 01898 01899 for (ii = 2; ii <= nfax; ++ii) { 01900 istop = nfax + 2 - ii; 01901 for (i = 2; i <= istop; ++i) { 01902 if (ifax[i] < ifax[i - 1]) { 01903 item = ifax[i - 1]; 01904 ifax[i - 1] = ifax[i]; 01905 ifax[i] = item; 01906 } 01907 } 01908 } 01909 return; 01910 } /* fax */ 01911 01912 static void fftrig(double *trigs, long n) 01913 { 01914 /* Local variables */ 01915 static long i, n2; 01916 static double angle; 01917 static long la, nh; 01918 static long nn; 01919 static double del; 01920 01921 /* Function Body */ 01922 nn = n / 2; 01923 del = (PI + PI) / (double) nn; 01924 n2 = nn + nn; 01925 for (i = 0; i < n2; i += 2) { 01926 angle = (double) i * 0.5 * del; 01927 trigs[i] = cos(angle); 01928 trigs[i + 1] = sin(angle); 01929 } 01930 del *= 0.5; 01931 nh = (nn + 1) / 2; 01932 n2 = nh + nh; 01933 la = nn + nn; 01934 for (i = 0; i < n2; i += 2) { 01935 angle = (double) i * 0.5 * del; 01936 trigs[la + i] = cos(angle); 01937 trigs[la + i + 1] = sin(angle); 01938 } 01939 01940 return; 01941 } /* fftrig */ 01942