NBSIR 73-303

FAST FOURIER TRANSFORM IMPLEMENTATION FOR THE

CALCULATION OF NETWORK FREQUENCY DOMAIN TRANSFER

FUNCTIONS FROM TIME DOMAIN WAVEFORMS

William L. Gans and N. S. Nahman

Electromagnetics Division

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado 80302

December 1972

Final Report

U.S. DEPARTMENT OF COMMERCE, Frederick Bo Dent, Secretary

NATIONAL BUREAU OF STANDARDS, Richard W. Roberts. Director

Prepared for DOD Calibration Coordination Group

72-65

CONTENTS

Chapter Page

1 DESCRIPTION OF FAST FOURIER TRANSFORM 2

1. 1 Definition of Discrete Fourier Transform 2

1. 2 Some Properties of the DFT 3

1. 3 Definition of the Fast Fourier Transform 4

2 . OPERATING SYSTEM SOFTWARE 8

2. 1 FFT Subroutines 8

2. 2 Input Data Programs 19

2. 3 Output Data Programs 22

3 EXPERIMENTAL CONCLUSIONS 25

3. 1 Problems Encountered Using the FFT 25

3. 2 Time Domain Measurement Example 27

3. 3 Future Developmental Considerations 28

REFERENCES 30

ILLUSTRATIONS 31

APPENDIX A 42

iii
LIST OF ILLUSTRATIONS

Number Page

2-1 Graphical interpretation of waveform averaging 31

process.

2-2 Typical noisy time domain waveform and various 32

results of waveform averaging.

2-3 Time domain waveform, both before and after 34

being made periodic.

2- 4 Typical displays of magnitude versus frequency 35

spectra for time domain waveforms.

3- 1 Two possible magnitude spectra of an ideal 36

rectangular pulse.

3-2 Magnitude versus frequency spectrum of an ideal 37

rectangular pulse with sin x variation removed.

3-3 Block diagram of time domain measurement system. 38

3-4 Input and output waveform from 10 MHz low pass 39

filter experiment.

3-5 Input and output spectra for 10 MHz low pass filter. 40

3-6 Plot of algebraic difference between input and output 41

spectra for 10 MHz low pass filter.

A-1 Eight-point representation of a time domain 47

waveform.

A-2 Two four-point representations of the same time 48

domain waveform.

A-3 Signal flow graph of FFT computation of an eight- 49

point time domain waveform.

IV
FAST FOURIER TRANSFORM IMPLEMENTATION

FOR THE CALCULATION OF NETWORK FREQUENCY

DOMAIN TRANSFER FUNCTIONS FROM TIME DOMAIN WAVEFORMS

by

William L. Cans and N. S. Nahman

ABSTRACT

This report is concerned with the software applications of the fast

Fourier transform algorithm to the relationship between time domain

waveforms and frequency domain spectra. The first chapter is devoted

to a description of the discrete Fourier transform and the fast Fourier

transform. Chapter 2 contains the text and a brief description of all

FORTRAN n programs utilized in connection with this work. All computation was performed on the in-house time share computing system in

the NBS facilities, Boulder, Colorado. In Chapter 3, problems en- countered using the fast Fourier transform algorithm are discussed, an

example of a time domain to frequency domain calculation is presented,

and future developmental considerations are mentioned. In addition

Appendix A contains a detailed example aimed at disclosing the inner

mechanisms of the fast Fourier transform algorithm.

Key Words: Discrete Fourier transform; Fast Fourier transform;

Frequency spectra, discrete; Network transfer function; Time domain

waveform; Transfer function.
1. DESCRIPTION OF FAST FOURIER TRANSFORM

1. 1 Definition of DFT. For many years, the Fourier integral transform, the Fourier

series, and the Laplace transform have been available to engineers for

the purposes of time domain- frequency domain correlations. More

recently, the sampling theorem and the discrete Fourier Transform

(DFT) have evolved as useful tools for handling sampled-data bandlimited time waveforms and frequency spectra. Thus, while the Fourier

series and the Fourier and Laplace transforms are well suited to operations involving continuous, or piece-wise continuous functions, the DFT

allows operations to be performed directly upon a series of discrete

data samples representing a continuous time fvinction.

Most recently, (1965) a method has been reported [ 1] whereby

the Fourier coefficients of the DFT may be machine calculated in a fraction of the time previously required. This algorithm is commonly

referred to as the Fast Fourier Transform, or FFT.

In order to discuss the FFT it is first necessary to define the DFT.

Definitions vary in the literature , however, this paper will follow the

definitions and notations used in [ 2] . The DFT is defined by

where Ar is the rth coefficient of the DFT and X^^ is the k sample

of the time series being transformed. j =/^l and there are N

samples in the total time series. Thus, a time waveform consisting of

N sequential samples will yield a transform consisting of N frequency

domain coefficients. For convenience, let

N-1

Ar r = 0,1,2. .., N-1 (1.1)

th

W

-27r j/N

(1.2)

2
then, in simplified form,

N-1

' k=o

As in the case of the Fourier integral pair there exists an inverse

transform of the DFT. This transform is written

N-1

X. = ^m^^W"''^ i = 0,l,2..., N-1 (1.4)

i N

^

r=o

and is called the inverse discrete Fourier transform or IDFT.

1. 2 Some Properties of the DFT

Both the Ar's and X^'s may be complex numbers; the Ar's, in

general, are almost always complex, but for a real time series wave - form, the X 's will be real.

i

In addition there exists a convolution relationship between the DFT

and IDFT, The IDFT of the product of two DFT's is the periodic mean

convolution of the two time series waveforms of the DFT's.

Likewise, most other properties of the Fourier integral transform may be shown to apply to the DFT. Asa further example, the

DFT of the sum of two time series waveforms is the s\im of the DFT's

of the two waveforms.

Thus, with few exceptions, ftie DFT lends itself to time- frequency

operations on discrete sampled-data functions with the same power and

versatility as the Fourier integral transform brings to operations on

continuous functions. One exception to be kept in mind, however, is that

data obtained by use of the DFT is valid only at the points in time or

frequency where samples were taken. Interpolation between the Ar's or

