A software for evaluating local accuracy in the fourier transform

Mathematics and Computers in Simulation North-Holland Publishing Company

XXII (1980) 141-150

A SOFTWARE FOR EVALUATING

LOCAL ACCURACY IN THE FOURIER TRANSFORM

P. BOIS and J. VICNES Utliversitc'Pierrc etMarieCurie.7.i230Paris Cedex 05,France

With the evolution of conrputerstoward ever increasing computing speed, the numerical aZgorithms implemented on such computers are more and more complex, and require a very great number of ::omputations.The numerical errors inherent in floating-point arithmetic of the computer generate errors in the results which may be fairly great. In this paper we propose a software based on the I,erturbationMethod for estimating the locaL accuracy in the Fast Fourier Transform, both in the case of exact data (computing errors aLone and in the ease of experimentaL data (data errors and computing errors).

1.

INTRODUCTION

Since 1965 when Cooley and Tukey 121 worked out an algorithm (FFT) for calculating the discrete Fourier transform (DFT), a great many researchers have been interested in studying its accuracy. Indeed, the DFT which involves the computing of a sum of products of numbers in CL? cannot exactly be implemented on a computer. Each multiplication and addition introduces an error caused by the limited-accuracy arithmetic of the computer. Rather than giving the upper bounds of these computer errors as Gentleman and Sande 131 did, many have preferred to use a statistical model of these errors to predict their variances. Several researchers have analyzed these roundoff errors in fixed-point arithmetic C8, 9, 10, 151, and others in floating-point arithmetic C5, 11, 141. More recently Alt Cl1 generalized to the case of complex numbers the formula presented in [61 that gives the error in the computer calculation of a scalar product and applied it to the FFT algorithm. The overall numerical error has been shown to be approximately Log N or (Log N)2 in rounding.zhmetic,wl .2h N being the off, or chopping at-1 number of samples of the function for which the Fourier Transform is being computed. It should be pointed out that the formulas worked out by Alt Cl1 giving the noise-to-signal ratio in the computing of the DFT by the FFT algorithm are in agreement with those of Weinstein Cl31 and Kaneko and Liu C51. The techniques used up to now have served only to establish standards for computing errors, i.e. overall estimates of the noise-to-signal ratio. The use of Vignes'sPermutation-Perturbation Method r121 in &'applied to any algorithm for computing the DFT enables the accuracy to be determined for the real part and the imaginary

part as well each component

for the module as of this DFT.

2.

ACCURACY

EVALUATING

and phase of

IN ALGORITHMS

Like any mathematical transform,the DFT is a direct method which, when performedwith infinite accuracy, gives the exact Fourier Transform in a finite number of operations. But when a Fourier transform algorithm is carried out on a computer, it is performedwith limited-accuracy arithmetic which 'generates an error at each operation. The propagation of these errors irremediably affects the OFT. The accuracy of each component of the DFT can be evaluated b,y the Permutation-Perturbation Method described in 1121. Here is a very succinct behing this method: 2.1 Let us consider defined by:

summary

a direct

of the basic idea

algebraic

procedure

procedure (d, r, +, -, x, :, funct.) in which d= [R is a set of data, r-6 R is a single result, +, -, x, :, funct. are exact mathematical operators.

(1)

To implement this algebraic procedure, it is writtenin the form of an equivalent syntaxic program which is defined by: Procedure

(D, R, 0,

0,

*, /,

FONCT)

(2)

in which DC Fis the set of data, R E pis the result, /, FUNCT are the data0, o,*, processing arithmetic operators, pis the set of values represented in the computer in the standardized floatingpoint mode.

142

P. Bob, J. l’ignes / Local accuracy

Since the rules of algebra, such as the associativitv of alqebraic addition. are not satisfied in the arithmetic of the computer, there is no one data-processing procedure image of the algebraic procedure but rather a set S of C dataprocessing procedures images. COPis thgPtotal number of combinations corresponding to all possible permutations of the permutable operations in the algebraic procedure. We have Card S = C 2.2 Perturbation

given by:

2.3 Permutation-Perturbation -

Method

By applying the Perturbation Method to each of the PI imaaes of S. we obtain a set of 2 dataprocessing-results-which are the images of the single result r of the algebraic procedure. Each element of $ just as legitimately represents the algebraic result r. The cardinal of !$ is defined by: = CopZq

.

(4)

As is shown in L121, from set 2 it is easy to compute the number C. of exact si nificant deciIf R. is mal figures of each 'element of J f? any element of fi , the corresponding number Ci is defined by:

/Ri I

R+tp &)I=

P,rr c (R-tp l&9

1 - &

(7)

in which tp is the value of the Student bution for (N-l) degrees of freedom.

distri-

For N = 3 and p = 5, tp = 4.3 (Student's Table). The confidence interval is then defined by:

with

Let us consider one of the data-processing images PI of the algebraic procedure, and let us perfor it on a computer. Since each dataprocessing arithmetic operator generates an error (by lack for a chopping arithmetic), we must consider that for each arithmetic operation on the computer there are two results, one by lack and the other by excess, which represent, as legitimately one as the other, the exact algebraic results. Therefore if the PI contains q elementary operations (assignment, arithmetic operations, functions), this procedure does not supply a single result R but a set of Zq results, each of which just as legitimately represent the algebraic result r.

with

transform

rcI op'

Method

Card 2

irl the Fourier

C

Ic = i:-2.48

6, ?+2.48

61.

The number of exact significant Ce is equal to: R

ce = log10

with

5

(8) decimal

figures

g

= 2.48 6

R then Ce = log10 6 - 0.69624.

(9)

We thus have: C, = C - 0.69624

.

(10)

The generating of three elements R. depends on the alaorithm used. These elements'are obtained by applying the Perturbation F&hod to an image of the data-processing procedure. In the following section, we will apply this method to the Cooley and Tukey algorithm FFT :21. 3.

ESTIMATING

LOCAL ACCURACY

IN FFl

Review of the Cooley-Tukey

3.1

Algorithm

The algorithm published in 1965 by Cooley and Tukey and known by the term FFT (Fast Fourier Transform) can be used to obtain the transform of a sampled function g t CN in the form of an equally sampled function G E(I~~: Gj =

N-l C e k=o

-2i& N gi , for j=O,l,

. . . . N-l

(11)

‘i = log10 q

in which

ci = d(Ri-

and g' is the vector having the samples of g as components.

Ti)' + 6'

in which 3 9

R is the mean value of the elements

elements

of 3.

(5) of

6 is the standard deviation of the

Reciprocally,

the inverse

-2injk 1 N-l N gi=w,C e

transform

is given by:

, for k=O,l,..., N-l.

(12)

J=O

kaille 171 has shown that the best estimate of R. is given by its mean R. Therefore E. = 6, ahd the exact number of significant decimal figures is equal to:

c=

i = m

Ti .

log10 1

Since the R.s can be considered as a Gaussian random variable with mean 3 and standard deviation d i 71, the interval of confidence at p% of the exact result r defined by Students's law is

Below we will give the FFT algorithm for N = rI.r2 in which rI and r2 are two integers. Then let the indices j=j,r,

in eq.(ll)

t jo, j,=O,l,...,rI-1,

k=kIr2 t k,, ko=0.1,...,r2-1,

be expressed:

j,=O,l,..., r 2-1, kI=O.l,...,rI-I. (13)

Eq. (11) gives

:

P. &is.

J

lipws

/ I,ocal

accurac)’

in the I+‘ourier

143

transform

errors and to the computing errors defined in the preceding section. jklr2

G. Jl'jo

I W

in the Case of 3.2.1. Estimating Local Accuracy _--_-___------_______ExaSt-Da~a__--___-______ -_-___-___

jko

(14) with W = e

-2in N

= cos %

2ll - i sinT.

As Wrlr2 = WN= 1, since:

j k r ,WOI2

Wjk?2

>

('5)

inner

the sum, over kl, depends solelylon k, and can be defined as a new array g .

j. and

r -1 1 joklr2 1 (16: 'jo,ko = k $, g:l,ko ' 1 and after having carried out the external sum, the result can then be written as an

2 gjo,jl

r -1 w(j lrl+jo)ko 2 1 = k :, 'jo,ko

('7)

0

;;S;?yParing

eq.(14) 0

and eq.(17),

we obtain

L

G. Jl.jo = gjo,j,

the

Implementing the Perturbation Applied toTFT

Method

As was said in the preceding section, to estimate the local accuracy of the DFT it suffices to perturb a data-processing image of the FFT algorithm. We are going to deal with two possible use of the FFT algorithm.

a)

b)

In this case, perturbation consists of randomly last order bit adding a zero or a one to the of the mantissa of any result supplied by each data-processing operation (assignment, addition, subtraction, division, function). The CooleyTukey algorithm uses both operators working with real values (computation of sin 2n/N and cos 2rr/N contained in eq.(l4)) and operators working with complex values (inner products contained in eq.(14)). Perturbations are done with two distinct which perturbs functions, one called PEPERR a real value X c B;(, and the other PEPERC(Z) which perturbs a complex value Z t 63 . Each of these two functions uses the function XV(Y), Y Em whose role is to randomly add 0 or 1 to bit of the mantissa of Y. This the last order function XV(Y) is specific to the computer used and depends on the coding of the numerical values in computer. Function XV(Y) described here can be used only with CDC computers of the 6000 and 7000 series. 3.2.2. Estirnating_AErurac~_~~_~~~_~~~~_~~ -------E@i?riieiitaT Data ontalnlng ------Errors __ _-______-___-___-__-----

In this case the experimental data g are known with a relative accuracy A which is sssumed to be constant. Therefore the errors of data A must be associated with the propagation of computing errors. In these conditions, computer data gm are defined by:

(18)

’

If r = r = 2 the results are in the reverse orde;. T 6. 1s algorithm can be extended to any number N which is a product of any number of integers. 3.2.

such as were

cases of

The initial sampled function g. for which the transform G is sought after, is known exactly, i.e. samples g are affected only by the error resulting from the assignment operator. In this way, the error on each value of G is due solely to the computing error generated by the limited-accuracy arithmetic of the computer, The initial sampled function g, for which the transform G is sought after, is known experimentally, and hence each sample g is assigned a relative error A. In this way the error for each value of G is due to the data

gm

= g,(i t 0A)

0 being a random number uniformly between the values -1 and +l.

(19) distributed

As was mentioned in section 2, local accuracy for each element G will be obtained from three runs on the computer of the Cooley-Tukey algo run being done with differithm, with each rent random perturbations. The number of exact significant decimal figures is computed by the ESN subroutine for which the listing is given in the Appendix.The ESN subroutine uses the formulas described in section 2. N.B.: The FFTPER subroutine (in the Appendix) thereal and imaginary parts of the DtT, but the modules and phases can'be computed very easily. 4.

LISTING

FOR FORTRAN

gives

SUBROUTINES

The algorithm which implements our method requires several subroutines presented in the Appendix:

. FFTPER perturbs FFT by using the functions

PEPERR and PEPERC which respectively perturb a real or a complex value. This perturbation .is performed by the function XV which randomly adds a zero or a one to the last order bit of the mantissa of the value considered. ESN computes

the exact number of significant

144

P. Bois, J. C’ignes / I.ocal

decimal

figures on each component

A main program

is

of the FFT.

also given in the Appendix.

5. EXAMPLES Two applications of the routines are given here, one concerning exact initial data and the other concerning initial data containing errors. Table 1 Results obtained 7024&s

from a complex

function with

accuracy

in the Fourier

5.1. Example

Since too large a table is required for presenting 1024 values of the Fourier transform, we will present only a certain number of them chosen at random.

Imaginary

Part

- 22.540676373160984080

27.362054314926023855 1

Using Exact Data

The initial data of the function for which we want to compute the Fourier transform are randomly generated complex numbers which are uniformly distributed between the values -1 and tl. Therefore the only errors affecting the Fourier transform are computing errors. We consider a function with 1024 complex samples.

Real Part

j

transform

l_l-

27.3620543149268

c = 14

0.30436576039091963768 31 0.304365760389297

c = 11

- 22.5406763731612 -

6.3199102026420999694

-

6.31991020264309

0.67643823378508592443 79

0.676438233781973

c = 13

c = 12

26.370038791154481990 c = 11

26.3700387911574

- 9.3062710329200841270

c = 13

10.360918417726311801

275 - 9.30627103292579

c : 12

39.225546346328048287 448

39.2255463463464

c = 12

20.315072867168419363 c = 13

1.0384628470167623661 628

~ - 0.489288362417199

c = 12

8.58485790920599

c = 12

- 20.3008439988596

c = 13

6.87236078082236

c = 12

c = 11

2.12679639642626

0.0814827802247082

c = 14

- 0.0292570411406814

c = 11

0.71015027929282147140

-

4.7055515901892891596

-

4.70555159019762

c=9

c=9

- 0.27666192684760097446 - 0.276661926841181

c = 11

0.86511479870731649062

962 0.710150279280430

c = 12

- 0.029257041129125689884

- 26.330955013083845979 - 26.3309550131019

c = 14

0.081482780233098124808

- 14.479042895715677390 - 14.4790428957155

25.3998573917352 2.1267963964204776780

6.8723607808165800158 799

878

-

25.399857391735793637

- 0.489288362417205760250

854

8.5848579092036735106

1.03846284701817

749

-

c = 13

- 20.300843998858786857

526 20.3150728671692

10.3609184177295

c = 10

0.865114798709573

c = 11

-15.543975041433886085

1024 c = 11

-15.5439750414376

c = 13

P. Bois, J. ViKnes J Local accuracy

Table 1 shows the agreement between the number of significant decimal figures C determined by our method and the number of exact decimal figures obtained by a run in double precision. In addition we checked eq.(8) by applying our method with approximately 100 functions made up of 128, 256, 512 and 1024 complex samples. As shown in Table 2, we found that, with a probability of 95% the true value of DFT (obtained by using the Cooley-Tukey algorithm with double precision) is always within the interval Ic defined by eq.(8). Table 2 NCS

N1

0 1C

%

N2 c Ic

128

6

250

256

10

502

2.0

512

29

483

2.8

1024

46

978

2.2

NCS: Number

of complex

in the Fourier

14s

transform

Table 3 gives only a certain number of values of the Fourier Transform with their number of exact significant decimal figures. 6.

CONCLUSION

The Perturbation Method applied to the FFT has enabled us to work out an effective software for determining the local accuracy of the DFT. This software can be applied to computing the DFT of a function using exact data. In this case the error in the real and imaginary parts of the DFT solely depends on the propagation of computing errors resulting from the limited-accuracy arithmetic of the computer. This software can also be applied for calculating the DFT of a function using experimental samples containing a known error. In this case the errors in the real and imaginary parts of the DFT depends on data errors and computer errors. This software is very easy to implement. It has been used with success in a great many cases and has never failed.

2.3

REFERENCES

Cl1

samples

R. Alt, Error Propagation Tra;;o;ms

in Fourier

Math. Comp. Sim. XX (1978)

: Number of real parts and imaginary parts {I,

Nl N2

: Number of real parts and imaginary parts

NE: The results given in this table correspond 5 the mean values of the approximately 100 DFTs computed. It can be seen that fewer than 3% of the points on the transform are outside of the interval I thus illustrating the perfect agreement befieen theoretical results and the results provided by the proposed method. 5.2.

Example

Using Data Containing

Errors

The initial data are 256 samples from a seismic trace taken from a survey made in the Paris Basin. The accurac of the data is taken as being equal to lo-Y . The errors affecting the Fourier Transform are both computing and data errors. Table 3 Results obtained from 256 real samples the seismic trace J

Real Part

C

Imaginary

10

0.125

3

- 2.744

4

2.605

65

0.004

1

0.242

0

0.009

2

- 0.086

120 154 I87 204 213 240

- 0: 0.016 - 0.01 0.055 0.16 - 2.944

254 0.149 * unsignificant value

c21

J.W. Cooley and J.W. Tukey, An Algorithm for the Machine Calculation of Complex Series, Math. Comput. vol. 19,(1965) 297-301.

131

W.M. Gentleman and G. Sande, Fast Fourier Transforms - for Fun and Profit, Proc. Fall Joint Computer Conf., (1966) 563-578.

c41

D.V. James, Quantization Errors in the Fast Fourier Transforms, IEEE Trans., Acoust., Speech and Signal Proc. (1975) VOI. ASSP - 23 IWO. 3.

I51

T. Kaneko and 6. Liu, Accumulation of Roundoff Error in Fast Fourier Transforms, J. Ass. Cornput. Mach., vol. 17 (1970) 637-654.

[61

M. La Porte and J. Vignes, Evaluation de l'incertitude sur la solution d'un systeme lineaire, Numer. Math., No. 24, (1965) 39-47.

c71

M. Maille, Methodes d'evaluation de la precision d'une mesure.ou d'un calcul numerique, LITP Report, Institut de Programmation, Universite Pierre et Marie Curie, (1979).

181

A.V. Oppenheim and C.J. Weinstein, Effects of Finite Reqister Lenqth in Diqital Filtering and Fast Fourier Transform, Proc. IEEE (invited paper), vol. 60, (1972) 957-976.

191

M. Sunoaramurthv

from -

Part

21

l1c

- 0.19

1

- 0.027

2

- 0.266

2

- 0.006

4

- 2.628

4

3

- 0.090

2

and Reddv V. Umaoathi. Some results in-Fixed-Point Fast'Fourier Transform Error Analysis, JEEE Trans. Cornput.,.C 26, No. 3, (1977).

1101

Tran-Thong and B. Liu, Fixed-Point Fast Fourier Transform Error Analysis, IEEE Trans. Acoust., Speech and Signal Proc., Vol. ASSP-24, No. 6, (1976) 563-573.

1111

Tran-Thong and B. Liu, Accumulation of Roundoff Errors in Floatina Point FFT. IEEE Trans. Circuits and S&t., Vol. CAS24, No. 3, (1977) 132-143.

I121

J. Vignes, New Methods for Evaluating the Validity of the Results of Mathematical Computations, Math. Comp. Sim., XX No. 4, (1978) 227-249.

1131

C.J. Weinstein, Quantization Effects in Digital Filters, Mass. Inst. Technol. Lincoln Lab. Lexington, Tech. Rep. 486, ASTIA Dot. DDCAD-70 6862, (1969).

I141

C.J. Weinstein, Roundoff Noise in Floating Point Fast Fourier Transform Computation,IEEE Trans. Audio Electroacoust., vol. AU17, No. 3, (1969) 209-215.

1151

P.D. Welch, A Fixed-Point Fast Fourier Transform Error Analysis, IEEE Trans. Audio Electroacoust. (special issue on Fast Fourier Transform), Vol. AU-17, (1969) 151-157.

P. Bois. J. l’ijyes

C C C

C C C C C C

C C C C C C : C C

acc!lracJ,

in the Fourier

tramform

PROGRAM TEST(INPUT,OUTPuT,TAPE5=INPUTtTAPE6=OUTPUT~ THIS MAIN PROGRAM CARRIES OUT THE COMPUTATION OF THE F F T AND GIVES THE EXACT NUMBER OF DECIMAL SIGNIFICANT FIGURES ON COMPONENTS OF THIS F F T IT USES THE SUBROUTINES FFTPER AND ESN FFTPER CARRIES OUT THE PERTURBATION METHOD ON COOLEY TUKEY ALGORITHM ESN COMPUTES THE EXACT NUMBER OF DECIMAL SIGNIFICANT DIGITS IN THE F F T INPUTS ARE NX,NExP,IK,EPS,DR,DI NX IS THE TOTAL NUMBER OF SAMPLES OF THE FUNCTDN TO BE TRANSFORMED NEXP IS DEFINED Ey NX=2*“NEXP IK*O:THE FUNCTION TO BE TRANSFORMED IS COMPLEX IK=OITHE FUNCTION To BE TRANSFORMED IS REAL IF THE DATA DO NOT CONTAIN ERROR EPS=O, EP%=X, IF THE DATA CONTAIN A RELATIVE ERROR OF x PER CENT DR-IS THE ARRAY OF THE REAL PART OF VALUES OF THE FUNCTION To BE TRANSFORMED DI-IS THE ARRAY OF THE IMAGINARY PARTS OF VALUES OF THE FUNCTION TO BE TRANSFORMED OUTPUTS ARE CXM,CR,CI CXM IS THE ARRAY OF THE COMPLEX VALUES OF THE F F T CR IS THE ARR4Y OF THE EXACT NUMBER OF DECIMAL SIGNIFICANT FIGURES IN THE REAL PARTS OF THE F F T CI IS THE ARRAY OF THE EXACT NUMBER OF DECIMAL SIGNIFICANT FIGURES IN THE IMAGINARY PARTS OF THE F F T DIMENSION DR(lO24),DI(lO24) COMPLEX CX(1024,41,A~1024),8(1024),CXM(10241 COMMON AA,HB DIMENSION ~~(10241,C1~1024)tC~1024) lOO,Nx,NEXPtIK,EPS 1 READ IF(NX.EQ.01 STOP READ lOl,(DR(I)vi~l,NX) IF(IK,NE.O)GO TO 2 READ lOlr(DI(I),I=lrNx) GO TO 4 00 3 I=l,NX DI(I)=O, DO 0 J’lr3 DO 5 I=leNX A(I)=CMPLX(DR(I),DI(I)) AA=O. BB=1,99999999999999 DO 6 I=l,NX 6 B(I)=A(I) CALL FFTPER(B,NEXP,EPS) DO 7 I=l,NX 7 CX(ItJ)=B(I) B CONTINUE CALL ESN(CX.NX,~,CR,CI,CXM) DO 9 I=l,NX 9 PRINT 200,CXM(I),CR(I1rCI(I) FORMAT(314rE6,3) 100 FORMAT (12F6.5) 101 FO~~AT(SX~~25.l5~23XIE25.l5~lOX,2FlO.O) 200 GO TO 1 END

E :C C

/ Local

147

148

P. Bois, .I. Vignes / Local

SUBROUTINE

accuracv

in the Fourier

transform

FFTPER(A,M,EPS)

C C

C C

C

THIS SUBROUTI,NE CARRIES OUT THE PERTUEATION IT USES THE ~ONCTIONS PEPERR AND PEPERC PEpERR PERTURBS A REAL VALUE PEpERC PERTUR3S A COMPLEX VALUE

METHOD

ON

CDOLEY

TUKEY

ALGORITHM

C

C C

INPUT

IS

A A

IS THE VECTOR OF COMPLEX SAMPLES OF THE DISCRETE FUNCTION “IUST BE TRANSFORMED BY F F T M IS DEFINED BY N=2**M,N IS THE TOTAL NUMBER OF SAMPLES EPS IS THE ACCURACY OF THE DATA (FOR EXACT DATA EPS=O.)

C C

C C

WHICH

C

E C

OUTPUTS

ARE A

C C

1 2 3

4 5

6

7 8

A IS THE I,NITIAL

VECTOR OF THE COMPONENTS SAMPLES ARE OF COURSE

COMPLEX A(l),U*W,T,PEPERC,CMPLX PI=PEPERR(4.*ATAN(l.O)) N=2**M Nv2=N/2 NMl=N-1 J-l IF(EPS.EQ.O.lGO TO 2 Do 1 I’lrN A(I)=A(I)*(l.*HASARD(-l.~l,)~EPS) DO 3 I=l,N A(I)=PEPERCfA(I)) 00 6 I=lrNMl IF(I.GE.J!!O TO 4 T3.A(J) A(J)=A(I) A(I)=T K=NV2 IF(K,GE,J)GO TO 6 J=J-K K-K/2 GO TO 5 J=J+K 00 8 L’lvM LE=2”*L LEl=LE/Z U=(l.*O,) ANG=PI/FLOAT(LEl) Gl=PEPERR(COS(ANG)) GZ=PEPERR(sIN(ANG)) W=cMPLX (Gl gG2) DO 8 J=l,LEl DO 7 I=JvN*LE Ip=I*LEl T=PEPERC(A(IP)*U) A(IP)=PEPE:C(A(I)-T) AtI)=PEPERc~A(I)+T) U=PEPERC (U’h) RETURN END

OF F DELETED)

F

T

(CONSEQUENTLY

THE

P. Bois. .I Vignes / Local

C C C

C

C C C C

accuracy

in the Fourier

COMPLEX FUNCTION PEPERC(PC1 THIS FUNCTION PERTURBS A COMPLEX VALUE PC USING XV IT USES A FUNCTION HASARD GENERATING A RANDOM NUMBER WITH A CONSTANT PROBABILITY DENSITY COMPLEX PC COMMON AA,!0 U=REAL (PC) V=AIMAG(PC) IPILE=HASARD(AA,BBI IF(IPILE.E~IOl GO TO 1 1u=u U2=IU IF~(U2.EQ.U~.OR.(U.EQ.O,OOl~ GO TO 1 PRPC=Xv~U,SIGN(FLOATfIPILE)~U)l GO TO 2 PRPC’U IPILE=HASARDfAA,BBl IF(IPILE.EQ.01 GO TO 3 IV=V Y2=IV IF~(V2.EQ.v)rOR.~V.EQ,O.O01) GD TO 3 PIPC=XV(V,~IGN(FLOAT(IPILE)tV)l GO TO 4 PIPCIV PEPERC=CMPLX(PRPC,PIPCl RETURN END FUNCTION Pt;PERR(PR) THIS FUNCTION PERTURBS A REAL VALUE IT USES A FUNCTION HASARD GENERATING WITH A CONSTA:;T PRDBAGILITY CENSiTY COMMON AA,HB IPILE=HASARD(AA,BR) IF(IPILE.Efj.01 GO TO 1 IPR=PR PRS=IPR IF(PRS.EQ.FR) GO TO 1 PEPERR=XV(PRISIGN(FLOAT(IPILE),PR)) GO TO 2 1 PEPERR=PR 2 RETURN EYD

RR USING A RANDOM

FUNCTION IN AN

XV FUNCTION NUMBER IN

FUNCTION XvfXtSl THIS FUNCTION RANDOMELY PERTURBS THE VALUE X DEPENDING IF S=O,NO PERTURBATION OF X IF S=+lePERTURBATiON OF POSITIVE X IF S=-I~PERTURBATION OF NEGATIVE X ~RE~Z~=OR~OO00400000OOOOOOOOOOB,AND~ARS~Z~,7777OOooOOOooOOOoooofl~~ DER~Z)=OR~O000000000OOOOOOOOOlB~AND~ABS~Z~,7777OOo~OOOOOOOOOOooB~) IF(X.EQ.O.lGO TO 1 Xv=X+SiGN(DER(X) rS) IFfPREtXl .gT.PRE(XV) lXv=Xv-SIGNtDERfXv) ,S) RETURN 1 ;;;;;;N”““O.‘,Sl END

149

transform

AN

ON

INTERVAL

AArBB

INTERVAL

AA,BB

THE

VALUE

OF

S

P. Bois. J. Vignes / Local accuracy in the Fourier trarlsform

: C C C C C C

c” C C

C

SUBROUTINE ESN(CX,N,M,CR,CI,CXM) THIS SUBROUT1,NE COMPUTES THE NUMBER OF DECIMAL SIGNIFICANT DIGITS ON THE FFT INPUTS ARE CX,N+M CX IS THE ARRAY OF THE F F T OBTAINED BY THREE SUCCESSIVE RUNS OF THE FFTPER SUBROUTINE N IS THE TOTAL NUMBER OF SAMPLES IN THE F F T M 1s THE NUMBER OF RUNS(GENEHALY M=3) OUTPUTS ARE CR,CI,CXM CR IS THE NUYBER OF DECIMAL SIGNIFICANT DIGITS IN THE REAL PART OF F F T CI IS THE NUMBER OF DECIMAL SIGNIFICANT DIGITS IN THE IMAGINARY PART OF F F T CXM IS THE MEAN VALUE OF CX DIMENSION CRf1024)~CI(1024),SIGR(lO24)rSIGI(lO24)rGR(lO24~ DIMENSION GI(1024) COMPLEX CX~1024,3),CXY(10241 ,SoM(1024) COMPUTATION OF CXM(I)eMEAN VALUE OF CX(II DO 1 I=lrN SOM(I)=fO.*O.) DO 1 J=lrM

1 SOM(I)=SOM(I)+CX(ItJ) DO 2 I=l,N 2 CXM(I)=SOM(Il/FLOAT(M) C

COMPUTATION

OF THE STANDART DEVIATIONS OF REAL AND IMAGINARY PARTS OF CX(1) DO 3 I=l,N SIGR(I)=O. SIGI(I)=O. DO 3 J=l,M SIGR(II=SI~R(I~*(REAL[CX(I,JII-REAL(CXM(I)l1*"2 3 SIGI(I~=SI~I(I~*(AIMAG~CX(I,J)~~AI~~AG(CXM~I~~l~~2 DO 4 I=lvN SIGR(I)=SQRT(SIGR(I)/FLOAT(M-11) 4 SIGI(l)=SQRT~SIGJ(I)/FLOATO) COMPCTATI3N OF CR(I) DO IO

I=lvN

IF(SIGRtI)eEQ.O.)GO TO 6 GR;I)=ABS(REAL(CXM(I,/SIGR(I)))

ANO,SIGR(I).NE.O,)GO TO 5 IF(GR(I).E'j.O.. AND.SIGRtI).EQ,O.)GO TO 6 IF(GR(I).E!j.O.. CR(I)=ALOGlO(GR(I)) GO TO 7 5 CR(I)=o, GO TO 7 6 CR(I)=lS. COMPUTATION OF CI(1) 7 IFISIGI(I).EP.O.)GO TO 9 GI(I)=ABS(AIMAG(CXM(Ij/SIGI(I))) IF(GI(I),EQ.O ..AND.SIGI(I),NE.O.~GOTO 8 IFIGI(I~.E~.O.. AND.SIGI(I).EQ.O,~GO TO 9 CI(I)=ALOGlO(GIfI)) GO TO 10 8 cI(I)=o. GO TO 10 9 CI(I)=15. 10 CONTINUE RETURN END