cosine transform for periodic functions with reflection symmetry

Computer Physics Communications 63 (1991) 569—577 North-Holland

569

Fast sine/cosine transform for periodic functions with reflection symmetry Anna Besprozvannaya and David J. Tannor Department of Chemistry and Biochemistry. University of Notre Dame, Notre Dame, IN 46556, USA Received 3 July 1990; in final form 28 July 1990

A modified fast cosine transform (FCT) algorithm is presented featuring the following three properties: (1) the entire calculation is performed using arrays half the size of what would be required using a common fast Fourier transform (FFT); (2) the result of the FCT is identical to that which would be obtained using an FFT and an array double the size; (3) the FCT is its own inverse. The algorithm requires a minor modification to the FCT algorithm in Numerical Recipes by W.H. Press et al., which satisfies only the first of the above three properties. An extension of fast sine/cosine transform algorithms to the case of complex data is included. Programs written in Fortran 77 for both real and complex versions are provided.

PROGRAM SUMMARY Titles of programs: ROFT/REFT/COFT/CEFT (for odd/ even real/complex fast transform) Catalogue number: ABTX Programs obtainable from: CPC Program Library, Queen’s University of Belfast, N.Ireland (see application form in this issue) Computer for which the program is designed and others on which it is operable: any computer with a Fortran 77 compiler and any FFT package Computer: Stellar GS2000; Installation: University of Notre Dame, Notre Dame, IN, USA

Nature of problem In many physical applications one is interested in obtaining a Fourier transform of a function having a particular reflection symmetry using only half of the grid points, thus saving both computation time and space without any loss of information. Method of solution The algorithm is based on fast sine/cosine transforms [1, p.401] with some modifications which make it applicable to a wider class of real and complex functions. Restrictions on the complexity of the problem The number of grid points should be a power of two in the odd case and a power of two plus one in the even case. We strongly recommend all the data to be double precision to achieve sufficient accuracy in the computations.

Operating system: Stellix (based on UNIX System V.3). Programming language used: Fortran 77. High speed storage required: 15 Kbytes words No. of bits in a word: 32

Typical running time It depends greatly on the computer type, the compiler optimization and the efficiency of the FFT routine used. For 5 forward/inverse transform operations example, to perform i0 of 33 elements using a highly-optiwith an even data array mized Stellar FFT library, took about 40 s for real and 80 s for complex double precision data.

No. of lines in programs (in total): 190 Keywords: sine/cosine fast transform, periodic real/complex function, odd/even symmetry, inverse transform

0010-4655/91/$03.50 © 1991

—

References [1] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes (Cambridge Univ. Press, Cambridge, MA, 1986).

Elsevier Science Publishers B.V. (North-Holland)

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transform for periodic functions

570

LONG WRITE-UP 1. Introduction In many physical applications one is interested in the Fourier transform of a function with a particular symmetry. Two very common examples are: f(x) has infinite translation symmetry, f(x + nX) =f(x) (n = cc,..., cc), and f(x) has even or odd reflection symmetry, f(—x) = ±f(x). In order to take full advantage of these symmetries in performing FFTs, to save both CPU time and storage, one must first understand the corresponding symmetry in reciprocal space. These corresponding symmetries, for the examples listed above, are as follows: —

—

—

f(x+nX) =f(x) ~-*F(k) where F(k) = fXf(x)e_2~11kx dx,

=F(k)~kk~,

f(—x)= ±f(x)4-~F(—k)=±F(k).

k’= —cc,..., —1,0, 1,...,cc,

(1.1)

(1.2)