X 's is neither valid nor wise.

i

3
1. 3 Definition of the Fast Fourier Transform

The fast Fourier transform is an algorithm which allows the

machine computation of a set of DFT coefficients to be accomplished

much faster and more efficiently. Previous straightforward

2 methods require approximately N arithmetic operations to

transform a time series waveform consisting of N = 2^ samples. By

use of the FFT, the same results can be obtained with only 2N log N

operations. A 1 024 point waveform, for example, would require about

2

(1024) or 1, 048, 576 arithmetic operations to obtain its 1024 DFT

coefficients by conventional means. Using the FFT, these same coefficients may be obtained after only (2) (1024X10) or 20, 480 operations.

Thus, the time and effort required to transform a 1024 point time series

using the FFT is about l/50th of that required using straight-forward

techniques. With larger time series, the savings increase.

Although many variations of the FFT algorithm now exist, only

Cooley and Ttikey's decimation in time algorithm shall be described. [ 1]

Suppose we wish to calculate the DFT of a time series waveform consisting of N samples. For simplicity, let N =

2''^ where n is an integer.

The DFT of this waveform is defined as:

A^ =5]^X^e-2'^j"^/^ r = 0, 1. 2..., N-1 (1.5)

^ k=o

If the original time sequence is now split into two sequences such that

the even-nxrmbered X^^'s become one sequence and the odd-nximbered

Xj^'s the other sequence, then we may define,

= X,, k 2k

N

and k = 0, 1, 2, . . . , -J

-1

4
N

Both Y and Z, are now -— -point time series and they each have

k k 2

DFT's as follows,

and

B = :=oN

-4 ujrk/N

C Z, e

r 4^ k

k=o

-4 7rjrk/N

N

r = 0, 1, 2, . . . , — -1

(1.7)

N

r = 0, 1, 2, . . - -1

Now, writing the original DFT in terms of these two new ones:

r = 0, 1, 2, . . . , N-1

or

k=o

k k

^ Y e-4 7rjrk/N

^ ^-2 TTjr/N ^-4 TTjrk/N

r Z-» k

^

k

^ k=o

using (1.7),

(1.8)

(1.9)

A =B +e-^^j^/^C

r r N

r = 0, 1, 2. . . , — - 1 (1. 10)

N

For values of r greater than -— , the DFT's of B and C repeat

2 r r N N

periodically their values for r < -r . Therefore, substituting r + -r for r in ( 1. 10) yields

5
= B +e N '• Z-' C

r r 0 ^r< - (1. 11)

or, since e

-J 77 1 .

- 2 ffjr/N _

e Li N 0^r< - (1. 12)

Therefore (1. 2). (1. 10) and (1. 12) may be combined to yield in terms of

B and C : r r A = B + W C

r r r

(1. 13)

A(r+^1 = B - W^C 0

\ 2; r N

(1. 14)

Thus, the DFT of X may be obtained by computing two smaller DFT's

N

of length — . It should be noted that

(1. 15)

If this decimation process is continued, B and C would each reduce to

r r N

two — -point transforms. The resulting four transforms would reduce

N

to eight — -point transforms, and eventually, the problem would be

^ N n reduced to solving —- two-point transforms. In particular, for a 2 -

point waveforn, n such reductions will be necessary to reduce the

problem to one of solving only two-point transforms.

The result of this reduction process is a problem consisting only

of complex multiplications and additions. Furthermore, three-fourths

of the multiplications may be eliminated because they either involve

multiplication by unity, or they are red\indant multiplications. Thus,

6
in general, an N-point transform will consist of N log N coraplex

1 additions and — N log^ N complex multiplications. See Appendix A

for a detailed example of an eight-point transform.

7
2. OPERATING SYSTEM SOFTWARE

2. 1 FFT Subroutines

Four subroutines are used to compute the DFT and IDFT

14

of a time waveform consisting of up to 2 , or 16, 384 sequential

data points. These subroutines, as well as all other programs in this

report, are written in FORTRAN II and are implemented on the in-hou

time share computing system at NBS, Boulder, Colorado.

In order to compute a "forward" DFT, or to go from the time

domain to the frequency domain, the following call sequence of subroutines is used:

CALL REORDER (A, B, MM)

CALL CFFTRC (A, B, MM, . 5>:^SC, NX)

CALL REALTRAN (A, B, MM, NX, INV).

The notation used in these call statements relates to the DFT as

follows:

N-1

S = A +jB = SCV X w|^^^^^^

r r r k N

k=o

where NX = ±1 (indicates exponent sign of W^^)

SC = real scaling factor

INV = ± 1 (indicates forward or reverse transform)

^(IvIM)+l , , r 2 = length of transform

As an example, for the case of a 1, 024 point DFT the specific call

statements would be:

8
CALL REORDER (A, B, 9)

CALL CFFTRC (A, B, 9, .S^H., 1)

CALL REALTRAN (A, B, 9, 1, 1)

The se quence of call statements for the IDFT to return from the

frequency domain to the original time domain waveform are:

CALL REALTRAN (A, B, MM, NX, INV)

CALL CFFTS (A, B, MM, l./(N-SC), NX)

CALL REORDER (A, B, MM)

For the previous example, the sequence of call statements for the IDFT

relating to the 1024 point original time series would be:

CALL REALTRAN (A, B, 9, -1, -1)

CALL CFFTS (A, B, 9, 1./1024. , -1)

CALL REORDER (A, B, 9)

For the 1024 point case, SC is assumed to be unity; NX and INV are +1

for the DFT and -1 for the IDFT.

A source of possible notational confusion exists with the two

arrays A and B since they are used to define both input and output data.

For the DFT, or time domain to frequency domain case, A and B contain the N values of time domain data points as follows:

N

0, 1, 2. . . , j-l (2. 2)

0, 1, 2. . . , j-1 (2. 3)

Thus the first half of the time domain waveform is contained in A

and the second half in B.

9
Once the DFT subroutines have run, the contents of A and B are

replaced with the output data as follows:

S = A +jB r = 0, 1. 2, --1 (2.4)

r r r where A and B are the real and imaginary parts, respectively, of

r r the coraplex frequency domain coefficients. For a 1024 point time domain waveform, then, the first 512 points would be sequentially stored

in A and the last 512 points in B. The resulting DFT solution would be

512 frequency coefficients each consisting of a real part. A, and an

imaginary part, B.

Therefore, for each frequency component.

Subroutine listings are given for REORDER CFFTS, CFFTRC, and

REALTRAN on pages 11-12, pages 13-14, pages 15-16, and page 17,

respectively.

These four subroutines are part of an FFT subroutine package

contained in the public user's subroutine library of the in- house time

share computing system at NBS, Boulder, Colorado. For additional

information concerning these subroutines see Ref. [5].

MAGNITUDE (r) (2. 5)

PHASE (r) (2.6)

10
01 MENS UN A(1^24)>d(liJ24)*LSr (lb)*LG(lb)*Sr(l:))

