00001 //# FFTPack.h: C++ wrapper functions for Fortran FFTPACK code 00002 //# Copyright (C) 1993,1994,1995,1997,1999,2000,2001 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 00027 //# $Id$ 00028 00029 #ifndef SCIMATH_FFTPACK_H 00030 #define SCIMATH_FFTPACK_H 00031 00032 #include <casacore/casa/aips.h> 00033 00034 //# The SGI compiler with -LANG:std has some trouble including both Complexfwd.h 00035 //# and Complex.h so we bypass the problem by include Complex.h only. 00036 #if defined(AIPS_USE_NEW_SGI) 00037 #include <casacore/casa/BasicSL/Complex.h> 00038 #else 00039 #include <casacore/casa/BasicSL/Complexfwd.h> 00040 #endif 00041 00042 namespace casacore { //# NAMESPACE CASACORE - BEGIN 00043 00044 // <summary>C++ interface to the Fortran FFTPACK library</summary> 00045 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="" demos=""> 00046 // </reviewed> 00047 // <synopsis> 00048 // The static functions in this class are C++ wrappers to the Fortran FFTPACK 00049 // library. This library contains functions that perform fast Fourier 00050 // transforms (FFT's) and related transforms. 00051 00052 // An additional purpose of these definitions is to overload the functions so 00053 // that C++ users can access the functions in either fftpak (single precision) 00054 // or dfftpack (double precision) with identical function names. 00055 00056 // These routines only do one-dimensional transforms with the first element of 00057 // the array being the "origin" of the transform. The <linkto 00058 // class="FFTServer">FFTServer</linkto> class uses some of these functions to 00059 // implement multi-dimensional transforms with the origin of the transform 00060 // either at the centre or the first element of the Array. 00061 00062 // You must initialise the work array <src>wsave</src> before using the forward 00063 // transform (function with a suffix of f) or the backward transform (with a 00064 // suffix of b). 00065 00066 // The transforms done by the functions in this class can be categorised as 00067 // follows: 00068 // <ul> 00069 // <li> Complex to Complex Transforms<br> 00070 // Done by the cttfi, cfftf & cfftb functions 00071 // <li> Real to Complex Transforms<br> 00072 // Done by the rffti, rfftf & rfftb functions. A simpler interface is 00073 // provided by the ezffti, ezfftf & ezfftb functions. The 'ez' functions 00074 // do not destroy the input array and provide the result in a slightly 00075 // less packed format. They are available in single precision only and 00076 // internally use the rfft functions. 00077 // <li> Sine Transforms<br> 00078 // Done by the sinti & sint functions. As the sine transform is its own 00079 // inverse there is no need for any distinction between forward and 00080 // backward transforms. 00081 // <li> Cosine Transforms<br> 00082 // Done by the costi & cost functions. As the cosine transform is its own 00083 // inverse there is no need for any distinction between forward and 00084 // backward transforms. 00085 // <li> Sine quarter wave Transforms<br> 00086 // Done by the sinqi, sinqf & sinqb functions. 00087 // <li> Cosine quarter wave Transforms<br> 00088 // Done by the cosqi, cosqf & cosqb functions. 00089 // </ul> 00090 00091 00092 // <note role=warning> 00093 // These functions assume that it is possible to convert between Casacore numeric 00094 // types and those used by Fortran. That it is possible to convert between 00095 // Float & float, Double & double and Int & int. 00096 // </note> 00097 00098 // <note role=warning> 00099 // These function also assume that a Complex array is stored as pairs of 00100 // floating point numbers, with no intervening gaps, and with the real 00101 // component first ie., <src>[re0,im0,re1,im1, ...]</src> so that the following 00102 // type casts work, 00103 // <srcblock> 00104 // Complex* complexPtr; 00105 // Float* floatPtr = (Float* ) complexPtr; 00106 // </srcblock> 00107 // and allow a Complex number to be accessed as a pair of real numbers. If this 00108 // assumption is bad then float Arrays will have to generated by copying the 00109 // complex ones. When compiled in debug mode mode the functions that require 00110 // this assumption will throw an exception (AipsError) if this assumption is 00111 // bad. Ultimately this assumption about Complex<->Float Array conversion 00112 // should be put somewhere central like Array2Math.cc. 00113 // </note> 00114 00115 // </synopsis> 00116 00117 00118 class FFTPack 00119 { 00120 public: 00121 // cffti initializes the array wsave which is used in both <src>cfftf</src> and 00122 // <src>cfftb</src>. The prime factorization of n together with a tabulation of 00123 // the trigonometric functions are computed and stored in wsave. 00124 // 00125 // Input parameter: 00126 // <dl compact> 00127 // <dt><b>n</b> 00128 // <dd> The length of the sequence to be transformed 00129 // </dl> 00130 // Output parameter: 00131 // <dl compact> 00132 // <dt><b>wsave</b> 00133 // <dd> A work array which must be dimensioned at least 4*n+15 00134 // The same work array can be used for both cfftf and cfftb 00135 // as long as n remains unchanged. Different wsave arrays 00136 // are required for different values of n. The contents of 00137 // wsave must not be changed between calls of cfftf or cfftb. 00138 // </dl> 00139 // <group> 00140 static void cffti(Int n, Float* wsave); 00141 static void cffti(Int n, Double* wsave); 00142 //Here is the doc from FFTPack 5.1 00143 //You can convert the linguo from fortran to C/C++ 00144 /* Input Arguments 00145 00146 L Integer number of elements to be transformed in the first 00147 dimension. The transform is most efficient when L is a 00148 product of small primes. 00149 00150 00151 M Integer number of elements to be transformed in the second 00152 dimension. The transform is most efficient when M is a 00153 product of small primes. 00154 00155 LENSAV Integer dimension of WSAVE array. LENSAV must be at least 00156 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. 00157 00158 00159 Output Arguments 00160 00161 WSAVE Real work array with dimension LENSAV, containing the 00162 prime factors of L and M, and also containing certain 00163 trigonometric values which will be used in routines 00164 CFFT2B or CFFT2F. 00165 00166 00167 WSAVE Real work array with dimension LENSAV. The WSAVE array 00168 must be initialized with a call to subroutine CFFT2I before 00169 the first call to CFFT2B or CFFT2F, and thereafter whenever 00170 the values of L, M or the contents of array WSAVE change. 00171 Using different WSAVE arrays for different transform lengths 00172 or types in the same program may reduce computation costs 00173 because the array contents can be re-used. 00174 00175 00176 IER Integer error return 00177 = 0 successful exit 00178 = 2 input parameter LENSAV not big enough 00179 = 20 input error returned by lower level routine 00180 */ 00181 00182 static void cfft2i(const Int& n, const Int& m, Float *& wsave, const Int& lensav, Int& ier); 00183 // </group> 00184 00185 // cfftf computes the forward complex discrete Fourier 00186 // transform (the Fourier analysis). Equivalently, cfftf computes 00187 // the Fourier coefficients of a complex periodic sequence. 00188 // the transform is defined below at output parameter c. 00189 // 00190 // The transform is not normalized. To obtain a normalized transform 00191 // the output must be divided by n. Otherwise a call of cfftf 00192 // followed by a call of cfftb will multiply the sequence by n. 00193 // 00194 // The array wsave which is used by <src>cfftf</src> must be 00195 // initialized by calling <src>cffti(n,wsave)</src>. 00196 // 00197 // Input parameters: 00198 // <dl compact> 00199 // <dt><b>n</b> 00200 // <dd> The length of the complex sequence c. The method is 00201 // more efficient when n is the product of small primes. 00202 // <dt><b>c</b> 00203 // <dd> A complex array of length n which contains the sequence to be 00204 // transformed. 00205 // <dt><b>wsave</b> 00206 // <dd> A real work array which must be dimensioned at least 4n+15 00207 // by the program that calls cfftf. The wsave array must be 00208 // initialized by calling <src>cffti(n,wsave)</src> and a 00209 // different wsave array must be used for each different 00210 // value of n. This initialization does not have to be 00211 // repeated so long as n remains unchanged thus subsequent 00212 // transforms can be obtained faster than the first. 00213 // The same wsave array can be used by cfftf and cfftb. 00214 // </dl> 00215 // Output parameters: 00216 // <dl compact> 00217 // <dt><b>c</b> 00218 // <dd> for j=1,...,n<br> 00219 // c(j)=the sum from k=1,...,n of<br> 00220 // c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)<br> 00221 // where i=sqrt(-1)<br> 00222 // <dt><b>wsave</b> 00223 // <dd> Contains initialization calculations which must not be 00224 // destroyed between calls of cfftf or cfftb 00225 // </dl> 00226 // <group> 00227 static void cfftf(Int n, Complex* c, Float* wsave); 00228 static void cfftf(Int n, DComplex* c, Double* wsave); 00229 00230 //Description from FFTPack 5.1 00231 /* 00232 Input Arguments 00233 00234 00235 LDIM Integer first dimension of two-dimensional complex array C. 00236 00237 00238 L Integer number of elements to be transformed in the first 00239 dimension of the two-dimensional complex array C. The value 00240 of L must be less than or equal to that of LDIM. The 00241 transform is most efficient when L is a product of small 00242 primes. 00243 00244 00245 M Integer number of elements to be transformed in the second 00246 dimension of the two-dimensional complex array C. The 00247 transform is most efficient when M is a product of small 00248 primes. 00249 00250 C Complex array of two dimensions containing the (L,M) subarray 00251 to be transformed. C's first dimension is LDIM, its second 00252 dimension must be at least M. 00253 00254 WSAVE Real work array with dimension LENSAV. WSAVE's contents 00255 must be initialized with a call to subroutine CFFT2I before 00256 the first call to routine CFFT2F or CFFT2B with transform 00257 lengths L and M. WSAVE's contents may be re-used for 00258 subsequent calls to CFFT2F and CFFT2B having those same 00259 transform lengths. 00260 00261 00262 00263 LENSAV Integer dimension of WSAVE array. LENSAV must be at least 00264 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. 00265 00266 00267 WORK Real work array. 00268 00269 00270 LENWRK Integer dimension of WORK array. LENWRK must be at least 00271 2*L*M. 00272 */ 00273 static void cfft2f (const Int& ldim, const Int& L, const Int& M, Complex*& C, Float*& WSAVE, const Int& LENSAV, 00274 Float *& WORK, const Int& LENWRK, Int& IER); 00275 // </group> 00276 00277 // cfftb computes the backward complex discrete Fourier 00278 // transform (the Fourier synthesis). Equivalently, cfftb computes 00279 // a complex periodic sequence from its Fourier coefficients. 00280 // The transform is defined below with output parameter c. 00281 // 00282 // A call of cfftf followed by a call of cfftb will multiply the 00283 // sequence by n. 00284 // 00285 // The array wsave which is used by <src>cfftb</src> must be 00286 // initialized by calling <src>cffti(n,wsave)</src>. 00287 // 00288 // Input parameters: 00289 // <dl compact> 00290 // <dt><b>n</b> 00291 // <dd> The length of the complex sequence c. The method is 00292 // more efficient when n is the product of small primes. 00293 // <dt><b>c</b> 00294 // <dd> A complex array of length n which contains the sequence to be 00295 // transformed. 00296 // <dt><b>wsave</b> 00297 // <dd> A real work array which must be dimensioned at least 4n+15 00298 // in the program that calls cfftb. The wsave array must be 00299 // initialized by calling <src>cffti(n,wsave)</src> 00300 // and a different wsave array must be used for each different 00301 // value of n. This initialization does not have to be 00302 // repeated so long as n remains unchanged thus subsequent 00303 // transforms can be obtained faster than the first. 00304 // The same wsave array can be used by cfftf and cfftb. 00305 // </dl> 00306 // Output parameters: 00307 // <dl compact> 00308 // <dt><b>c</b> 00309 // <dd> for j=1,...,n<br> 00310 // c(j)=the sum from k=1,...,n of<br> 00311 // c(k)*exp(i*(j-1)*(k-1)*2*pi/n)<br> 00312 00313 // <dt><b>wsave</b> 00314 // <dd> Contains initialization calculations which must not be 00315 // destroyed between calls of cfftf or cfftb 00316 // </dl> 00317 // <group> 00318 static void cfftb(Int n, Complex* c, Float* wsave); 00319 static void cfftb(Int n, DComplex* c, Double* wsave); 00320 00321 00322 //Documentation from FFTPack 5.1 00323 /* 00324 Input Arguments 00325 00326 LDIM Integer first dimension of two-dimensional complex array C. 00327 00328 00329 L Integer number of elements to be transformed in the first 00330 dimension of the two-dimensional complex array C. The value 00331 of L must be less than or equal to that of LDIM. The 00332 transform is most efficient when L is a product of small 00333 primes. 00334 00335 00336 M Integer number of elements to be transformed in the second 00337 dimension of the two-dimensional complex array C. The 00338 transform is most efficient when M is a product of small 00339 primes. 00340 00341 C Complex array of two dimensions containing the (L,M) subarray 00342 to be transformed. C's first dimension is LDIM, its second 00343 dimension must be at least M. 00344 00345 00346 WSAVE Real work array with dimension LENSAV. WSAVE's contents 00347 must be initialized with a call to subroutine CFFT2I before 00348 the first call to routine CFFT2F or CFFT2B with transform 00349 lengths L and M. WSAVE's contents may be re-used for 00350 subsequent calls to CFFT2F and CFFT2B with the same 00351 transform lengths L and M. 00352 00353 LENSAV Integer dimension of WSAVE array. LENSAV must be at least 00354 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. 00355 00356 WORK Real work array. 00357 00358 LENWRK Integer dimension of WORK array. LENWRK must be at least 00359 2*L*M. 00360 00361 Output Arguments 00362 00363 00364 C Complex output array. For purposes of exposition, 00365 assume the index ranges of array C are defined by 00366 C(0:L-1,0:M-1). 00367 00368 For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given 00369 in the traditional aliased form by 00370 00371 L-1 M-1 00372 C(I,J) = SUM SUM C(L1,M1)* 00373 L1=0 M1=0 00374 00375 EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) 00376 00377 And in unaliased form, the C(I,J)'s are given by 00378 00379 LF MF 00380 C(I,J) = SUM SUM C(L1,M1,K1)* 00381 L1=LS M1=MS 00382 00383 EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) 00384 00385 where 00386 00387 LS= -L/2 and LF=L/2-1 if L is even; 00388 LS=-(L-1)/2 and LF=(L-1)/2 if L is odd; 00389 MS= -M/2 and MF=M/2-1 if M is even; 00390 MS=-(M-1)/2 and MF=(M-1)/2 if M is odd; 00391 00392 and 00393 00394 C(L1,M1) = C(L1+L,M1) if L1 is zero or negative; 00395 C(L1,M1) = C(L1,M1+M) if M1 is zero or negative; 00396 00397 The two forms give different results when used to 00398 interpolate between elements of the sequence. 00399 00400 IER Integer error return 00401 = 0 successful exit 00402 = 2 input parameter LENSAV not big enough 00403 = 3 input parameter LENWRK not big enough 00404 = 5 input parameter L > LDIM 00405 = 20 input error returned by lower level routine 00406 */ 00407 00408 00409 00410 static void cfft2b (const Int& LDIM, const Int& L, const Int& M, Complex* & C, Float *& WSAVE, const Int& LENSAV, 00411 Float*& WORK, const Int& LENWRK, Int& IER); 00412 // </group> 00413 00414 // rffti initializes the array wsave which is used in both <src>rfftf</src> and 00415 // <src>rfftb</src>. The prime factorization of n together with a tabulation of 00416 // the trigonometric functions are computed and stored in wsave. 00417 // 00418 // Input parameter: 00419 // <dl compact> 00420 // <dt><b>n</b> 00421 // <dd> The length of the sequence to be transformed. 00422 // </dl> 00423 // Output parameter: 00424 // <dl compact> 00425 // <dt><b>wsave</b> 00426 // <dd> A work array which must be dimensioned at least 2*n+15. 00427 // The same work array can be used for both rfftf and rfftb 00428 // as long as n remains unchanged. Different wsave arrays 00429 // are required for different values of n. The contents of 00430 // wsave must not be changed between calls of rfftf or rfftb. 00431 // </dl> 00432 // <group> 00433 static void rffti(Int n, Float* wsave); 00434 static void rffti(Int n, Double* wsave); 00435 // </group> 00436 00437 // rfftf computes the Fourier coefficients of a real perodic sequence (Fourier 00438 // analysis). The transform is defined below at output parameter r. 00439 // 00440 // Input parameters: 00441 // <dl compact> 00442 // <dt><b>n</b> 00443 // <dd> The length of the array r to be transformed. The method 00444 // is most efficient when n is a product of small primes. 00445 // n may change so long as different work arrays are provided 00446 // <dt><b>r</b> 00447 // <dd> A real array of length n which contains the sequence 00448 // to be transformed 00449 // <dt><b>wsave</b> 00450 // <dd> A work array which must be dimensioned at least 2*n+15 00451 // in the program that calls rfftf. The wsave array must be 00452 // initialized by calling <src>rffti(n,wsave)</src> and a 00453 // different wsave array must be used for each different 00454 // value of n. This initialization does not have to be 00455 // repeated so long as n remains unchanged thus subsequent 00456 // transforms can be obtained faster than the first. 00457 // The same wsave array can be used by rfftf and rfftb. 00458 // </dl> 00459 // output parameters 00460 // <dl compact> 00461 // <dt><b>r</b> 00462 // <dd> r(1) = the sum from i=1 to i=n of r(i)<br> 00463 // if n is even set l = n/2 , if n is odd set l = (n+1)/2<br> 00464 // then for k = 2,...,l<br> 00465 // r(2*k-2) = the sum from i = 1 to i = n of<br> 00466 // r(i)*cos((k-1)*(i-1)*2*pi/n)<br> 00467 // r(2*k-1) = the sum from i = 1 to i = n of<br> 00468 // -r(i)*sin((k-1)*(i-1)*2*pi/n)<br> 00469 // if n is even<br> 00470 // r(n) = the sum from i = 1 to i = n of<br> 00471 // (-1)**(i-1)*r(i)<br> 00472 // 00473 // note: 00474 // this transform is unnormalized since a call of rfftf 00475 // followed by a call of rfftb will multiply the input 00476 // sequence by n. 00477 // <dt><b>wsave</b> 00478 // <dd> Contains results which must not be destroyed between 00479 // calls of rfftf or rfftb. 00480 // </dl> 00481 // <group> 00482 static void rfftf(Int n, Float* r, Float* wsave); 00483 static void rfftf(Int n, Double* r, Double* wsave); 00484 // </group> 00485 00486 00487 // rfftb computes the real perodic sequence from its Fourier coefficients 00488 // (Fourier synthesis). The transform is defined below at output parameter r. 00489 // 00490 // Input parameters: 00491 // <dl compact> 00492 // <dt><b>n</b> 00493 // <dd> The length of the array r to be transformed. The method 00494 // is most efficient when n is a product of small primes. 00495 // n may change so long as different work arrays are provided 00496 // <dt><b>r</b> 00497 // <dd> A real array of length n which contains the sequence 00498 // to be transformed 00499 // <dt><b>wsave</b> 00500 // <dd> A work array which must be dimensioned at least 2*n+15 00501 // in the program that calls rfftb. The wsave array must be 00502 // initialized by calling <src>rffti(n,wsave)</src> and a 00503 // different wsave array must be used for each different 00504 // value of n. This initialization does not have to be 00505 // repeated so long as n remains unchanged thus subsequent 00506 // transforms can be obtained faster than the first. 00507 // The same wsave array can be used by rfftf and rfftb. 00508 // </dl> 00509 // Output parameters: 00510 // <dl compact> 00511 // <dt><b>r</b> 00512 // <dd> for n even and for i = 1,...,n<br> 00513 // r(i) = r(1)+(-1)**(i-1)*r(n)<br> 00514 // plus the sum from k=2 to k=n/2 of<br> 00515 // 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)<br> 00516 // -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)<br> 00517 // for n odd and for i = 1,...,n<br> 00518 // r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of<br> 00519 // 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)<br> 00520 // -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)<br> 00521 // 00522 // note: 00523 // this transform is unnormalized since a call of rfftf 00524 // followed by a call of rfftb will multiply the input 00525 // sequence by n. 00526 // <dt><b>wsave</b> 00527 // <dd> Contains results which must not be destroyed between 00528 // calls of rfftb or rfftf. 00529 // </dl> 00530 // <group> 00531 static void rfftb(Int n, Float* r, Float* wsave); 00532 static void rfftb(Int n, Double* r, Double* wsave); 00533 // </group> 00534 00535 // ezffti initializes the array wsave which is used in both <src>ezfftf</src> 00536 // and <src>ezfftb</src>. The prime factorization of n together with a 00537 // tabulation of the trigonometric functions are computed and stored in wsave. 00538 // 00539 // Input parameter: 00540 // <dl compact> 00541 // <dt><b>n</b> 00542 // <dd> The length of the sequence to be transformed. 00543 // </dl> 00544 // Output parameter: 00545 // <dl compact> 00546 // <dt><b>wsave</b> 00547 // <dd> A work array which must be dimensioned at least 3*n+15. 00548 // The same work array can be used for both ezfftf and ezfftb 00549 // as long as n remains unchanged. Different wsave arrays 00550 // are required for different values of n. 00551 // </dl> 00552 static void ezffti(Int n, Float* wsave); 00553 00554 // ezfftf computes the Fourier coefficients of a real 00555 // perodic sequence (Fourier analysis). The transform is defined 00556 // below at output parameters azero, a and b. ezfftf is a simplified 00557 // but slower version of rfftf. 00558 // 00559 // Input parameters: 00560 // <dl compact> 00561 // <dt><b>n</b> 00562 // <dd> The length of the array r to be transformed. The method 00563 // is most efficient when n is the product of small primes. 00564 // <dt><b>r</b> 00565 // <dd> A real array of length n which contains the sequence 00566 // to be transformed. r is not destroyed. 00567 // <dt><b>wsave</b> 00568 // <dd> A work array which must be dimensioned at least 3*n+15 00569 // in the program that calls ezfftf. The wsave array must be 00570 // initialized by calling <src>ezffti(n,wsave)</src> and a 00571 // different wsave array must be used for each different 00572 // value of n. This initialization does not have to be 00573 // repeated so long as n remains unchanged thus subsequent 00574 // transforms can be obtained faster than the first. 00575 // The same wsave array can be used by ezfftf and ezfftb. 00576 // </dl> 00577 // Output parameters: 00578 // <dl compact> 00579 // <dt><b>azero</b> 00580 // <dd> The sum from i=1 to i=n of r(i)/n 00581 // <dt><b>a,b</b> 00582 // <dd> Real arrays of length n/2 (n even) or (n-1)/2 (n odd)<br> 00583 // for n even<br> 00584 // b(n/2)=0, and <br> 00585 // a(n/2) is the sum from i=1 to i=n of (-1)**(i-1)*r(i)/n<br> 00586 // 00587 // for n even define kmax=n/2-1<br> 00588 // for n odd define kmax=(n-1)/2<br> 00589 // then for k=1,...,kmax<br> 00590 // a(k) equals the sum from i=1 to i=n of<br> 00591 // 2./n*r(i)*cos(k*(i-1)*2*pi/n)<br> 00592 // b(k) equals the sum from i=1 to i=n of<br> 00593 // 2./n*r(i)*sin(k*(i-1)*2*pi/n)<br> 00594 // </dl> 00595 static void ezfftf(Int n, Float* r, Float* azero, Float* a, Float* b, 00596 Float* wsave); 00597 00598 // ezfftb computes a real perodic sequence from its 00599 // Fourier coefficients (Fourier synthesis). The transform is 00600 // defined below at output parameter r. ezfftb is a simplified 00601 // but slower version of rfftb. 00602 // 00603 // Input parameters: 00604 // <dl compact> 00605 // <dt><b>n</b> 00606 // <dd> The length of the output array r. The method is most 00607 // efficient when n is the product of small primes. 00608 // <dt><b>azero</b> 00609 // <dd> The constant Fourier coefficient 00610 // <dt><b>a,b</b> 00611 // <dd> Arrays which contain the remaining Fourier coefficients 00612 // these arrays are not destroyed. 00613 // The length of these arrays depends on whether n is even or 00614 // odd. 00615 // If n is even n/2 locations are required, 00616 // if n is odd (n-1)/2 locations are required. 00617 // <dt><b>wsave</b> 00618 // <dd> A work array which must be dimensioned at least 3*n+15. 00619 // in the program that calls ezfftb. The wsave array must be 00620 // initialized by calling <src>ezffti(n,wsave)</src> and a 00621 // different wsave array must be used for each different 00622 // value of n. This initialization does not have to be 00623 // repeated so long as n remains unchanged thus subsequent 00624 // transforms can be obtained faster than the first. 00625 // The same wsave array can be used by ezfftf and ezfftb. 00626 // </dl> 00627 // Output parameters: 00628 // <dl compact> 00629 // <dt><b>r</b> 00630 // <dd> if n is even define kmax=n/2<br> 00631 // if n is odd define kmax=(n-1)/2<br> 00632 // then for i=1,...,n<br> 00633 // r(i)=azero plus the sum from k=1 to k=kmax of<br> 00634 // a(k)*cos(k*(i-1)*2*pi/n)+b(k)*sin(k*(i-1)*2*pi/n)<br> 00635 // where<br> 00636 // c(k) = .5*cmplx(a(k),-b(k)) for k=1,...,kmax<br> 00637 // c(-k) = conjg(c(k))<br> 00638 // c(0) = azero<br> 00639 // and i=sqrt(-1)<br> 00640 // </dl> 00641 static void ezfftb(Int n, Float* r, Float* azero, Float* a, Float* b, 00642 Float* wsave); 00643 00644 // sinti initializes the array wsave which is used in 00645 // <src>sint</src>. The prime factorization of n together with a tabulation of 00646 // the trigonometric functions are computed and stored in wsave. 00647 // 00648 // Input parameter: 00649 // <dl compact> 00650 // <dt><b>n</b> 00651 // <dd> The length of the sequence to be transformed. the method 00652 // is most efficient when n+1 is a product of small primes. 00653 // </dl> 00654 // Output parameter: 00655 // <dl compact> 00656 // <dt><b>wsave</b> 00657 // <dd> A work array with at least int(2.5*n+15) locations. 00658 // Different wsave arrays are required for different values 00659 // of n. The contents of wsave must not be changed between 00660 // calls of sint. 00661 // </dl> 00662 // <group> 00663 static void sinti(Int n, Float* wsave); 00664 static void sinti(Int n, Double* wsave); 00665 // </group> 00666 00667 // sint computes the discrete Fourier sine transform 00668 // of an odd sequence x(i). The transform is defined below at 00669 // output parameter x. 00670 // sint is the unnormalized inverse of itself since a call of sint 00671 // followed by another call of sint will multiply the input sequence 00672 // x by 2*(n+1). 00673 // The array wsave which is used by sint must be 00674 // initialized by calling <src>sinti(n,wsave)</src>. 00675 // 00676 // Input parameters: 00677 // <dl compact> 00678 // <dt><b>n</b> 00679 // <dd> The length of the sequence to be transformed. The method 00680 // is most efficient when n+1 is the product of small primes. 00681 // <dt><b>x</b> 00682 // <dd> An array which contains the sequence to be transformed 00683 // <dt><b>wsave</b> 00684 // <dd> A work array with dimension at least int(2.5*n+15) 00685 // in the program that calls sint. The wsave array must be 00686 // initialized by calling <src>sinti(n,wsave)</src> and a 00687 // different wsave array must be used for each different 00688 // value of n. This initialization does not have to be 00689 // repeated so long as n remains unchanged thus subsequent 00690 // transforms can be obtained faster than the first. 00691 // </dl> 00692 // Output parameters: 00693 // <dl compact> 00694 // <dt><b>x</b> 00695 // <dd> for i=1,...,n<br> 00696 // x(i) = the sum from k=1 to k=n<br> 00697 // 2*x(k)*sin(k*i*pi/(n+1))<br> 00698 // 00699 // a call of sint followed by another call of 00700 // sint will multiply the sequence x by 2*(n+1). 00701 // Hence sint is the unnormalized inverse 00702 // of itself. 00703 // 00704 // <dt><b>wsave</b> 00705 // <dd> Contains initialization calculations which must not be 00706 // destroyed between calls of sint. 00707 // </dl> 00708 // <group> 00709 static void sint(Int n, Float* x, Float* wsave); 00710 static void sint(Int n, Double* x, Double* wsave); 00711 // </group> 00712 00713 // costi initializes the array wsave which is used in 00714 // <src>cost</src>. The prime factorization of n together with a tabulation of 00715 // the trigonometric functions are computed and stored in wsave. 00716 // 00717 // Input parameter: 00718 // <dl compact> 00719 // <dt><b>n</b> 00720 // <dd> The length of the sequence to be transformed. The method 00721 // is most efficient when n-1 is a product of small primes. 00722 // </dl> 00723 // Output parameter: 00724 // <dl compact> 00725 // <dt><b>wsave</b> 00726 // <dd> A work array which must be dimensioned at least 3*n+15. 00727 // Different wsave arrays are required for different values 00728 // of n. The contents of wsave must not be changed between 00729 // calls of cost. 00730 // </dl> 00731 // <group> 00732 static void costi(Int n, Float* wsave); 00733 static void costi(Int n, Double* wsave); 00734 // </group> 00735 00736 // cost computes the discrete Fourier cosine transform 00737 // of an even sequence x(i). The transform is defined below at output 00738 // parameter x. 00739 // cost is the unnormalized inverse of itself since a call of cost 00740 // followed by another call of cost will multiply the input sequence 00741 // x by 2*(n-1). The transform is defined below at output parameter x. 00742 // The array wsave which is used by <src>cost</src> must be 00743 // initialized by calling <src>costi(n,wsave)</src>. 00744 // 00745 // Input parameters: 00746 // <dl compact> 00747 // <dt><b>n</b> 00748 // <dd> The length of the sequence x. n must be greater than 1. 00749 // The method is most efficient when n-1 is a product of 00750 // small primes. 00751 // <dt><b>x</b> 00752 // <dd> An array which contains the sequence to be transformed 00753 // <dt><b>wsave</b> 00754 // <dd> A work array which must be dimensioned at least 3*n+15 00755 // in the program that calls cost. The wsave array must be 00756 // initialized by calling <src>costi(n,wsave)</src> and a 00757 // different wsave array must be used for each different 00758 // value of n. This initialization does not have to be 00759 // repeated so long as n remains unchanged thus subsequent 00760 // transforms can be obtained faster than the first. 00761 // </dl> 00762 // Output parameters: 00763 // <dl compact> 00764 // <dt><b>x</b> 00765 // <dd> for i=1,...,n<br> 00766 // x(i) = x(1)+(-1)**(i-1)*x(n)<br> 00767 // + the sum from k=2 to k=n-1<br> 00768 // 2*x(k)*cos((k-1)*(i-1)*pi/(n-1))<br> 00769 // 00770 // a call of cost followed by another call of 00771 // cost will multiply the sequence x by 2*(n-1) 00772 // hence cost is the unnormalized inverse 00773 // of itself. 00774 // <dt><b>wsave</b> 00775 // <dd> Contains initialization calculations which must not be 00776 // destroyed between calls of cost. 00777 // </dl> 00778 // <group> 00779 static void cost(Int n, Float* x, Float* wsave); 00780 static void cost(Int n, Double* x, Double* wsave); 00781 // </group> 00782 00783 // sinqi initializes the array wsave which is used in both <src>sinqf</src> and 00784 // <src>sinqb</src>. The prime factorization of n together with a tabulation of 00785 // the trigonometric functions are computed and stored in wsave. 00786 // 00787 // Input parameter: 00788 // <dl compact> 00789 // <dt><b>n</b> 00790 // <dd> The length of the sequence to be transformed. The method 00791 // is most efficient when n is a product of small primes. 00792 // </dl> 00793 // Output parameter: 00794 // <dl compact> 00795 // <dt><b>wsave</b> 00796 // <dd> A work array which must be dimensioned at least 3*n+15. 00797 // The same work array can be used for both sinqf and sinqb 00798 // as long as n remains unchanged. Different wsave arrays 00799 // are required for different values of n. The contents of 00800 // wsave must not be changed between calls of sinqf or sinqb. 00801 // </dl> 00802 // <group> 00803 static void sinqi(Int n, Float* wsave); 00804 static void sinqi(Int n, Double* wsave); 00805 // </group> 00806 00807 // sinqf computes the fast Fourier transform of quarter wave data. That is, 00808 // sinqf computes the coefficients in a sine series representation with only 00809 // odd wave numbers. The transform is defined below at output parameter x. 00810 // 00811 // sinqb is the unnormalized inverse of sinqf since a call of sinqf followed by 00812 // a call of sinqb will multiply the input sequence x by 4*n. 00813 // 00814 // The array wsave which is used by sinqf must be initialized by calling 00815 // <src>sinqi(n,wsave)</src>. 00816 // 00817 // Input parameters: 00818 // <dl compact> 00819 // <dt><b>n</b> 00820 // <dd> The length of the array x to be transformed. The method 00821 // is most efficient when n is a product of small primes. 00822 // <dt><b>x</b> 00823 // <dd> An array which contains the sequence to be transformed 00824 // <dt><b>wsave</b> 00825 // A work array which must be dimensioned at least 3*n+15. 00826 // in the program that calls sinqf. The wsave array must be 00827 // initialized by calling <src>sinqi(n,wsave)</src> and a 00828 // different wsave array must be used for each different 00829 // value of n. This initialization does not have to be 00830 // repeated so long as n remains unchanged thus subsequent 00831 // transforms can be obtained faster than the first. 00832 // </dl> 00833 // Output parameters: 00834 // <dl compact> 00835 // <dt><b>x</b> 00836 // <dd> for i=1,...,n<br> 00837 // x(i) = (-1)**(i-1)*x(n)<br> 00838 // + the sum from k=1 to k=n-1 of<br> 00839 // 2*x(k)*sin((2*i-1)*k*pi/(2*n))<br> 00840 // 00841 // a call of sinqf followed by a call of 00842 // sinqb will multiply the sequence x by 4*n. 00843 // therefore sinqb is the unnormalized inverse 00844 // of sinqf. 00845 // <dt><b>wsave </b> 00846 // <dd> Contains initialization calculations which must not 00847 // be destroyed between calls of sinqf or sinqb. 00848 // </dl> 00849 // <group> 00850 static void sinqf(Int n, Float* x, Float* wsave); 00851 static void sinqf(Int n, Double* x, Double* wsave); 00852 // </group> 00853 00854 // sinqb computes the fast Fourier transform of quarter 00855 // wave data. that is, sinqb computes a sequence from its 00856 // representation in terms of a sine series with odd wave numbers. 00857 // the transform is defined below at output parameter x. 00858 // 00859 // sinqf is the unnormalized inverse of sinqb since a call of sinqb 00860 // followed by a call of sinqf will multiply the input sequence x 00861 // by 4*n. 00862 // 00863 // The array wsave which is used by <src>sinqb</src> must be 00864 // initialized by calling <src>sinqi(n,wsave)</src>. 00865 // 00866 // Input parameters: 00867 // <dl compact> 00868 // <dt><b>n</b> 00869 // <dd> The length of the array x to be transformed. The method 00870 // is most efficient when n is a product of small primes. 00871 // <dt><b>x</b> 00872 // <dd> An array which contains the sequence to be transformed 00873 // <dt><b>wsave</b> 00874 // A work array which must be dimensioned at least 3*n+15. 00875 // in the program that calls sinqb. The wsave array must be 00876 // initialized by calling <src>sinqi(n,wsave)</src> and a 00877 // different wsave array must be used for each different 00878 // value of n. This initialization does not have to be 00879 // repeated so long as n remains unchanged thus subsequent 00880 // transforms can be obtained faster than the first. 00881 // </dl> 00882 // Output parameters: 00883 // <dl compact> 00884 // <dt><b>x</b> 00885 // <dd> for i=1,...,n<br> 00886 // x(i)= the sum from k=1 to k=n of<br> 00887 // 4*x(k)*sin((2k-1)*i*pi/(2*n))<br> 00888 // 00889 // a call of sinqb followed by a call of 00890 // sinqf will multiply the sequence x by 4*n. 00891 // Therefore sinqf is the unnormalized inverse 00892 // of sinqb. 00893 // <dt><b>wsave</b> 00894 // <dd> Contains initialization calculations which must not 00895 // be destroyed between calls of sinqb or sinqf. 00896 // </dl> 00897 // <group> 00898 static void sinqb(Int n, Float* x, Float* wsave); 00899 static void sinqb(Int n, Double* x, Double* wsave); 00900 // </group> 00901 00902 // <group> 00903 static void cosqi(Int n, Float* wsave); 00904 static void cosqi(Int n, Double* wsave); 00905 // </group> 00906 // <group> 00907 static void cosqf(Int n, Float* x, Float* wsave); 00908 static void cosqf(Int n, Double* x, Double* wsave); 00909 // </group> 00910 // <group> 00911 static void cosqb(Int n, Float* x, Float* wsave); 00912 static void cosqb(Int n, Double* x, Double* wsave); 00913 // </group> 00914 }; 00915 00916 } //# NAMESPACE CASACORE - END 00917 00918 #endif