The first of these symmetries leads to a Fourier series, rather than a Fourier transform representation in reciprocal space. Since the FF1’ is a Fourier series representation, it is therefore capable of exactly describing the k-space representation of a translationally infinite periodic function. The function f(x) need only be represented on a single unit cell, and no information is lost in either x- or k-space. The use of this periodic boundary condition of the FFT has been extensively applied in time-dependent molecule— surface scattering calculations [2]. It also emerges in the use of the FFT in non-cartesian coordinate systems, where typically f(x) is periodic in one or more angular variables [3]. The second of these symmetries, odd or even reflection synmletry, arises in many physical applications, often in conjunction with the translational symmetry discussed above. It stands to reason that when reflection symmetry is present, one may reduce the x-space range by a factor of two; although this will reduce the number of values in k-space by a factor of two, this need not lead to any loss of information, since the values in k-space are all duplicated, as may be seen from eq. (1.2). The transformation between f(x) and F(k) is accomplished via a fast sine/cosine transform (FST/FCT) for odd/even reflection symmetry, respectively. A variety of FST/FCT routines may be found in standard computer libraries and in numerical methods books. In appraising these routines the following simple criteria were examined: (1) the entire calculation be performed on a grid of half the range; (2) the result of the FST/FCT, up to a normalization factor, be identical to that which would be obtained using an FFT on a grid double the size. Although a suitable FST was easily found [1], we were unable to find an FCT that satisfied these simple criteria. The routine in IMSL, for example, doubles the range of x and proceeds to employ an FFT on the doubled function, hence providing no grid savings. The routine described in ref. [1] apparently assumes that f 0 = fx/2 = 0, where f~and fx/2 are values of the function f at the beginning and at the middle points of the periodicity interval X. This is unnecessarily restrictive and leads to incorrect answers if half the range of f(x) is used and f(x) does not satisfy this symmetry requirement. In this paper, a modified FCT algorithm is presented which satisfies both ~f the above criteria. It has a third attractive feature, in that it is its own inverse (whereas the routine in ref. [1] is not). The algorithm involves only a minor change in the FCT of ref. [1]. The crucial, albeit simple, algorithmic modification on which the routine is based is not new [4]. The main object of this paper is, first, to alert researchers who may be major users of FFTs to the potential savings, but also to the deficiencies that are connected with many of the existing FCT routines; and second, to provide the algorithm and the subroutine for a safe and

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transform for periodicfunctions

571

flexible (unrestricted boundary condition on f(x)) FCT satisfying all the desirable criteria described above. 2. Algorithm Our task was to obtain the fast Fourier transform (FFT) of an odd/even real periodic function only half of the grid points. Let us define the forward FFT as 2N—1

f

using

2N—1

~

~

f~e~”~~=

j=0

f~(cos(’rrjk/N)—isin(’rrjk/N)),

k=0,...,2N—1,

(2.1)

j=0

provided that function f is sampled in 2 N points of the periodicity interval X (N must be a power of 2). Unlike ref. [1, ch.12], we use a negative exponent so eq. (2.1) corresponds to the built-in forward FFT of most computers and software packages. Then for odd f N-I

f,

=-2i ~

F’k

sin(.~ijk/N).

F2Nk—

k=0,...,N—1,

_~dd,

(2.2)

j—0

and for even

f

~

F

(2.3)

k=0,...,N.

2N_kFk,