CJMMJN M*LC*ST

M = MM

CALL FFTC

533 CJNriNJE

IF(M)301 #173*331

3131 CJNriNUE

JB=M-1

K8=3

1=3

KJ=LC( JA)-l

0 J 23 KA = 1 >KJ*2

T = A<KA+l ) A(KA+1 ) = 3 (KA)

3(KA) = T

23 CJNTINJE

IF (M-1 )332> I 73>332

332 CONTINUE

LIM=M/2+l

33 KS = LC( JB+1 )+Ka

KJ = KS

JJ=LC( JA-J8)

KK = K3 + JJ

43 K = KK + JJ

53 r = A <KK + l ) A(KK+1 ) = A(KS+1

)

A(KS+1 ) = T

T = B(KK+1

)

B<KK + l ) = B (KS + l ) B<K3+1 ) = T

= KK \ KS = KS + 1 IF(KK-K)53*333*333

333 CJNTINJE

KK=KK+JJ

KS = KS + JJ

IF (KK-KU )40*334>334

334 CJNTINJE

IF(J8-LIM)93*93>335

11
306 CONTINUE

J8=JB-1

I s I-^l ^ LST< I) = J3

GO TJ 3^

90 IF(I)17a*170#306

306 CONTINUE

JB=LST(I)

I - I - I KB = KS

GJ TJ 30

170 RETURN

END

12
SJQ'RJJTIiME CFFrS (A*3*MM*SCALii:*iMc:XP)

DIM£t\JSI JiM A(163 84)*a(16384)# JC(l:>)>SNr<lb)

CJMMJN M>JC*SNT

M=MM

MA = MM

CALL FFTC

IF(M-23i<3)53'3*l 73>5^3

6 3 3 CJiMTINJE

iM=JC(M+l ) NH=iM/2

NQ=NH/2

SC = SCALE

IF(M-2 )301 *332*3a2

331 IF(M)123>123> 101

332 CJNTIiNiJE

EXPS = ISIGiMCl *N£XP)

NiM = N-1

K = 1 03 100 JA = 2*M

C£=SNT<MA)

I * CE

33

43

MA MA - I CD 2. * CE

SO -SNTCMA)

R = -2' * CD

CN 1 • CM 0.

S.M 0.

JJ 0

KK 1 SM ^EXPS

JSPAN = NH

JSPAiM/2

KS KK + JSP

RE s A(KK) - A(KS)

ACKK) = A(KK) + A(KS)

FIM « 8(KK) - a(KS)

8CKK) = 8(KK) + a<KS)

A(KS) = CN*R£ -SN*FIM

8<KS) = SN*R£ CN*FIM

KK = KK + NH

KS = KS + NH

13
R£ = A(KK) - A(KS)

A<KK) s A(KK) + A(KS)

FIM = aCKK) -8<KS>

^ B<KK) = 3(KK) + 8(KS)

8(KS) = SM * RE + CM * FIM

KK = KS NH

IF<KK-N)40*333#3a3

3 33 CONTINUE

50 KK s KK - NN

JJ * JJ + K

IF( JJ-NQ)304*90#93

304 CONTINUE

60 CD = R CN + CO

CN = CD CN

SO s R * CM + SD

CM * CM + SO

SNa-CM*EXPS

SM=CN*EXPS

GO TO 40

90 K s K + K

100 CONTINUE

101 DO 1 10 KK s 1 *N#2

KS = KK + I RE = A(KK) - A(KS)

ACKK) s (ACKK) + ACKS))

ACKS) = RE

FIM = 8CKK) - 8CKS)

8CKK) = C8CKK) + SCKS))

8CKS) = FIM

110 CONTINUE

120 IF CABS CSC -1 • )-l .E-1 3)13b*335>305

305 CJNTINUE

DO 130 Ja=l^N

A<J3) = SC * ACJ3)

BCJB) = SC * 8CJ3)

130 CONTINUE

135 RETURN

173 RETURN

END

14
"FFTS"

SJB^JtJTINE GFrT^C ( A * 3 * MM, SCALd * ) 131M£NSI JN A(l /)^4)>3(l/)34)* JlJCl is ) > S ( 1 3 > GJMMJ.vJ M*JD*S

M = MM

CALL FFTC

IF(M-2^^'4)bi3<3*17^>oa^

l^J = JL)(M^l ) K = iM/4 = K

iMiM = iNj-l

J3PAN = 1 SC = SCALE

IF CASS (SC- 1 . )-l .£-1 ^)7*30l *3J1

3(dl CONTINUE

6 OJ i> JC = 1 *iM A(JC) = SC * A(JC)

8<JC) = SC * a(jc)

5 CJi^JTIi^aE

7 IF (M)332* 170*302

332 CJNriNUE

OJ 13 KK = 1 >iM,2

KS = KK+1

>^£ = A(KK) - A(KS)

A(KK) = A(KK) + A(KS)

A(KS) = RE

FIM = 8(KK) - aCKS)

8(KK) = B(r<K) + 3(KS)

aCKS) = FIM

' 10 CJiMTIiMJE

IF(M-1 )303*170,303

303 CJNTIiMJE

EXPS = ISIGiMCl ,N£XP)

i>iP = 1 DJ 90 J3 = 2*M

SD = -S< Ja-1

)

CD = 2.* scja) * s(ja)

R = -2» * CD

