casacore
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FFTPack.h
Go to the documentation of this file.
1 //# FFTPack.h: C++ wrapper functions for Fortran FFTPACK code
2 //# Copyright (C) 1993,1994,1995,1997,1999,2000,2001
3 //# Associated Universities, Inc. Washington DC, USA.
4 //#
5 //# This library is free software; you can redistribute it and/or modify it
6 //# under the terms of the GNU Library General Public License as published by
7 //# the Free Software Foundation; either version 2 of the License, or (at your
8 //# option) any later version.
9 //#
10 //# This library is distributed in the hope that it will be useful, but WITHOUT
11 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 //# License for more details.
14 //#
15 //# You should have received a copy of the GNU Library General Public License
16 //# along with this library; if not, write to the Free Software Foundation,
17 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 //#
19 //# Correspondence concerning AIPS++ should be addressed as follows:
20 //# Internet email: aips2-request@nrao.edu.
21 //# Postal address: AIPS++ Project Office
22 //# National Radio Astronomy Observatory
23 //# 520 Edgemont Road
24 //# Charlottesville, VA 22903-2475 USA
25 //#
26 
27 //# $Id$
28 
29 #ifndef SCIMATH_FFTPACK_H
30 #define SCIMATH_FFTPACK_H
31 
32 #include <casacore/casa/aips.h>
33 
34  //# The SGI compiler with -LANG:std has some trouble including both Complexfwd.h
35  //# and Complex.h so we bypass the problem by include Complex.h only.
36 #if defined(AIPS_USE_NEW_SGI)
38 #else
40 #endif
41 
42 #warning "FFTPack is deprecated and will be removed in a future version of casacore. Please use FFTW."
43 
44 namespace casacore { //# NAMESPACE CASACORE - BEGIN
45 
46 // <summary>C++ interface to the Fortran FFTPACK library</summary>
47 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="" demos="">
48 // </reviewed>
49 // <synopsis>
50 // The static functions in this class are C++ wrappers to the Fortran FFTPACK
51 // library. This library contains functions that perform fast Fourier
52 // transforms (FFT's) and related transforms.
53 
54 // An additional purpose of these definitions is to overload the functions so
55 // that C++ users can access the functions in either fftpak (single precision)
56 // or dfftpack (double precision) with identical function names.
57 
58 // These routines only do one-dimensional transforms with the first element of
59 // the array being the "origin" of the transform. The <linkto
60 // class="FFTServer">FFTServer</linkto> class uses some of these functions to
61 // implement multi-dimensional transforms with the origin of the transform
62 // either at the centre or the first element of the Array.
63 
64 // You must initialise the work array <src>wsave</src> before using the forward
65 // transform (function with a suffix of f) or the backward transform (with a
66 // suffix of b).
67 
68 // The transforms done by the functions in this class can be categorised as
69 // follows:
70 // <ul>
71 // <li> Complex to Complex Transforms<br>
72 // Done by the cttfi, cfftf & cfftb functions
73 // <li> Real to Complex Transforms<br>
74 // Done by the rffti, rfftf & rfftb functions. A simpler interface is
75 // provided by the ezffti, ezfftf & ezfftb functions. The 'ez' functions
76 // do not destroy the input array and provide the result in a slightly
77 // less packed format. They are available in single precision only and
78 // internally use the rfft functions.
79 // <li> Sine Transforms<br>
80 // Done by the sinti & sint functions. As the sine transform is its own
81 // inverse there is no need for any distinction between forward and
82 // backward transforms.
83 // <li> Cosine Transforms<br>
84 // Done by the costi & cost functions. As the cosine transform is its own
85 // inverse there is no need for any distinction between forward and
86 // backward transforms.
87 // <li> Sine quarter wave Transforms<br>
88 // Done by the sinqi, sinqf & sinqb functions.
89 // <li> Cosine quarter wave Transforms<br>
90 // Done by the cosqi, cosqf & cosqb functions.
91 // </ul>
92 
93 
94 // <note role=warning>
95 // These functions assume that it is possible to convert between Casacore numeric
96 // types and those used by Fortran. That it is possible to convert between
97 // Float & float, Double & double and Int & int.
98 // </note>
99 
100 // <note role=warning>
101 // These function also assume that a Complex array is stored as pairs of
102 // floating point numbers, with no intervening gaps, and with the real
103 // component first ie., <src>[re0,im0,re1,im1, ...]</src> so that the following
104 // type casts work,
105 // <srcblock>
106 // Complex* complexPtr;
107 // Float* floatPtr = (Float* ) complexPtr;
108 // </srcblock>
109 // and allow a Complex number to be accessed as a pair of real numbers. If this
110 // assumption is bad then float Arrays will have to generated by copying the
111 // complex ones. When compiled in debug mode mode the functions that require
112 // this assumption will throw an exception (AipsError) if this assumption is
113 // bad. Ultimately this assumption about Complex<->Float Array conversion
114 // should be put somewhere central like Array2Math.cc.
115 // </note>
116 
117 // </synopsis>
118 
119 
120 class FFTPack
121 {
122 public:
123 // cffti initializes the array wsave which is used in both <src>cfftf</src> and
124 // <src>cfftb</src>. The prime factorization of n together with a tabulation of
125 // the trigonometric functions are computed and stored in wsave.
126 //
127 // Input parameter:
128 // <dl compact>
129 // <dt><b>n</b>
130 // <dd> The length of the sequence to be transformed
131 // </dl>
132 // Output parameter:
133 // <dl compact>
134 // <dt><b>wsave</b>
135 // <dd> A work array which must be dimensioned at least 4*n+15
136 // The same work array can be used for both cfftf and cfftb
137 // as long as n remains unchanged. Different wsave arrays
138 // are required for different values of n. The contents of
139 // wsave must not be changed between calls of cfftf or cfftb.
140 // </dl>
141 // <group>
142 static void cffti(Int n, Float* wsave);
143 static void cffti(Int n, Double* wsave);
144 //Here is the doc from FFTPack 5.1
145 //You can convert the linguo from fortran to C/C++
146 /* Input Arguments
147 
148  L Integer number of elements to be transformed in the first
149  dimension. The transform is most efficient when L is a
150  product of small primes.
151 
152 
153  M Integer number of elements to be transformed in the second
154  dimension. The transform is most efficient when M is a
155  product of small primes.
156 
157  LENSAV Integer dimension of WSAVE array. LENSAV must be at least
158  2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8.
159 
160 
161  Output Arguments
162 
163  WSAVE Real work array with dimension LENSAV, containing the
164  prime factors of L and M, and also containing certain
165  trigonometric values which will be used in routines
166  CFFT2B or CFFT2F.
167 
168 
169  WSAVE Real work array with dimension LENSAV. The WSAVE array
170  must be initialized with a call to subroutine CFFT2I before
171  the first call to CFFT2B or CFFT2F, and thereafter whenever
172  the values of L, M or the contents of array WSAVE change.
173  Using different WSAVE arrays for different transform lengths
174  or types in the same program may reduce computation costs
175  because the array contents can be re-used.
176 
177 
178  IER Integer error return
179  = 0 successful exit
180  = 2 input parameter LENSAV not big enough
181  = 20 input error returned by lower level routine
182 */
183 
184 static void cfft2i(const Int& n, const Int& m, Float *& wsave, const Int& lensav, Int& ier);
185 // </group>
186 
187 // cfftf computes the forward complex discrete Fourier
188 // transform (the Fourier analysis). Equivalently, cfftf computes
189 // the Fourier coefficients of a complex periodic sequence.
190 // the transform is defined below at output parameter c.
191 //
192 // The transform is not normalized. To obtain a normalized transform
193 // the output must be divided by n. Otherwise a call of cfftf
194 // followed by a call of cfftb will multiply the sequence by n.
195 //
196 // The array wsave which is used by <src>cfftf</src> must be
197 // initialized by calling <src>cffti(n,wsave)</src>.
198 //
199 // Input parameters:
200 // <dl compact>
201 // <dt><b>n</b>
202 // <dd> The length of the complex sequence c. The method is
203 // more efficient when n is the product of small primes.
204 // <dt><b>c</b>
205 // <dd> A complex array of length n which contains the sequence to be
206 // transformed.
207 // <dt><b>wsave</b>
208 // <dd> A real work array which must be dimensioned at least 4n+15
209 // by the program that calls cfftf. The wsave array must be
210 // initialized by calling <src>cffti(n,wsave)</src> and a
211 // different wsave array must be used for each different
212 // value of n. This initialization does not have to be
213 // repeated so long as n remains unchanged thus subsequent
214 // transforms can be obtained faster than the first.
215 // The same wsave array can be used by cfftf and cfftb.
216 // </dl>
217 // Output parameters:
218 // <dl compact>
219 // <dt><b>c</b>
220 // <dd> for j=1,...,n<br>
221 // c(j)=the sum from k=1,...,n of<br>
222 // c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)<br>
223 // where i=sqrt(-1)<br>
224 // <dt><b>wsave</b>
225 // <dd> Contains initialization calculations which must not be
226 // destroyed between calls of cfftf or cfftb
227 // </dl>
228 // <group>
229 static void cfftf(Int n, Complex* c, Float* wsave);
230 static void cfftf(Int n, DComplex* c, Double* wsave);
231 
232 //Description from FFTPack 5.1
233 /*
234 Input Arguments
235 
236 
237  LDIM Integer first dimension of two-dimensional complex array C.
238 
239 
240  L Integer number of elements to be transformed in the first
241  dimension of the two-dimensional complex array C. The value
242  of L must be less than or equal to that of LDIM. The
243  transform is most efficient when L is a product of small
244  primes.
245 
246 
247  M Integer number of elements to be transformed in the second
248  dimension of the two-dimensional complex array C. The
249  transform is most efficient when M is a product of small
250  primes.
251 
252  C Complex array of two dimensions containing the (L,M) subarray
253  to be transformed. C's first dimension is LDIM, its second
254  dimension must be at least M.
255 
256  WSAVE Real work array with dimension LENSAV. WSAVE's contents
257  must be initialized with a call to subroutine CFFT2I before
258  the first call to routine CFFT2F or CFFT2B with transform
259  lengths L and M. WSAVE's contents may be re-used for
260  subsequent calls to CFFT2F and CFFT2B having those same
261  transform lengths.
262 
263 
264 
265  LENSAV Integer dimension of WSAVE array. LENSAV must be at least
266  2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8.
267 
268 
269  WORK Real work array.
270 
271 
272  LENWRK Integer dimension of WORK array. LENWRK must be at least
273  2*L*M.
274 */
275 static void cfft2f (const Int& ldim, const Int& L, const Int& M, Complex*& C, Float*& WSAVE, const Int& LENSAV,
276  Float *& WORK, const Int& LENWRK, Int& IER);
277 // </group>
278 
279 // cfftb computes the backward complex discrete Fourier
280 // transform (the Fourier synthesis). Equivalently, cfftb computes
281 // a complex periodic sequence from its Fourier coefficients.
282 // The transform is defined below with output parameter c.
283 //
284 // A call of cfftf followed by a call of cfftb will multiply the
285 // sequence by n.
286 //
287 // The array wsave which is used by <src>cfftb</src> must be
288 // initialized by calling <src>cffti(n,wsave)</src>.
289 //
290 // Input parameters:
291 // <dl compact>
292 // <dt><b>n</b>
293 // <dd> The length of the complex sequence c. The method is
294 // more efficient when n is the product of small primes.
295 // <dt><b>c</b>
296 // <dd> A complex array of length n which contains the sequence to be
297 // transformed.
298 // <dt><b>wsave</b>
299 // <dd> A real work array which must be dimensioned at least 4n+15
300 // in the program that calls cfftb. The wsave array must be
301 // initialized by calling <src>cffti(n,wsave)</src>
302 // and a different wsave array must be used for each different
303 // value of n. This initialization does not have to be
304 // repeated so long as n remains unchanged thus subsequent
305 // transforms can be obtained faster than the first.
306 // The same wsave array can be used by cfftf and cfftb.
307 // </dl>
308 // Output parameters:
309 // <dl compact>
310 // <dt><b>c</b>
311 // <dd> for j=1,...,n<br>
312 // c(j)=the sum from k=1,...,n of<br>
313 // c(k)*exp(i*(j-1)*(k-1)*2*pi/n)<br>
314 
315 // <dt><b>wsave</b>
316 // <dd> Contains initialization calculations which must not be
317 // destroyed between calls of cfftf or cfftb
318 // </dl>
319 // <group>
320 static void cfftb(Int n, Complex* c, Float* wsave);
321 static void cfftb(Int n, DComplex* c, Double* wsave);
322 
323 
324 //Documentation from FFTPack 5.1
325 /*
326  Input Arguments
327 
328  LDIM Integer first dimension of two-dimensional complex array C.
329 
330 
331  L Integer number of elements to be transformed in the first
332  dimension of the two-dimensional complex array C. The value
333  of L must be less than or equal to that of LDIM. The
334  transform is most efficient when L is a product of small
335  primes.
336 
337 
338  M Integer number of elements to be transformed in the second
339  dimension of the two-dimensional complex array C. The
340  transform is most efficient when M is a product of small
341  primes.
342 
343  C Complex array of two dimensions containing the (L,M) subarray
344  to be transformed. C's first dimension is LDIM, its second
345  dimension must be at least M.
346 
347 
348  WSAVE Real work array with dimension LENSAV. WSAVE's contents
349  must be initialized with a call to subroutine CFFT2I before
350  the first call to routine CFFT2F or CFFT2B with transform
351  lengths L and M. WSAVE's contents may be re-used for
352  subsequent calls to CFFT2F and CFFT2B with the same
353  transform lengths L and M.
354 
355  LENSAV Integer dimension of WSAVE array. LENSAV must be at least
356  2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8.
357 
358  WORK Real work array.
359 
360  LENWRK Integer dimension of WORK array. LENWRK must be at least
361  2*L*M.
362 
363  Output Arguments
364 
365 
366  C Complex output array. For purposes of exposition,
367  assume the index ranges of array C are defined by
368  C(0:L-1,0:M-1).
369 
370  For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given
371  in the traditional aliased form by
372 
373  L-1 M-1
374  C(I,J) = SUM SUM C(L1,M1)*
375  L1=0 M1=0
376 
377  EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M))
378 
379  And in unaliased form, the C(I,J)'s are given by
380 
381  LF MF
382  C(I,J) = SUM SUM C(L1,M1,K1)*
383  L1=LS M1=MS
384 
385  EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M))
386 
387  where
388 
389  LS= -L/2 and LF=L/2-1 if L is even;
390  LS=-(L-1)/2 and LF=(L-1)/2 if L is odd;
391  MS= -M/2 and MF=M/2-1 if M is even;
392  MS=-(M-1)/2 and MF=(M-1)/2 if M is odd;
393 
394  and
395 
396  C(L1,M1) = C(L1+L,M1) if L1 is zero or negative;
397  C(L1,M1) = C(L1,M1+M) if M1 is zero or negative;
398 
399  The two forms give different results when used to
400  interpolate between elements of the sequence.
401 
402  IER Integer error return
403  = 0 successful exit
404  = 2 input parameter LENSAV not big enough
405  = 3 input parameter LENWRK not big enough
406  = 5 input parameter L > LDIM
407  = 20 input error returned by lower level routine
408 */
409 
410 
411 
412 static void cfft2b (const Int& LDIM, const Int& L, const Int& M, Complex* & C, Float *& WSAVE, const Int& LENSAV,
413  Float*& WORK, const Int& LENWRK, Int& IER);
414 // </group>
415 
416 // rffti initializes the array wsave which is used in both <src>rfftf</src> and
417 // <src>rfftb</src>. The prime factorization of n together with a tabulation of
418 // the trigonometric functions are computed and stored in wsave.
419 //
420 // Input parameter:
421 // <dl compact>
422 // <dt><b>n</b>
423 // <dd> The length of the sequence to be transformed.
424 // </dl>
425 // Output parameter:
426 // <dl compact>
427 // <dt><b>wsave</b>
428 // <dd> A work array which must be dimensioned at least 2*n+15.
429 // The same work array can be used for both rfftf and rfftb
430 // as long as n remains unchanged. Different wsave arrays
431 // are required for different values of n. The contents of
432 // wsave must not be changed between calls of rfftf or rfftb.
433 // </dl>
434 // <group>
435 static void rffti(Int n, Float* wsave);
436 static void rffti(Int n, Double* wsave);
437 // </group>
438 
439 // rfftf computes the Fourier coefficients of a real perodic sequence (Fourier
440 // analysis). The transform is defined below at output parameter r.
441 //
442 // Input parameters:
443 // <dl compact>
444 // <dt><b>n</b>
445 // <dd> The length of the array r to be transformed. The method
446 // is most efficient when n is a product of small primes.
447 // n may change so long as different work arrays are provided
448 // <dt><b>r</b>
449 // <dd> A real array of length n which contains the sequence
450 // to be transformed
451 // <dt><b>wsave</b>
452 // <dd> A work array which must be dimensioned at least 2*n+15
453 // in the program that calls rfftf. The wsave array must be
454 // initialized by calling <src>rffti(n,wsave)</src> and a
455 // different wsave array must be used for each different
456 // value of n. This initialization does not have to be
457 // repeated so long as n remains unchanged thus subsequent
458 // transforms can be obtained faster than the first.
459 // The same wsave array can be used by rfftf and rfftb.
460 // </dl>
461 // output parameters
462 // <dl compact>
463 // <dt><b>r</b>
464 // <dd> r(1) = the sum from i=1 to i=n of r(i)<br>
465 // if n is even set l = n/2 , if n is odd set l = (n+1)/2<br>
466 // then for k = 2,...,l<br>
467 // r(2*k-2) = the sum from i = 1 to i = n of<br>
468 // r(i)*cos((k-1)*(i-1)*2*pi/n)<br>
469 // r(2*k-1) = the sum from i = 1 to i = n of<br>
470 // -r(i)*sin((k-1)*(i-1)*2*pi/n)<br>
471 // if n is even<br>
472 // r(n) = the sum from i = 1 to i = n of<br>
473 // (-1)**(i-1)*r(i)<br>
474 //
475 // note:
476 // this transform is unnormalized since a call of rfftf
477 // followed by a call of rfftb will multiply the input
478 // sequence by n.
479 // <dt><b>wsave</b>
480 // <dd> Contains results which must not be destroyed between
481 // calls of rfftf or rfftb.
482 // </dl>
483 // <group>
484 static void rfftf(Int n, Float* r, Float* wsave);
485 static void rfftf(Int n, Double* r, Double* wsave);
486 // </group>
487 
488 
489 // rfftb computes the real perodic sequence from its Fourier coefficients
490 // (Fourier synthesis). The transform is defined below at output parameter r.
491 //
492 // Input parameters:
493 // <dl compact>
494 // <dt><b>n</b>
495 // <dd> The length of the array r to be transformed. The method
496 // is most efficient when n is a product of small primes.
497 // n may change so long as different work arrays are provided
498 // <dt><b>r</b>
499 // <dd> A real array of length n which contains the sequence
500 // to be transformed
501 // <dt><b>wsave</b>
502 // <dd> A work array which must be dimensioned at least 2*n+15
503 // in the program that calls rfftb. The wsave array must be
504 // initialized by calling <src>rffti(n,wsave)</src> and a
505 // different wsave array must be used for each different
506 // value of n. This initialization does not have to be
507 // repeated so long as n remains unchanged thus subsequent
508 // transforms can be obtained faster than the first.
509 // The same wsave array can be used by rfftf and rfftb.
510 // </dl>
511 // Output parameters:
512 // <dl compact>
513 // <dt><b>r</b>
514 // <dd> for n even and for i = 1,...,n<br>
515 // r(i) = r(1)+(-1)**(i-1)*r(n)<br>
516 // plus the sum from k=2 to k=n/2 of<br>
517 // 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)<br>
518 // -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)<br>
519 // for n odd and for i = 1,...,n<br>
520 // r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of<br>
521 // 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)<br>
522 // -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)<br>
523 //
524 // note:
525 // this transform is unnormalized since a call of rfftf
526 // followed by a call of rfftb will multiply the input
527 // sequence by n.
528 // <dt><b>wsave</b>
529 // <dd> Contains results which must not be destroyed between
530 // calls of rfftb or rfftf.
531 // </dl>
532 // <group>
533 static void rfftb(Int n, Float* r, Float* wsave);
534 static void rfftb(Int n, Double* r, Double* wsave);
535 // </group>
536 
537 // ezffti initializes the array wsave which is used in both <src>ezfftf</src>
538 // and <src>ezfftb</src>. The prime factorization of n together with a
539 // tabulation of the trigonometric functions are computed and stored in wsave.
540 //
541 // Input parameter:
542 // <dl compact>
543 // <dt><b>n</b>
544 // <dd> The length of the sequence to be transformed.
545 // </dl>
546 // Output parameter:
547 // <dl compact>
548 // <dt><b>wsave</b>
549 // <dd> A work array which must be dimensioned at least 3*n+15.
550 // The same work array can be used for both ezfftf and ezfftb
551 // as long as n remains unchanged. Different wsave arrays
552 // are required for different values of n.
553 // </dl>
554 static void ezffti(Int n, Float* wsave);
555 
556 // ezfftf computes the Fourier coefficients of a real
557 // perodic sequence (Fourier analysis). The transform is defined
558 // below at output parameters azero, a and b. ezfftf is a simplified
559 // but slower version of rfftf.
560 //
561 // Input parameters:
562 // <dl compact>
563 // <dt><b>n</b>
564 // <dd> The length of the array r to be transformed. The method
565 // is most efficient when n is the product of small primes.
566 // <dt><b>r</b>
567 // <dd> A real array of length n which contains the sequence
568 // to be transformed. r is not destroyed.
569 // <dt><b>wsave</b>
570 // <dd> A work array which must be dimensioned at least 3*n+15
571 // in the program that calls ezfftf. The wsave array must be
572 // initialized by calling <src>ezffti(n,wsave)</src> and a
573 // different wsave array must be used for each different
574 // value of n. This initialization does not have to be
575 // repeated so long as n remains unchanged thus subsequent
576 // transforms can be obtained faster than the first.
577 // The same wsave array can be used by ezfftf and ezfftb.
578 // </dl>
579 // Output parameters:
580 // <dl compact>
581 // <dt><b>azero</b>
582 // <dd> The sum from i=1 to i=n of r(i)/n
583 // <dt><b>a,b</b>
584 // <dd> Real arrays of length n/2 (n even) or (n-1)/2 (n odd)<br>
585 // for n even<br>
586 // b(n/2)=0, and <br>
587 // a(n/2) is the sum from i=1 to i=n of (-1)**(i-1)*r(i)/n<br>
588 //
589 // for n even define kmax=n/2-1<br>
590 // for n odd define kmax=(n-1)/2<br>
591 // then for k=1,...,kmax<br>
592 // a(k) equals the sum from i=1 to i=n of<br>
593 // 2./n*r(i)*cos(k*(i-1)*2*pi/n)<br>
594 // b(k) equals the sum from i=1 to i=n of<br>
595 // 2./n*r(i)*sin(k*(i-1)*2*pi/n)<br>
596 // </dl>
597 static void ezfftf(Int n, Float* r, Float* azero, Float* a, Float* b,
598  Float* wsave);
599 
600 // ezfftb computes a real perodic sequence from its
601 // Fourier coefficients (Fourier synthesis). The transform is
602 // defined below at output parameter r. ezfftb is a simplified
603 // but slower version of rfftb.
604 //
605 // Input parameters:
606 // <dl compact>
607 // <dt><b>n</b>
608 // <dd> The length of the output array r. The method is most
609 // efficient when n is the product of small primes.
610 // <dt><b>azero</b>
611 // <dd> The constant Fourier coefficient
612 // <dt><b>a,b</b>
613 // <dd> Arrays which contain the remaining Fourier coefficients
614 // these arrays are not destroyed.
615 // The length of these arrays depends on whether n is even or
616 // odd.
617 // If n is even n/2 locations are required,
618 // if n is odd (n-1)/2 locations are required.
619 // <dt><b>wsave</b>
620 // <dd> A work array which must be dimensioned at least 3*n+15.
621 // in the program that calls ezfftb. The wsave array must be
622 // initialized by calling <src>ezffti(n,wsave)</src> and a
623 // different wsave array must be used for each different
624 // value of n. This initialization does not have to be
625 // repeated so long as n remains unchanged thus subsequent
626 // transforms can be obtained faster than the first.
627 // The same wsave array can be used by ezfftf and ezfftb.
628 // </dl>
629 // Output parameters:
630 // <dl compact>
631 // <dt><b>r</b>
632 // <dd> if n is even define kmax=n/2<br>
633 // if n is odd define kmax=(n-1)/2<br>
634 // then for i=1,...,n<br>
635 // r(i)=azero plus the sum from k=1 to k=kmax of<br>
636 // a(k)*cos(k*(i-1)*2*pi/n)+b(k)*sin(k*(i-1)*2*pi/n)<br>
637 // where<br>
638 // c(k) = .5*cmplx(a(k),-b(k)) for k=1,...,kmax<br>
639 // c(-k) = conjg(c(k))<br>
640 // c(0) = azero<br>
641 // and i=sqrt(-1)<br>
642 // </dl>
643 static void ezfftb(Int n, Float* r, Float* azero, Float* a, Float* b,
644  Float* wsave);
645 
646 // sinti initializes the array wsave which is used in
647 // <src>sint</src>. The prime factorization of n together with a tabulation of
648 // the trigonometric functions are computed and stored in wsave.
649 //
650 // Input parameter:
651 // <dl compact>
652 // <dt><b>n</b>
653 // <dd> The length of the sequence to be transformed. the method
654 // is most efficient when n+1 is a product of small primes.
655 // </dl>
656 // Output parameter:
657 // <dl compact>
658 // <dt><b>wsave</b>
659 // <dd> A work array with at least int(2.5*n+15) locations.
660 // Different wsave arrays are required for different values
661 // of n. The contents of wsave must not be changed between
662 // calls of sint.
663 // </dl>
664 // <group>
665 static void sinti(Int n, Float* wsave);
666 static void sinti(Int n, Double* wsave);
667 // </group>
668 
669 // sint computes the discrete Fourier sine transform
670 // of an odd sequence x(i). The transform is defined below at
671 // output parameter x.
672 // sint is the unnormalized inverse of itself since a call of sint
673 // followed by another call of sint will multiply the input sequence
674 // x by 2*(n+1).
675 // The array wsave which is used by sint must be
676 // initialized by calling <src>sinti(n,wsave)</src>.
677 //
678 // Input parameters:
679 // <dl compact>
680 // <dt><b>n</b>
681 // <dd> The length of the sequence to be transformed. The method
682 // is most efficient when n+1 is the product of small primes.
683 // <dt><b>x</b>
684 // <dd> An array which contains the sequence to be transformed
685 // <dt><b>wsave</b>
686 // <dd> A work array with dimension at least int(2.5*n+15)
687 // in the program that calls sint. The wsave array must be
688 // initialized by calling <src>sinti(n,wsave)</src> and a
689 // different wsave array must be used for each different
690 // value of n. This initialization does not have to be
691 // repeated so long as n remains unchanged thus subsequent
692 // transforms can be obtained faster than the first.
693 // </dl>
694 // Output parameters:
695 // <dl compact>
696 // <dt><b>x</b>
697 // <dd> for i=1,...,n<br>
698 // x(i) = the sum from k=1 to k=n<br>
699 // 2*x(k)*sin(k*i*pi/(n+1))<br>
700 //
701 // a call of sint followed by another call of
702 // sint will multiply the sequence x by 2*(n+1).
703 // Hence sint is the unnormalized inverse
704 // of itself.
705 //
706 // <dt><b>wsave</b>
707 // <dd> Contains initialization calculations which must not be
708 // destroyed between calls of sint.
709 // </dl>
710 // <group>
711 static void sint(Int n, Float* x, Float* wsave);
712 static void sint(Int n, Double* x, Double* wsave);
713 // </group>
714 
715 // costi initializes the array wsave which is used in
716 // <src>cost</src>. The prime factorization of n together with a tabulation of
717 // the trigonometric functions are computed and stored in wsave.
718 //
719 // Input parameter:
720 // <dl compact>
721 // <dt><b>n</b>
722 // <dd> The length of the sequence to be transformed. The method
723 // is most efficient when n-1 is a product of small primes.
724 // </dl>
725 // Output parameter:
726 // <dl compact>
727 // <dt><b>wsave</b>
728 // <dd> A work array which must be dimensioned at least 3*n+15.
729 // Different wsave arrays are required for different values
730 // of n. The contents of wsave must not be changed between
731 // calls of cost.
732 // </dl>
733 // <group>
734 static void costi(Int n, Float* wsave);
735 static void costi(Int n, Double* wsave);
736 // </group>
737 
738 // cost computes the discrete Fourier cosine transform
739 // of an even sequence x(i). The transform is defined below at output
740 // parameter x.
741 // cost is the unnormalized inverse of itself since a call of cost
742 // followed by another call of cost will multiply the input sequence
743 // x by 2*(n-1). The transform is defined below at output parameter x.
744 // The array wsave which is used by <src>cost</src> must be
745 // initialized by calling <src>costi(n,wsave)</src>.
746 //
747 // Input parameters:
748 // <dl compact>
749 // <dt><b>n</b>
750 // <dd> The length of the sequence x. n must be greater than 1.
751 // The method is most efficient when n-1 is a product of
752 // small primes.
753 // <dt><b>x</b>
754 // <dd> An array which contains the sequence to be transformed
755 // <dt><b>wsave</b>
756 // <dd> A work array which must be dimensioned at least 3*n+15
757 // in the program that calls cost. The wsave array must be
758 // initialized by calling <src>costi(n,wsave)</src> and a
759 // different wsave array must be used for each different
760 // value of n. This initialization does not have to be
761 // repeated so long as n remains unchanged thus subsequent
762 // transforms can be obtained faster than the first.
763 // </dl>
764 // Output parameters:
765 // <dl compact>
766 // <dt><b>x</b>
767 // <dd> for i=1,...,n<br>
768 // x(i) = x(1)+(-1)**(i-1)*x(n)<br>
769 // + the sum from k=2 to k=n-1<br>
770 // 2*x(k)*cos((k-1)*(i-1)*pi/(n-1))<br>
771 //
772 // a call of cost followed by another call of
773 // cost will multiply the sequence x by 2*(n-1)
774 // hence cost is the unnormalized inverse
775 // of itself.
776 // <dt><b>wsave</b>
777 // <dd> Contains initialization calculations which must not be
778 // destroyed between calls of cost.
779 // </dl>
780 // <group>
781 static void cost(Int n, Float* x, Float* wsave);
782 static void cost(Int n, Double* x, Double* wsave);
783 // </group>
784 
785 // sinqi initializes the array wsave which is used in both <src>sinqf</src> and
786 // <src>sinqb</src>. The prime factorization of n together with a tabulation of
787 // the trigonometric functions are computed and stored in wsave.
788 //
789 // Input parameter:
790 // <dl compact>
791 // <dt><b>n</b>
792 // <dd> The length of the sequence to be transformed. The method
793 // is most efficient when n is a product of small primes.
794 // </dl>
795 // Output parameter:
796 // <dl compact>
797 // <dt><b>wsave</b>
798 // <dd> A work array which must be dimensioned at least 3*n+15.
799 // The same work array can be used for both sinqf and sinqb
800 // as long as n remains unchanged. Different wsave arrays
801 // are required for different values of n. The contents of
802 // wsave must not be changed between calls of sinqf or sinqb.
803 // </dl>
804 // <group>
805 static void sinqi(Int n, Float* wsave);
806 static void sinqi(Int n, Double* wsave);
807 // </group>
808 
809 // sinqf computes the fast Fourier transform of quarter wave data. That is,
810 // sinqf computes the coefficients in a sine series representation with only
811 // odd wave numbers. The transform is defined below at output parameter x.
812 //
813 // sinqb is the unnormalized inverse of sinqf since a call of sinqf followed by
814 // a call of sinqb will multiply the input sequence x by 4*n.
815 //
816 // The array wsave which is used by sinqf must be initialized by calling
817 // <src>sinqi(n,wsave)</src>.
818 //
819 // Input parameters:
820 // <dl compact>
821 // <dt><b>n</b>
822 // <dd> The length of the array x to be transformed. The method
823 // is most efficient when n is a product of small primes.
824 // <dt><b>x</b>
825 // <dd> An array which contains the sequence to be transformed
826 // <dt><b>wsave</b>
827 // A work array which must be dimensioned at least 3*n+15.
828 // in the program that calls sinqf. The wsave array must be
829 // initialized by calling <src>sinqi(n,wsave)</src> and a
830 // different wsave array must be used for each different
831 // value of n. This initialization does not have to be
832 // repeated so long as n remains unchanged thus subsequent
833 // transforms can be obtained faster than the first.
834 // </dl>
835 // Output parameters:
836 // <dl compact>
837 // <dt><b>x</b>
838 // <dd> for i=1,...,n<br>
839 // x(i) = (-1)**(i-1)*x(n)<br>
840 // + the sum from k=1 to k=n-1 of<br>
841 // 2*x(k)*sin((2*i-1)*k*pi/(2*n))<br>
842 //
843 // a call of sinqf followed by a call of
844 // sinqb will multiply the sequence x by 4*n.
845 // therefore sinqb is the unnormalized inverse
846 // of sinqf.
847 // <dt><b>wsave </b>
848 // <dd> Contains initialization calculations which must not
849 // be destroyed between calls of sinqf or sinqb.
850 // </dl>
851 // <group>
852 static void sinqf(Int n, Float* x, Float* wsave);
853 static void sinqf(Int n, Double* x, Double* wsave);
854 // </group>
855 
856 // sinqb computes the fast Fourier transform of quarter
857 // wave data. that is, sinqb computes a sequence from its
858 // representation in terms of a sine series with odd wave numbers.
859 // the transform is defined below at output parameter x.
860 //
861 // sinqf is the unnormalized inverse of sinqb since a call of sinqb
862 // followed by a call of sinqf will multiply the input sequence x
863 // by 4*n.
864 //
865 // The array wsave which is used by <src>sinqb</src> must be
866 // initialized by calling <src>sinqi(n,wsave)</src>.
867 //
868 // Input parameters:
869 // <dl compact>
870 // <dt><b>n</b>
871 // <dd> The length of the array x to be transformed. The method
872 // is most efficient when n is a product of small primes.
873 // <dt><b>x</b>
874 // <dd> An array which contains the sequence to be transformed
875 // <dt><b>wsave</b>
876 // A work array which must be dimensioned at least 3*n+15.
877 // in the program that calls sinqb. The wsave array must be
878 // initialized by calling <src>sinqi(n,wsave)</src> and a
879 // different wsave array must be used for each different
880 // value of n. This initialization does not have to be
881 // repeated so long as n remains unchanged thus subsequent
882 // transforms can be obtained faster than the first.
883 // </dl>
884 // Output parameters:
885 // <dl compact>
886 // <dt><b>x</b>
887 // <dd> for i=1,...,n<br>
888 // x(i)= the sum from k=1 to k=n of<br>
889 // 4*x(k)*sin((2k-1)*i*pi/(2*n))<br>
890 //
891 // a call of sinqb followed by a call of
892 // sinqf will multiply the sequence x by 4*n.
893 // Therefore sinqf is the unnormalized inverse
894 // of sinqb.
895 // <dt><b>wsave</b>
896 // <dd> Contains initialization calculations which must not
897 // be destroyed between calls of sinqb or sinqf.
898 // </dl>
899 // <group>
900 static void sinqb(Int n, Float* x, Float* wsave);
901 static void sinqb(Int n, Double* x, Double* wsave);
902 // </group>
903 
904 // <group>
905 static void cosqi(Int n, Float* wsave);
906 static void cosqi(Int n, Double* wsave);
907 // </group>
908 // <group>
909 static void cosqf(Int n, Float* x, Float* wsave);
910 static void cosqf(Int n, Double* x, Double* wsave);
911 // </group>
912 // <group>
913 static void cosqb(Int n, Float* x, Float* wsave);
914 static void cosqb(Int n, Double* x, Double* wsave);
915 // </group>
916 };
917 
918 } //# NAMESPACE CASACORE - END
919 
920 #endif
static void cosqb(Int n, Float *x, Float *wsave)
static void cost(Int n, Float *x, Float *wsave)
cost computes the discrete Fourier cosine transform of an even sequence x(i).
int Int
Definition: aipstype.h:50
static void sint(Int n, Float *x, Float *wsave)
sint computes the discrete Fourier sine transform of an odd sequence x(i).
static void cfft2f(const Int &ldim, const Int &L, const Int &M, Complex *&C, Float *&WSAVE, const Int &LENSAV, Float *&WORK, const Int &LENWRK, Int &IER)
Description from FFTPack 5.1.
static void ezfftf(Int n, Float *r, Float *azero, Float *a, Float *b, Float *wsave)
ezfftf computes the Fourier coefficients of a real perodic sequence (Fourier analysis).
static void rfftb(Int n, Float *r, Float *wsave)
rfftb computes the real perodic sequence from its Fourier coefficients (Fourier synthesis).
static void cfftb(Int n, Complex *c, Float *wsave)
cfftb computes the backward complex discrete Fourier transform (the Fourier synthesis).
static void rffti(Int n, Float *wsave)
rffti initializes the array wsave which is used in both rfftf and rfftb.
static void sinqb(Int n, Float *x, Float *wsave)
sinqb computes the fast Fourier transform of quarter wave data.
C++ interface to the Fortran FFTPACK library.
Definition: FFTPack.h:120
static void cfft2i(const Int &n, const Int &m, Float *&wsave, const Int &lensav, Int &ier)
Here is the doc from FFTPack 5.1 You can convert the linguo from fortran to C/C++.
static void cfft2b(const Int &LDIM, const Int &L, const Int &M, Complex *&C, Float *&WSAVE, const Int &LENSAV, Float *&WORK, const Int &LENWRK, Int &IER)
Documentation from FFTPack 5.1.
static void sinqf(Int n, Float *x, Float *wsave)
sinqf computes the fast Fourier transform of quarter wave data.
static void ezfftb(Int n, Float *r, Float *azero, Float *a, Float *b, Float *wsave)
ezfftb computes a real perodic sequence from its Fourier coefficients (Fourier synthesis).
double Double
Definition: aipstype.h:55
static void cffti(Int n, Float *wsave)
cffti initializes the array wsave which is used in both cfftf and cfftb.
static void ezffti(Int n, Float *wsave)
ezffti initializes the array wsave which is used in both ezfftf and ezfftb.
float Float
Definition: aipstype.h:54
static void cosqi(Int n, Float *wsave)
static void cfftf(Int n, Complex *c, Float *wsave)
cfftf computes the forward complex discrete Fourier transform (the Fourier analysis).
static void sinqi(Int n, Float *wsave)
sinqi initializes the array wsave which is used in both sinqf and sinqb.
static void costi(Int n, Float *wsave)
costi initializes the array wsave which is used in cost.
static void cosqf(Int n, Float *x, Float *wsave)
static void sinti(Int n, Float *wsave)
sinti initializes the array wsave which is used in sint.
const Double c
Fundamental physical constants (SI units):
static void rfftf(Int n, Float *r, Float *wsave)
rfftf computes the Fourier coefficients of a real perodic sequence (Fourier analysis).