(Note that eq. (2.3) requires N + 1, rather than N, sample points, the Nth one being the point of symmetry of the even function.) Up to the factor —2i, eq. (2.2) can be obtained using the sine transform [1, eqs. (12.3.12)—(12.3.16)], the only difference being the negative sign of sin(2~rjk/N)due to eq. (2.1). To calculate eq. (2.3) we use a modified algorithm of the cosine transform [1, eqs. 2(12.3.17)—(12.3.19)]. instead of f~. Namely, constructing the auxiliary function Then, applying y[l, (12.3.17)1. a real-to-complex we set Yo toFFT be (f~ routine, +JN)/we get real and imaginary parts of the result as N—i

R,

N—i

~ y~cos(2irfl/N)

N—i

Yo + ~ y 1 cos(2’rrjl/N)

j=0

j=1

=

~f0+ ~

+

~

f, cos(2rrjl/N)

=

j=1

(2.4) and N—i

11=

~ y~,sin(21Tj//N),

—

1= 0,..., N/2

—

1

j=0

and by analogy with ref. [1, eq.(l2.3.l8)], Ii

I~,

—

—

even

I

even \

,,~F21~1 —

—

1

—

~F21_1 ~

even

1 —

even

(2.5)

Thus .peven2R even

— —

(2.6) even

F21÷1 F21.1

—

2I~,

/=0,...,N/2—1.

(2.7)

572

A. Besprozvannaya, D.J. Tannor

/

Fast sine/cosine transform forperiodic functions

The initial value of the recursion (2.7) is calculated separately as N-i

~fo~~f~+

~

N-i

~ f1cos(~~j/N)=2

Yo~fN~

j=1

f1cos(lTj/N)

~

(2.7’)

.

J=i

1,we took advantage of the fact that standard real-to-complex FFT To obtain last+ 1, value FJ~’ routines returnthe N/2 rather than N/2, complex components: the first and the last components are real and, therefore, are stored as real and imaginary parts of the first complex element, respectively. Thus, in our notation, RN/i,

therefore E’eVefl_’~D

‘N

_‘~T

Li~N/2~ui.

Unlike the cosine transform algorithm given in ref. [1], this FCT is its own inverse, i.e. applying it twice, we get the initial data scaled upwards by a factor of 2N. The FST/FCT algorithm can be used in the case of complex data. The auxiliary complex array y~ (j = 0,..., N 1) is constructed as before, but after applying complex-to-complex FFT we have to make one more computational step to obtain the sums —

N—i

N—i

y 1 cos(2’rrjl/N)

and

j=O

~ y1 sin(2’TTjl/N) j=O

separately. Namely, let y1

=

N—i

a1 + ib~,then the result of FFT N—i

~ y1e~”~”~=~ y~(cos(2’rrjl/N)—isin(2iiJ1/N)) 1=0

j=O

b~sin(2lTjl/N)] comes

Out

+

i~

[-ai

sin(2iij//N)

+

b~cos(211j1/N)1,

“mixed”; however, it can be easily verified that

N—i

N—i

y1 cos(2’rrjl/N)

=

j=O N—i

~ j=0

(a~+ ib1.) cos(2’rrjl/N)

=

~(Z,

+ ZN_i),

N—i

y~sin(2’rrjl/N) = 1=0

~ 1=0

(a~+ ui1) sin(2’rrjl/N) = ~i(Z1

—

ZN_i).

Thus, according to eqs. (2.2)—(2.8), for odd functions (2.9)

j~O~ddZl_ZN/, =

F~±’~ i(Z1+ —

ZN_i),

(2.10)

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transform for periodicfunctions

573

and for even ones, F

(2.11)

21=Zl+ZN1,

F2~’i~ =F~’~ +i(ZI—ZN_,);

(2.12)

the recursion starts from F~’= —iZ0

(2.13)

in the first case and from ~ also be noted that

which is calculated according to eq. (2.7’), in the second case. It should (2.14)

and 2ZN/ F~/~=

2.

(2.15)

As in the case of real data, the complex FST/FCT is its own inverse. However, the scaling factor is 2 N for the FCT and —2 N for the FST.

3. Program description The algorithms were implemented in programs ROFT, REFT, COFT and CEFT, standing for real odd, real even, complex odd and complex even FT respectively. They are provided in the form of subroutines with a specific set of parameters, assuming that the main program is supplied by the user. The calling sequences are the following: CALL ROFT(Y,N,SINE); CALL REFT(Y,N,SINE,COSINE); CALL COFT(Z,N,IP,F,SINE); CALL CEFT(Z,N,IP,F,SINE,COSINE). The main program should contain the declaration of the one-dimensional initial data array Y of REAL*8 or Z of COMPLEX *16 type and length N in the odd case or Ni = N + 1 in the even case (N must be a power of two: N = 2~’;in routines COFT and CEFT IP is also included in the parameter list). The initial array corresponding to the half of the function 1 (plus the additional middle point in even case) is replaced by the transform itself (for real even functions) or by the transform multiplied by i (for real odd functions). In the complex case another array F of the same length should be declared to contain the transform. Two additional arrays, SINE and COSINE, both of length M = N/2 and of the same type as the initial data, are intended for calculating sin(jrr/N) and cos(j’rr/N), j = 0,... ,N/2 1, [1, eqs. (12.3.12), (12.3.17) and (12.3.19)] in the main program, rather than inside the subroutine (as in ref. [1]). Calculating them in the subroutine slows down the performance, especially in the case of a large number of subroutine calls. —

4. Test run As an example we chose functions 11 = sin x + sin 2x + sin 3x, 12 = cos x + cos 2x + cos 3x, f~ =11 + if~and 14 = 12 + if2 in order to test ROFT, REFT, COFT and CEFT respectively. The results of the test

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transformfor periodicfunctions

574

runs are placed in tables, at the end of this paper, together with the results obtained from the FFT routines on arrays twice as big. In all cases they are identical (after normalization to the same scale) *~

Acknowledgements This work was supported in part by the National Science Foundation under grant No. NSF-CHE8996226-01.

References [1] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes (Cambridge Univ. Press, Cambridge, MA, 1986). [2] RB. Gerber, R. Kosloff and M. Berman, Comput. Phys. Rep. 5 (1986) 59. [3] Yun Shi and D.J. Tannor, J. Chem. Phys. 92 (1990) 2517. [4] M. Pickering, An Introduction to Fast Fourier Transform Methods for Partial Differential Equations, with Applications (Wiley, New York, 1986).

*

One should be careful with the normalization factor since it depends on the particular FFT routine used.

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transform for periodic functions

TEST RUN OUTPUT FORWARD

TRANSFORM

real odd function fl=sirt(x)+sin(2x)+sin(3x) COMMON FFT

SINE TRANSFORM 0,0000000000000 8.0000000000000 8.0000000000000 8.0000000000000

0.0000000000000 0.0000000000000 0.0000000000000

I I

( ( ( C

I

0.0000000000000 0.0000000000000

0.0000000000000 0.0000000000000 I ( 0.0000000000000 I ( 0.0000000000000 I C 0.0000000000000

0.0000000000000 —8.0000000000000 —8.0000000000000 —8.0000000000000

, , , ,

0.0000000000000 0.0000000000000

, ,

, 0.0000000000000 0.0000000000000 ( 0.0000000000000 , 0.0000000000000 * Note that the result of the sine transform is exactly the result of the FFT multiplied by i.

real even function f2=cos(x)+cos(2x)+cos(3x) MOD. COSINE TRANSFORM I 0.0000000000000 8.0000000000000

8.0000000000000 8.0000000000000 0.0000000000000

I I

0.0000000000000

I

0.0000000000000 0.0000000000000 0.0000000000000

I

COMMON FFT ( ( C C C C ( C

0.0000000000000 8.0000000000000 8.0000000000000 8.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000

, , , ,

, , ,

0.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000 0.0000000000000

I

INVERSE

TRANSFORM (NORMALIZED)

SINE TRANSFORM

COMMON FFT

0.0000000000000

I

0.0000000000000

2.0136697460629

I

2.0136697460629 I —0.5994561836898

2.4142135623731 1.2483028813327 0.0000000000000 —0.1659106810404 0.4142135623731 0.5994561836898

I I I I

2.4142135623731 —0.4142135623731 1.2483028813327 0.1659106810404 0.0000000000000 I 0.0000000000000 —0.1659106810404 —1.2483028813327 0.4142135623731 —2.4142135623731 0.5994561836898 I —2.0136697460629

I

MOD. COSINE TRANSFORNI

I

0.0000000000000

COMMON FFT

3.0000000000000

3.0000000000000

2.0136697460629

2.0136697460629 I

—0.5994561836898

0.0000000000000

0.0000000000000 —1.2483028813327 —1.0000000000000

I

0.0000000000000 —0.1659106810404 —1.2483028813327 I —1.0000000000000 —1.0000000000000 I —1.2483028813327

—0.1659106810404

I

—0.1659106810404 I

0.0000000000000 —0.5994561836898 —1.0000000000000

0.0000000000000 I I

0.0000000000000

I

2.0136697460629

—0.5994561836898

I

3.0000000000000

—1.0000000000000

I

) 1

575

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transform for periodicfunctions

576

FORWARD

TRANSFORM

complex odd function f3 SINE TRANSFORM

I

0.0000000000000, 0.0000000000000) 8.0000000000000, —8.0000000000000) 8.0000000000000, —8.0000000000000) 8.0000000000000, —8.0000000000000) 0.0000000000000, 0.0000000000000) 0.0000000000000, 0.0000000000000) 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000)

I I I I I I I I I I I I I I I I

COMMON FFT ( 0.0000000000000, 0.0000000000000) ( 8.0000000000000, —8.0000000000000) ( 8.0000000000000, —8.0000000000000) ( 8.0000000000000, —8.0000000000000) C 0.0000000000000, 0.0000000000000) ( 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000)

complex even function f4 COSINE TRANSFORM 0.0000000000000, 0.0000000000000) 8.0000000000000, 8.0000000000000) 8.0000000000000, 8.0000000000000) 8.0000000000000, 8.0000000000000) 0.0000000000000, 0.0000000000000) 0.0000000000000, 0.0000000000000) ( 0.0000000000000, 0.0000000000000) ( 0.0000000000000, 0.0000000000000) 0.0000000000000, 0.0000000000000)