CN = 1 . CM = 0.

SlM = 0»

JJ = 0

KK = I SM = +EXPS

15
12 JSPAi>JH = JSPAN

JSPAiM = JSPA.M + JS-^AiM

23 KS = KK + JSPAiNJ

RE = CiM A<KS) - * BCKS)

FIM = Si^ * A(KS) + Gi>J * B(KS)

A(KS) = A(KK) - R£. A(KK) = A(KK) + RE

3(KS) = aCKK) - FIM

B(KK) = a(KK) FIM

KK s KK + JSPANH

KS = KS JSPANH

FIM = SM * A(KS) + CM * B(KS)

RE * CM * A(KS) - SM * t3<KS)

ACKS) = A(KK) - RE

ACKK) = ACKK) + RE

3CKS) = 3CKK) - FIM

3CKK) = 3CKK) + FIM

KK « KS + JSPAiNJH

IFCKK-N)2a*304*304

334 CONTINUE

30 KK = KK - NN

JJ = JJ • K

IFC JJ-NQ)305*80>80

305 CJNTINUE

40 CD s * CN + CD

CN = CD + CN

SM=CN*EXPS

SD = R CM + SD

CM = SD + CM

SN=-CM*£XPS

GJ TO 20

30 K s K/2

RAO=.b*RAD

90 CONTINUE

170 RETURN

END

16
SUBRJLiriNE REALTRAN(A#9>MM*;M£* IN^)

DIMEiMSIJi^ ACia24)*B(ia2 4)*JC<lb)*SrClb)

CJMMJiM M*JC*ST

1 MsMM

CALL FFTC

IF(M-20^ )i>30>l 70*b3i3

500 CJNTINUe

K=JC(M+l)

N«K

iMH=N/2

NK=MH+l

CN = ISIGf^(l > IN\/)

S»^=ISIGiN)(l *iM£)

IF(CN)301 #3^3

301 CONTINUE

2 f IM=A(1 )-A<N+l ) Ad )=A(1 )+A(N+l ) 8(1 )=FIM

GJ TJ 4

3 A<N+1 )s2.*(A(l )-B(l )

)

ACl )«2.*(A<1 )+B(l )

)

B(N+1 )=0

8(1 )=0

4 IF(M)302*170*302

302 CJNTINUE

SD*CN*SN*ST (M)

F?a2 .*sr (M+l ) R=-R*R

CD=-.5*CN*R

SN=0.

DJ 5 J=:2*NK

CD = R * CN + CD

CN = CD + CN

SO = R * SN + SD

SN = SD + SN

AA = A (J) + A(K)

A8 = A(J) - A(K)

BA=3(J)+8(K)

3B = B(J) - B(K)

RE = CN * BA + SN AB

FIM = SN * BA - CN * AB

B(K) = FIM - BB

8(J) = FIM + 88

A(K) = AA - RE

A(J) = AA + RE

5 K=K-1

170 RETURN

END

17
Two FORTRAN II programs are described in this section. They

are used to modify the time domain waveform data stored in local

memory, making these data suitable for frequency domain transformation.

The first program, stored in a user's file in the time share

system under file name /XXI/, performs three functions. However,

before a detailed description of these functions is presented, it is

necessa^ry to discuss the format of the data retrieval from local

memory.

The time domain waveform data is stored locally in a 1000 point

BCD memory. Each Y-axis data point consists of three decimal digits.

The X-axis data is not needed since the Y data is stored and transmitted sequentially. This 1000 point waveform is converted to ASCII

teletype code and put on paper tape. It is then read into the time share

system from a paper tape reader and placed in file name /RAW/,

Program /XXI/ is then implemented as follows:

1. The first f\inction of /XXI/ is to shift the array indexing

of /RAW/, An idiosyncracy of the hardware is that the

first Y data point appears on paper tape three times in

succession. Thus it is necessary to delete the first two

data points in file /RAW/. This is accomplished by index

shifting such that M(1) = M(3), M(2) = M(4), etc,

2, At this point, /RAW/ consists of 1000 three-digit Y data

points representing a time domain waveform. The FFT

algorithm requires the waveform data array to be of length

N

2 , where N is an integer. The closest length is 1024

data points so 24 points must be added to the array. This

may be accomplished in any suitable manner, but in order

18
to introduce a minimum of error to the data, the following

method has been chosen. The average value of the last ten

data points is computed. This value is then used to fill the

array by placing it in positions 1001 thru 1024. In other

words, the averaged final value of the waveform is simply

repeated 24 times.

3. The third program function is a form of signal averaging.

This portion of the program is optional. If an unaveraged

waveform array is desired it is stored in file name /DON/.

If the averaging process is required, the data output appears

in file name /DONER/. Waveform averaging is accomplished

by comparing the value of each data point with the average

value of its neighbors on either side. If the value of the

tested data point exceeds certain limits, that value is changed

to equal the neighbors' averaged value. See Figure 2. 1.

Program /XXI/ receives values for NR and NL from the TTY

keyboard at run time. From Figure 2. 1 it is apparent that NR is the

total number of neighbors nearest the test point. In this example NR =

12, or six on either side of the test point. LOAV is the computed average value of the six data points to the left of X , and LIAV, the average

of the six points to the right of X . The average value of LOAV and LIAV

then is computed and stored as LOT. Therefore, LOT is the average

value of the NR nearest neighbors of X . The absolute value of X -LOT

is then compared with the value chosen for NL. If Ix^-LOT | exceeds

NL, then X^ is reset equal to LOT. If |X^-LOT| is less than or equal

to NL, then X is unchanged. NR must be a positive even integer and

NL, a positive integer.

Figure 2-2 shows a typical time domain waveform along with the

results of averaging for various values of NR and NL using Program

19
/XXI/ The program listing is as follows:

NOTE: Because a detailed explanation of both the input data and the output data programs is contained in the text, comment statements

are not included in the program listings.

/XXI/

r DIM£i>4SIJ>J M(1024)

QIMEi^JSIJN L(1324)

JPEN (3* INPJT*/i^AW/)

JPEN (4*0jrPJr*/0JNER/)

READ 3* 1 3* (M(K) *K=l * 1 3 33

)

