Calculating the DFT

The most obvious procedure for computing the DFT is a direct numerical evaluation of Eq. (4.3.6).

N — l

F’(j(ok) = (l/N) £ xpexp(-j<oktp) (4.4.1)

P= o or

F’(j(ok) = (l/N) £ *pexp(-jpOk) (4.4.2)

P= 0

where 8k = 2kn/N. The real and imaginary parts are

Rk = (l/N)l xpcosPek (4.4.3)

P= o N — 1

Ik = (- l/N) £ Xp sin pOk (4.4.4)

P-0

This analysis would require 2N multiplications for each frequency, and the sine and cosine functions would have to be calculated or looked up in a stored table. Of course, these trigonometric functions can be obtained by up-dating prior values using the following identities:

cos(p + l)6t = cos pdk cos вк — sin pdk sin 0k sin(p + l)6t = sin pdk cos вк + cos pdk sin 6k

A less obvious but more efficient procedure (3, 5) may be developed by using a simple recursion relation that gives a result that is easily transformed into the desired Fourier integrals. The recursion relationship involved the calculation of a new quantity, C(p), defined as follows

C(p) = x(p) + 2 cos вк C(p — 1) — C(p — 2) (4.4.5)

and with C( —2) = C( — 1) = 0. This defined relationship may be solved for x(p) and substituted into Eqs. (4.4.3) and (4.4.4) to give

Rk = (l/N) Y [C(p) — F(k)C(p — 1) + C(p — 2)] cos рвк (4.4.6)

P = o

Ik = -(l/N) Y [C(p) — F(k)C(p — 1) + C(p — 2)] sin P6k

P = 0

where F(k) = 2 cos 6k. These series may be expanded and rearranged to give

cn-з

Я* = (1 /N)< Y C(p)[cosp0* — F(k) cos(p + 1)0* + cos(p + 2)0*]

Ip = o

+ C(N — 2)[cos(N — 2)0* — F(k) cos(N — 1)0*]

+ C(N — l)cos(N — 1)0* j (4.4.8)

J

rN-3

Ik = -(l/NH Y C(p)[sinp0* — F(k) sin(p + 1)0* + sin(p + 2)0*] lp = 0

+ C(N — 2)[sin(N — 2)0* — F(k)sm(N — 1)0*]

+ C(N — 1) sin(iV — 1)0*| (4.4.9)

The advantage of casting the equations in this form is that the summation term in each expression is zero because of the following trigonometric

Подпись: identities: cos тф — 2 cos ф cos(m + 1)ф + cos(m + 2)ф = 0 (4.4.10) sin тф — 2 cos ф sin(m + 1)ф + sin(m + 2)ф = 0 (4.4.11) Thus Eqs. (4.4.8) and (4.4.9) become Rk = (l/N)C(N - 2)[cos(N - 2)0* - 2cos0*cos(lV - 1)0*] + (1 /N)C(N - l)cos(N - 1)0* (4.4.12) /* = -(1 /N)C(N - 2)[sin(N - 2)0* - 2 cos 0*sin(lV - 1)0*] - (l/N)C(N - l)sin(N - 1)0* (4.4.13)

These may be further simplified using Eqs. (4.4.10) and (4.4.11) and the following relations:

sin(N — 1)0* = sin Л/0* cos 0* — cos Л/0* sin 0* cos(N — 1)0* = cos N9k cos 0* + sin Л/0* sin 0* cos Л/0* = cos 2kn = 1 sin Л/0* = sin 2kn = 0 The final results are

Rk = (l/N)[-C(N — 2) + C(N — 1)cos 0*] (4.4.14)

/* = (l/N)[C(N — l)sin0*] (4.4.15)

We observe that the algorithm only requires the evaluation of C(p) for each data point for each frequency using Eq. (4.4.5). This involves a single

multiplication and two additions. After all the data points have been treated, the sine and cosine integrals are obtained from the last and next-to-last terms in the C(p) calculation using Eqs. (4.4.14) and (4.4.15). The advantage of this method for on-line analysis using a dedicated digital computer is obvious. The procedure is :

1. Set up a table of the F(k) function. This will require storage of one number per frequency to be analyzed.

2. Begin sampling the signals. The samping rate must be slow enough to allow the evaluation of C(p) between sampling times, using

C(p) = x(p) + F(k)C(p — 1) — C(p — 2)

3. When the data sampling is complete, convert the values of C(N — 2) and C(N — 1) into the desired Fourier integrals, using Eqs. (4.4.14) and (4.4.15). This has no bearing on the allowable sampling rate in an analyzer, since it is done after the sampling is completed.

A BASIC language computer program for the algorithm of this section is given in Appendix B. The simplicity of the programming required is illus­trated by this program.