fl + i*fl

=

=

I I I I I I I I I I I I I I I I f2

C 0.0000000000000, 0.0000000000000) ( 0.0000000000000, 0.0000000000000) ( 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000) (—8.0000000000000, 8.0000000000000) (—8.0000000000000, 8.0000000000000) (—8.0000000000000, 8.0000000000000) +

i*f2

I

COMMON FFT I I I I

I

C 0.0000000000000, 0.0000000000000) C 8.0000000000000, 8.0000000000000) I C 8.0000000000000, I 8.0000000000000) I ( 8.0000000000000, I 8.0000000000000) I ( 0.0000000000000, I 0.0000000000000) I C 0.0000000000000, I 0.0000000000000) I C 0.0000000000000, I 0.0000000000000) I C 0.0000000000000, I 0.0000000000000) I C 0.0000000000000, I 0.0000000000000) I

I I I I I I I I I I I I I I I I

( 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000) C 0.0000000000000, 0.0000000000000) C 0.O000000000000, 0.0000000000000) C 8.0000000000000, 8.0000000000000) C 8.0000000000000, 8.0000000000000) C 8.0000000000000, 8.0000000000000)

A. Besprozvannaya, D.J. Tannor / Fast sine/cosine transform forperiodic functions INVERSE

SINE TRANSFORM C 0.0000000000000, 0.0000000000000) (—2.0136697460629, —2.0136697460629) (—2.4142135623731, —2.4142135623731) (—1.2483028813327, —1.2483028813327) C 0.0000000000000, 0.0000000000000) C 0.1659106810404, 0.1659106810404) (—0.4142135623731, —0.4142135623731) (—0.5994561836898, —0.5994561836898) *