TYPE 25

25 FJRMAT (10HNL AND NR ?

)

ACCEPT lb*NL*NR

15 FORMAT (212)

0 J 5 1 = 1 *1000

:> i-U I )=i^( 1+2 )

D J 7 1 = 991 »\d'dd

7 SJM=SJM+M(I)

A\/G=SJM/10

OJ 8 1=1^31*1324

8 M(I)=Ay/G

CLJSE (3)

I JPEN (5* Jjr?>JT*/D JiM/)

DJ 17 J=l*1324

17 WRITE 5*23*(M(J))

13 FJRMAT (1333(13/))

23 FJRMAT (1324(13/))

0 J 35 K = l * 1324

35 L(K)=M(K)

DJ 4 3 J=l * ( 1 324-NR

)

LJ=3

LI =3

OJ 60 K = J* ( J+(.MR/2 )-l ) 60 LJ=LJ+L(K)

OJ 73 K=( J+(NR/2 )+l ) * J+NR

73 LI=LI+L(K)

LJA\/=LJ/ (NR/2 ) LIA\/ = LI/(i>JR/2 ) LJT = (LJA\/+LIA\/)/2

IF ( IA8S (L( J+(NR/2 )

) -L JT ) -NL) 43 * 43 * 2

21 L( J+(iMR/2 ) ) = LJT

43 CJiMTINJE

OJ 83 J=W1324

83 WRITE 4*23* (L( J)

)

ST JP

END

20
The second data input program is called /MODI/. The input

data to this program is listed in file name /DONER/ and is the output

from program /XXI/. Since the FFT algorithm will produce erroneous

results for non-periodic waveform inputs, this program insures that the

input waveform will appear periodic. This is accomplished by first

deleting all the odd-indexed data points in the input waveform. Then

the resulting 512 point waveform is inverted, vertically shifted, and

placed in indices 513 through 1024. See Figure 2. 3 for a sample waveform both before and after being made periodic.

The output of this program is then in proper form to be transformed using the FFT subroutines. The /MODI/ listing is as follows:

/MJOl/

DIMENSION K<1024)

DIMENSION L(1024)

DIMENSION M<ia24)

DIMENSION N<1024)

OPEN <3* INPUT */OONE!^/)

READ 3*10* CMC J)# Jsl > I 324)

10 FORMAT CI3)

DO 5 J=l*5l2

5 NCJ)=:M<J*2)

ISJM=0

DO 20 J=l * 10

20 IS'JM=IS'JM-»-N( J)

IA^G=ISUM/10

• ISHIFT = N(5l I )-IA\/G

DO 30 J=l *512

30 L(J)=N(J)

DO 40 J=l *512

40 K< J)=2*N(51 1 )-L( J)-lSHIFr

DO 50 J=l *512

50 N( J+512)=KC J)

OPZi^ (4*OJTPUT*/O0N£/

)

00 70 J=l * 1 024

70 WRITE 4*60*N<J)

60 FORMAT (13/)

STOP

END

21
2. 3 Output Data Programs

The first data output prograna is called /TRI/. This program

accepts the output of /MODI/, a 1024 point averaged, periodic timedomain waveform stored in file name /DONE/ and performs the FFT of

this waveform. The data output from this program is available in the

following files:

/DONEA/ = A(r)

/DONEB/ = B(r)

/DONEX/ =

I [A(r)]^ +[B(r)]^ .^^^

where S(r) = A(r) + jB(r)

and S(r) are the 512 xinscaled complex frequency coefficients associated

with the 1024 point input waveform. /DONEX/ is really a scaled discreet magnitude spectrxxm of the input waveform. Scaling is accomplished

by multiplying each A(r) and B(r) (with the exception of A(0)) by 2/N where

N is the length of the transform. In this case N = 1024. A(0), or the

D. C. term, must be multiplied by 1/N for proper scaling.

The second data output program is called /LSP/ and is used to

scale the output magnitude spectrum /DONEX/ so it may be returned to

local memory and displayed on a CRT. The scaling equation used to yield

an integer spectrum in the range 0 ^ M(r) < 1000 is:

M(r) = B'(r) = 200 1ogjQ[l00 |S(r)|]

where |S(r)| are the scaled magnitudes of the complex frequency coefficients, B (r) are the rescaled decibel values of |S(r)|, andM(r) are

the nearest integer values of B'(r). See Figure 2. 4 for some sample

displays of magnitude spectra in dB versus frequency.

22
The third data output program, called /UNTRI/, is used to perform the inverse FFT on frequency domain data. The inputs to this

program are the unsealed complex frequency coefficients from /TRI/.

They are contained in file names /DONEA/ and /DONEB/. The outputs

of this program, contained in /Yl/ and /Y2/, are the first and last half,

respectively, of the time domain waveform associated with /DONEA/ and

/DONEB/. The /UNTRI/ program was used only to check the validity of

the FFT program package. The program listings for /TRI/, /UNTRI/

and /LSP/ are as follows:

OIMiNSIJN MC1324)>A(513),"3(513)>C(1.3:'4).0(513)

R£AD 3>13* CM<K)^K=1 > 1024)

13 FORMAT (13)

OJ 43 J=l > I 324

40 C(J)=M(J)

DJ 53 J=l ,512

53 A<J)=C<J)

DO 60 J=l *512

63 8CJ)=C< J+512)

30 FORMAT <F1 1 .2)

CALL REORDER (A>a#9)

CALL CFFTRC <A,3*9# .5*1 .* 1

)

CALL REALTRAN CA*8^9*1,1)

CLOSE (3)

OPEN <4>0UrPUT> /DONEA/)

WRITE 4*33* (A(J)> J=l >513)

OPEN <5*0UrPUT>/00NEa/)

WRITE 5>33* (3(J), J=l *513)

DO 100 J=l,513

133 DC J) = CSQRT< CA<J)**2) + (8<J)**P) ))/5iP,

CLOSE (4)

CLOSE <5)

OPEN (6,0UTPUT^/00NEX/) • WRITE 6*120> (DCJ)* J=l >513)

123 FJRMAT (F7.2)

STOP

END

23
OlMENSIJiNJ A<5l3):.3(bl3)

