FFTPack.h

Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1