I

TRANSFORM

(NORMALIZED)

COMMON FFT

I C 0.0000000000000, I 0.0000000300000) I C 2.0136697460629, I 2.0136697460629) I C 2.4142135623731, I 2.4142135623731) I ( 1.2483028813327, I 1.2483028813327) I C 0.0000000000000, I 0.0000000000000) I (—0.1659106810404, I —0.1659106810404) I C 0.4142135623731, I 0.4142135623731) I ( 0.5994561836898, I 0.5994561836898)

I I I I I I I I I I I I I I I I

C 0.0000000000000, 0.0000000000000) (—0.5994561836898, —0.5994561836898) (—0 .4142135623731, —0.4142135623731) C 0.1659106810404, 0.1659106810404) C 0.0000000000000, 0.0000000000000) (—1.2483028813327, —1.2483028813327) (—2.4142135623731, —2.4142135623731) (—2.0136697460629, —2.0136697460629)

Note that the inverse sine transform comes out with the negative sign

COSINE TRANSFORM

I

C 3.0000000000000, 3.0000000000000) C 2.0136697460629, 2.0136697460629) C 0.0000000000000, 0.0000000000000)

I

COMMON FFT C 3.0000000000000, 3.0000000000000) C 2.0136697460629, 2.0136697460629) ( 0.0000000000000, 0.0000000000000) I I I I I

I

C—0.5994561836898, —0.5994561836898) I ( 0.0000000000000, I 0.0000000000000) I (—0.1659106810404, I —0.1659106810404) I

(—1.2483028813327, I (—1.2483028813327, I (—1.0000000000000, —1.2483028813327) I —1.2483028813327) I —1.0000000000000) (—1.0000000000000, —1.0000000000000) (—0.1659106810404, —0.1659106810404) C 0.0000000000000, 0.0000000000000) (—0.5994561836898, —0.5994561836898) (—1.0000000000000, —1.0000000000000)

I I I I I I

I

I

I

I

(—1.0000000000000, —1.0000000000000) (—0.1659106810404, —0.1659106810404) ( 0.0000000000000, 0.0000000000000) (—0.5994561836898, —0.5994561836898) (—1.0000000000000, —1.0000000000000)

I (—1.2483028813327, I —1.2483028813327) I C 0.0000000000000, I 0.0000000000000) I C 2.0136697460629, I 2.0136697460629) I I I I

577