J-^EN C3 > INPJr^/D JiMEA/)

J-^EN (4> lNPJT*/0 JiNJES/)

^EAO 4*1 3> (8CK)*K=1 *513)

13 FJ-Ri-JAr CF11.3)

CALL REALr.^AiM (A*B* 9> -1 , -I ) CALL CFFrS CA>8# 9^ 1 ./ 1 354. > -I

)

CALL REJ^OEr^ (A*3*9)

CLOSE (3)

CLJSE (4)

J-'EiM (5> JUrPUT^/yi /)

JPEN (6* JjrPUr*/Y3/)

W^ITE 5>1 0> (A<K)>K=1 *D13

)

w^irE 6>l'a* (a(K)jK=l^:5l3)

STJP

END

NOTE: In the listing which follows B(K) represents the term

B'(r) as defined on page 22.

/LS-^/

0 I MEiNJS 1 A ( 5 1 3 ) * 3 ( 1 4 ) * H ( 1 J a 4 ) OPEN <3* IiMPiJT>/D JiMEX/

)

READ 3*l'^> CA<K)*K=1 *:>13)

A(l)=3.

OJ 15 K=l >256

15 AC2*K)=2k)0.*<ALJG ( 1 3<3 .*A (2*K) ) )/ CALJGC 1 J . J) ) D J 25 K=l ^512

3 (2*K-l )=0.

25 8C2*K)=ACK)

DO 35 K=l >1024

35 M(K)=8CK)

OPEN <4*0UTPUT>/SPEC/)

WRITE 4#22l* <M(K)*K=sl ^ 1 324)

OPE,i^ (5*0UTPLfT*/SHJRr/)

WRITE 5*20* CM<4*K)*K=1 *256)

10 FORMAT <F7.a)

20 FORMAT <I3)

STOP

END '

24
3. EXPERIMENTAL CONCLUSIONS

3. 1 Problems Encountered Using the FFT

The primary goal of the work upon which this report is based is to

develop both hardware and software for a particular system. At this

point in time, the main desired function of the system is to produce the

gain versus frequency characteristics of certain electrical networks

using time domain techniques. As a result, only specific FFT application problems, pertinent to this particular system, will be discussed.

Problems of a more theoretical nature are discussed in detail in the

literature, e.g., [3] or [4].

Basically, the system works in the following manner. A small

transition- time voltage step-waveform is generated, observed with a sampling oscilloscope, and the oscilloscope analog output recorded.

The FFT of the recorded waveform is then computed, and a magnitude

versus frequency graph is recorded. The voltage step-waveform is then

applied to the input terminals of the network under test (e. g. , a coaxial

attenuator); the network output waveform is observed with the sampling

oscilloscope, recorded and transformed. The logarithmic difference

of these two transforms, then, is computed to yield the gain versus

frequency characteristics of the network.

Actually, four major problems or limitations were encountered

which had a detrimental effect on network measurements. Each problem

shall be discussed individually.

The first problem was due to the choice of a voltage step as an

input waveform. The FFT algorithm assximes the input waveform to be

periodic. Thus an unwanted discontinuity will appear between the

beginning and end points of the step waveform, implying that the voltage

step was instantaneously turned off. The FFT will transform the instantaneous transition as well as the valid voltage step waveform data,

yielding a grossly invalid frequency spectrum.

25
The objective, then, was to make the step waveform appear

periodic while introducing neither discontinuities nor other non- causal

voltage changes. This was accomplished through the use of FORTRAN

- II program /MOD 1/. This program modifies the input step waveform

in such a manner that it appears to the FFT as if the voltage step were

turned off in exactly the same manner in which it was turned on. (See

Figure 2.3). Consequently, the modified waveform appears periodic

with no discontinuities and no new non-causal data added to it. The second problem was also solved using program /MOD 1/.

Figure 3. 1 consists of two photographs of the magnitude spectra in dB

of two ideal rectangular pulses of arbitrary width. These spectra are

both difficult to interpret since they are made up of two curves superimposed upon one another. In fact, one curve on each photograph

represents the magnitude spectr\im of the even harmonics while the

other curve represents the magnitude spectrum of the odd harmonics.

It was found that the manner in which program /MOD 1/ produces

a periodic waveform from a voltage step waveform also ensures that

the even harmonic magnitude spectra will be forced to zero and the

odd harmonic magnitude spectra will decrease monotonically as in

Figure 3. 2. In terms of Fourier series analysis, program /MOD 1/

generates a single period of a periodic waveform in which

f(t) = - f (t + ^ ) (3. 1)

where T is the period of the program- generated waveform. From

Fourier series analysis, any periodic waveform that fulfills equation

(3. 1) will contain only odd harmonics.

The third problem is really a limitation inherent in the use of a step waveform as a driving function. The step waveform must be made

periodic in order to satisfy the requirements of the FFT. Thus, the

26
modified waveform will always resemble a rectangular pulse. But the

Fourier integral transform of a rectangular pulse is of the form

six x/x where x is a function of both frequency and the width of the

time window. Consequently, there is an inherent 1/x decay of

signal input level as a function of frequency. This signal level decay

may be overcome somewhat by careful selection of the time window

used, and possibly by pulsed RF techniques.

The fourth major problem is that of noise. This problem is not

caused by use of the FFT, but rather, it is caused by the data acquisition

system. It is worth mentioning because several software problems

result in processing the resultant noisy signals. The magnitude spectra

shown in Figure 2. 4 are noisy spectra typically encountered. As may

be seen, only the very beginning of each spectrum appears to show any

continuity. The rest is much too noisy to be interpreted.

Aside from attempts to diminish the noise sources within the hardware, two software techniques were used to bring the noise level down.

One was signal averaging as in FORTRAN II program /XXI/. The other

was taking transforms over several decreasingly smaller time windows

to obtain usable frequency domain data over an increasingly wider

frequency band.

3. 2 Time Domain Measurement Example

In this section the measurement of the gain versus frequency

characteristics of a coaxial 10 MHz low-pass filter is presented. Figure

3. 3 is a block diagram of the waveform data acquisition system used to

perform this measurement. For purposes of this report, this data

system simply acquires and stores in local memory a 1000 point representation of a time domain waveform. The input step waveform and the

low-pass filter output waveform are shown in Figure 3. 4. In the example,

the input waveform is unsmoothed while the output waveform has been

smoothed using FORTRAN II program /XXI/.

27
The FFT software programs were then applied to both waveforms

and the logarithmic magnitude versus frequency spectrum of each waveform was recorded in local memory. An X-Y recording of these two

spectra is shown in Figure 3. 5. This figure clearly shows the 1/x

decay characteristics of the transformed step waveform. Also, it

exhibits a relatively noiseless gain window of about 60 dB below the

fvmdamental frequency amplitude.

Figure 3. 6 is an X-Y recording of the algebraic difference of the

input and output spectra shown in Figure 3. 5. This figure shows a fairly

flat response of the low -pass filter out to approximately 9 MHz. The

gain then drops at a rate of about 20 dB/MHz into the noise.

The results of this measurement appear to be promising. The

purpose of this measurement was merely a qualitative check of the basic

system hardware and software, but the resulting curve of Figure 3. 6

shows a noise level and a frequency resolution that indicates the feasibility of the time- domain calibration system.

3. 3 Future Developmental Considerations

With the goal in mind of developing a time domain network calibration system competitive with current capabilities of frequency domain

calibration systems, much more work is necessary. Outlined below is a schedule of tasks deemed necessary to meet this goal.

1. The first task, that of implementing a dedicated minicomputer

into the system, is already underway. The minicomputer will

allow much higher speeds of data handling and processing as

well as increased system automation capabilities.

2. After implementation of the minicomputer, a more detailed

and exhaustive set of experiments aimed at error analysis

may be performed. At present, use of the time share computing system is much too time consuming to permit an indepth study of systematic or random errors.

28
3, With sources of errors pinpointed, software may then be

developed as far as possible both to determine the magnitude

of and to reduce the effects of those errors.

4. Concurrent with software development, the fourth task would

be hardware development. This would involve development

of fast- rise noiseless step generators, improved sampling

systems, increased data acquisition resolution, and increased

system automation.

The general concept of using the FFT to compute the gain versus

frequency characteristics of electrical networks appears quite feasible.

Most probably, the real limitations on measurement precision and

accuracy lie in system hardware, and not in FFT programming or other

software.

This report has been concerned with three main topics: a definition

of the fast Fourier transform; a description of software programs used to

compute the FFT on the in-house time share computing system; and a brief description of problems encountered in the data reduction process

used.

The next report on this continuing work will be addressed to the

actual use of the system for quantitative measurements of amplitude and

phase characteristics for microwave attenuators, filters, and mismatches.

These measurements will also be compared to those made at other

laboratories using time domain techniques (e, g, , see Ref. [6] and [7] ).

29
REFERENCES

1. Cooley, J. W. and Tvikey, J. W. , "An Algorithm for the Machine

Calculation of Complex Fourier Series", Mathematics of Computers, Vol. 19., pp. 297-301, April, 1965.

2. Cochran, W. T. , et al. , "What is the Fast Fourier Transform? ",

IEEE Trans, on Audio and Electroacoustics, Vol. AU- 15, No. 2,

June, 1967.

3. Cooley, J. W. , et al. , "Application of the Fast Fourier Transform to Computation of Fourier Integrals, Fourier Series, and

Convolution Integrals", IEEE Trans, on Audio and Electroacoustics, Vol. AU-15. No. 2, pp. 79-84, June 1967.

4. Bergland, G. D. , "A Guided Tour of the Fast Fourier Transform", IEEE Spectrum, pp. 41-51, July, 1969.

5. FFT subroutine package written by David L. Lewis, National

Oceanic and Atmospheric Administration, Space Environment

Laboratory, Boulder, Colorado 80302.

6. Nicholson, M. A. , et al. , "Applications of Time-Domain Metrology

to the Automation of Broad-Band Microwave Measurements", IEEE

Trans, on Microwave Theory and Techniques, Vol. MTT-20, No. 1,

January 197 2.

7. Cronson, H. M. , et al. , "Time Domain Measurement Study",

Sperry Rand Research Center, Sudbury, Massachusetts, August,

1972.

30

inCO

o"

II

•-Iu>uo

(U>

M-to

cn

•+-> I—

:it»

<u

COO>auo>o

CO

.1-1 O1—

au

•i-i

bio

32

34
cr

oa

uo

<u

too

?J

<u

?scr

<i)

<+^

CO

>

<i;

lis

-I-l fl

two

o

01

CO

• r-l

t—

I

d

o

•CO

*0

-1-1 F=4

M

Om

>

.1-1

ao

35

37
W

H

OMPU J

oI—I >-

I— UJ q;

llj ct: o

r— ^ O 2:

C3 O LjJ

LlJ \ U-l O

O

Q LU

CD Q

1—1 o

Q CQ

0 0

00 0 1—

1 UJ cc:

LjJ 0 ai 1— ^ 1—

cc: 1—f

1—

, Q (— Q- ZD Q <c 0

<: CO 00 1—

X 0

_i Q q:

CQ UJ UJ <: >- o •—' ct o

q; _i >—

ct UJ q: > Q I—

C3

I—

_l

CO

UJ X UJ

00

CO

PE CD 0

0 1—

LT) _J

Q_ _l

2: _l <c 1—

00

OSC

Q

O

CD

CD

UJ >- UJ Q _J Q 0 a. 0 1—1 Q_ 1—

Q ID Q

00

_l _J UJ 00 UJ

2:Z 1—

ZD CQ ZD 1— 1—

0 on

u_ 0

UJ h- 3 q; > < 0 0

2: 3

3 UJ 1—

UJ

• UJ

U. CD

UJ

CD

cc: _J

(—

UJ <c UJ

_i CD CJ _i

Q- CD t—

Qs: 1—

1— 2:

<: Ql <

00 1— UJ

1

f >

CD

1—

Q

_J <C

Q_ UJ

2: n: <c

00

38
AtXI 00^

39

41
APPENDIX A

To yield some insight into the mechanisms of the FFT, an eightpoint DFT will be solved below using the FFT algorithm.

Given a time domain waveform represented by eight equally timespaced data points, as illustrated in Figure A-1, we wish to compute

the DFT for this waveform. In particular, we wish to calculate the

eight A 's associated with the eight assumed X, 's. The DFT for the

r k

set of data points is.

:=o

2 TT jrk

8

r = 0, 1, 2, . . . , 7 (A-1)

For the first reduction of transform length we let

B = > Y, e

r k

8

r = 0, 1, 2, 3 {A-2)

and

k=o

4 ff j rk

8

r = 0, 1, 2, 3 (A-3)

where

and

k 2k

\ - ^2k+l

(A-4)

(A-5)

Using the simplifying notation that

-2jrl

N W = e we have

(A-6)

. = > X, W r k

rk

r = 0, 1, 2, . . . , 7 (A-7)

42
^ 2rk

and B =^ Y^W ^ r = 0, 1, 2, 3. (A-8)

k=o

=J2 Zj^^^""^ r = 0, 1, 2, 3. (A-9)

k=o

Thus, we have now reduced the problem to that of finding two four-point

DFT's rather than one eight-point DFT. (See Figure A- 2).

For the second reduction, B and C may each be reduced to

r r two two-point transforms as follows:

Let

^ D R W^^^ r = 0, 1 (A-10)

k=o

and E = VsW^^^ r = 0, 1 (A- 11)

r k

k=o

where R^^ = Y^^ = X^^ (A-U)

and \ = ^2k+l = ^4k+2 '^"l^'

and let

4rk

F = V W r = 0, 1 (A-14)

r ^

"Z) ^^^""^ r = 0, 1 (A-15)

k=o

where T^^ = Z^^ = X^^^^ (A-16)

43
Then D is the two point transform of data points X and X^, E , r o 4 r the transform of X^ and X., F , that of X, and X^, and G , that of

2 6 r 15 r X^ and X^.

3 7

The third reduction, which would yield eight one-point transforms, is really unnecessary for the DFT of a single point function is

the sample itself. Therefore, substitution of the data points into (A- 10),

(A- 11), (A- 14), and (A- 15) yields the following set of equations.

D = X +W°X F = X, +W°X

0 o 4 o 1 5

D, = X - W°X, F, = X - W°X^ (A-18)

1 o 4 115

E = X^ +W°X, G = X^ +W°X^

o 2 6 o 3 7

As was shown in the text,

and

A = B +W^C r = 0, 1, 2, 3 (A-19)

r r r

A, ^. = B - W^C r = 0, 1, 2, 3 (A-20)

r r

Likewise, it may be shown that

B = D +W^^E r = 0, 1 (A-21)

r r r

44
and

= D - W^'^E

r r

C = F +W^^G r r r

Recalling that

W = - Wh)

r = 0, 1

r = 0, 1

r 0, 1

{A-22)

(A-23)

{A-24)

(A-25)

we may substitute (A- 18) into (A-21) and (A-22) and then into (A- 19)

and (A- 20) to yield the following matrix equation:

A

o

1 1 1 1 1 1 1 1 X

o ^1 1

. 1 ^1

1 - 1 1 - 1

^3 1 + - 1 -w' ^3

A

4

1 -1 + 1 - 1 + 1 - 1 + 1 -

X

4

^5 1 -w' - 1

^6 1 -W^ -1 1 -1

^7 1 -w^ - 1 + 3

(A-26)

45
Examining (A- 26) we see that half of the multiplying matrix terms are

vinity. In addition, the matrix is symmetrical about the main diagonal.

Consequently, about three -fourths of the multiplications are unnecessary. Figure A-3 is a signal flow graph illustrating this computation.

46
CO

<J3

un

CO

O0

<p

CO

OkoaI

•r-l w

J

L o J

X I • r-l

47
o

Fig, A- 2 Two four-point representations of the same time domain

waveform.

48

FORM NB5-114A (1.7«) U.S. DEPT. OF COMM. 1. PUBLICATION OR REPORT NO. 2. Gov't fV: c;DssiGn BIBLIOGRAPHIC DATA NBSIR 73-303 ***** 5n c E T

3. Recipient's Accession No.

4. TITLE AND SUBTITLE

Fast Fourier Transform Implementation for the Calculation of

Network Frequency Domain Transfer Functions from Time

Domain Waveforms

5, Publication Date December 1972

6. Performing Oganizatim Code

/• AuinuK(.a^ w liiiam Xj. vjans ana iN. o. iNanman 8« Performing Organization

NBSIR 73-303

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS, Boulder Labs

DEPARTMENT OF COMMERCE

Boulder. Colo, ovovd

10. Project/Task/Work Unit No. 2722390

11. Contract/Grant No.

12. Sponsoring Organization Name and Address

Department of Defense

Calibration Coordination Group

CCG 72-65

13. Type of Report & Period

Covered

Final

14. ^onsoring Agency Code

15. SUPPLEMENTARY NOTES

16. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here.) This report is concerned with the software applications of the fast Fourier

transform algorithm to the relationship between time domain waveforms and frequency

domain spectra. The first chapter is devoted to a description of the discrete Fourier

transform and the fast Fourier transform. Chapter 2 contains the text and a brief

description of all FORTRAN II programs utilized in connection with this work. All

computation was performed on the in-house time share computing system in the NBS

facilities, Boulder, Colorado. In Chapter 3, problems encountered using the fast

Fourier transform algorithm are discussed, an example of a time domain to frequency

domain calculation is presented, and future developmental considerations are

mentioned. In addition Appendix A contains a detailed example aimed at disclosing

the inner mechanisms of the fast Fourier transform algorithm.

17. KEY WORDS (Alphabetical order, separated by semicolons) Discrete Fourier transform; fast Fouricr

transform; frequency spectra, discrete; network transfer function; time domain

waveform; transfer function.

18. AVAILABILITY STATEMENT 19. SECURITY CLASS

(THIS REPORT)

21. NO. OF PAGES

[5 UNLIMITED.

UNCLASSIFIED

1 1 FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE

TO NTIS.

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

22. Price

USCOMM-DC 66244-P7I

IISCOMM - ERL