This chapter gives a short description of the frequency response functions that describe nuclear reactor dynamics. A power reactor has inherent feed­back due to temperature (and possibly pressure) changes that accompany power changes. This closed-loop system may be modeled using a forward loop to describe the power response to a reactivity change and a feedback loop to describe the power-to-reactivity feedback processes.

The forward-loop transfer function gives the response to a reactivity input that would occur in the absence of any feedback. This transfer function is called the zero-power transfer function.

The most general form of model that has been used to calculate the zero — power frequency response includes the energy and spatial dependence of the neutron flux as well as the frequency dependence. In this case, the multi­group, time-dependent neutron diffusion equations are used. The space dependence of the neutron flux can be important at higher frequencies, but this usually does not influence tests on power reactors because the main emphasis in these tests is to identify feedback effects, and these influence the response at lower frequencies. In special measurements, the higher-frequency response may be measured to provide reactor physics data (such as neutron generation time). Space dependence is also important at very low fre­quencies (less than 10“5 Hz) for xenon spatial transients.

There is also some interest in the dynamics of subcritical reactors. The degree of subcriticality can be determined in tests in which the reactivity or a neutron source is modulated.

lie

The pseudo-random binary sequence was the first periodic, discrete-level signal used to measure the frequency response at a number of harmonic fre­t See the literature (1-55).

quencies simultaneously. In addition to numerous tests on nuclear reactors, this signal has been used in tests for distillation columns (5,19, 39), oil refineries (44,47), evaporators (19), heat exchangers (7), steam generators (21), diesel engines (43), machine tools (29), electric furnaces (36, 37), paper mills (45), and steel mills (45).

Two different types of PRBS signals, called the quadratic residue sequences and the maximal-length or m sequences, are available. Quadratic residue sequences with Z bits may be generated for Z = 4/c — 1 for integer values of к if Z is a prime number; m sequences with Z bits may be generated for Z = 2n — 1 for integer values of n. Since m sequences may be generated easily using a shift-register algorithm, only these will be considered here. A description of the method for generating m sequences appears later in this section.

The two main characteristics of interest for a PRBS are its power spectrum and its autocorrelation function. These are as follows.

(1) Power Spectrum

Pk

Pk

where Pk is the power in the fcth harmonic, A the signal amplitude (the input is A or — A), Z the number of bits, and к the harmonic number. For example, two periods of a 7-bit PRBS appear in Fig. 3.1. Power spectra for several

PRBS signals appear in Fig. 3.2. In order to make the cases comparable, all are based on the same period and the same input amplitude.

where Cі i(t) is the autocorrelation function, т the lag time, and T the period. The general form of the autocorrelation appears in Fig. 3.3.

From these results we can deduce the following:

1. The autocorrelation function has a spike at 0, T, 2T,… and a negative bias of 1/Z between spikes. As the number of bits increases, the spike becomes sharper and the bias becomes smaller. Thus the autocorrelation function looks more like a series of delta functions as Z increases.

2. The power spectrum becomes flatter as Z increases. This may be useful in tests in which it is desired to measure the frequency response at a number of harmonic frequencies.

3. The absolute magnitude of the harmonics decreases as Z increases. This is expected, since the same total power is distributed rather evenly over more harmonics for the longer sequences. This suggests that a compromise may be required in selecting a sequence for a particular application. Long sequences give the desired flat spectra, but the signal-to-noise ratio at each harmonic may be too low because of inadequate signal strength.

We observe that the power spectrum of the PRBS is smaller than the maxi­mum possible power spectrum of a binary pulse chain (Eq. 2.12.7) by a constant factor, (Z + 1)/Z2. Thus the bandwidth relation derived in Section 2.12 applies for the PRBS. That is,

fc1/2 = 0.44Z (3.1.4)

where fc1/2 is the harmonic at which the signal power is down to half of its maximum value. A rule of thumb in estimating the useful frequency range of a signal is

where At is the bit duration (in seconds). This only says that if there is enough energy in the first harmonic, then there is probably a good chance of obtain­ing results out to harmonic number 0.44Z because it has half as much energy as the first harmonic.

4. The PRBS is not antisymmetric. Because of this, the PRBS does not have the tendency to discriminate against nonlinear contamination possessed by the antisymmetric signals. (See Section 2.13.)

The pseudo-random binary sequence may be generated with a digital shift register with modulo-2 adder feedback. Modulo-p addition is defined for any two digits m and n, where m = 0,1,…, p — 1 and n = 0,1,…, p — 1. The modulo-p sum of m and n is the same as in common arithmetic if (m + n) < p. If (m — I — n) is equal to p, then the modulo-p sum is given by m + n — p. Thus modulo-2 addition is defined as shown in Table 3.1.

A three-stage digital shift register with feedback from stages 2 and 3 is shown in Fig. 3.4. The operation of this digital shift register with modulo-2

—H 11213

Fig. 3.4. A three-stage shift register with feedback.

adder feedback proceeds as follows:

1. Set initial values (0 or 1) on the stages of the digital shift register. They may be selected arbitrarily except for the requirement that at least one of them must be 1.

2. Generate a new set of digits for the stages of the shift register as follows.

a. Add modulo-2 the values in specified stages of the shift register.

b. Shift all existing values in the shift register one stage to the right and insert the feedback term from step (a) into stage 1. Note that the value in the last stage is removed when the shifting is accomplished.

3. Repeat until the initial values reappear in all of the stages of the shift register.

This set of operations produces a series of digits (0 or 1) for each stage in a pattern that repeats periodically. The pattern in each stage is the same, but the patterns are offset relative to one another. The 0 and 1 indicate which value of the two-level input signal to use. We may interpret either 0 or 1 as the high-level input or the low-level input. For instance, the three-bit sequence (1 0 1) may be experimentally implemented as (high, low, high) or (low, high, low). The correlation functions and power spectra are identical in each case.

It is known that the maximum number of bits in such a periodic sequence generated by an n-stage shift register with modulo-2 adder feedback is 2й — 1. Other shorter periodic sequences are obtained for certain feedback connections, but only the maximum-length sequence has the desired pseudo­random character. The arrangement of the feedback connections needed to generate maximum-length sequences are known for a number of sequences. In some cases, several different feedback connections will work. Some of the sequences may be generated with one modulo-2 adder and others require several. Table 3.2 gives feedback connections that require only one modulo-2 adder if a single adder connection is known. For others, connections are given for the two additions from the numbers in three stages.

As an example, consider the arrangement for a three-stage digital shift register shown in Fig. 3.5. Take the initial conditions on the three stages as 10 0. Then the following sequence is obtained.

1 0 0 Initial

1 1 0 1 1 1 0 1 1 1 0 1 0 1 0 0 0 1

1 0 0 Cycle repeats

Thus a 7-bit pseudo-random binary sequence is (1 1 1 0 1 00).

Two properties of pseudo-random binary sequences are:

(1) Zero Crossings. For a Z bit sequence, there are (Z + l)/2 zero crossings. This gives the number of shifts required of the input hardware per period and may be useful for assessing possible wear on this hardware during a test.

(2) Runs. Successive occurrence of one of the states is called a run. There are (Z + l)/2 runs per period. Half of the runs are of one-bit length; one — fourth are of two-bit length; one-eighth are of three-bit length; and so on, provided that the number of bits generated is greater than 1. In addition, there is one run of length n, where Z = 2" — 1.

The surface has only been scratched in exploiting the possibilities for extracting useful information from dynamics tests. They contain all of the information in static tests (steady state is only one point in a transient) in addition to other information not available from static tests. Also, it is easier to obtain certain types of information from dynamics tests than from static tests. This is because it is possible to separate effects due to processes with different time constants in a dynamics test, but not in a static test.

The development of power-reactor testing methods that are easy and inexpensive to implement along with the development of system identifica­tion techniques (primarily from work in nonnuclear industries) should result in common use of these techniques.

The calculation of Fourier transforms on a digital computer requires the numerical evaluation of

F(ja>) = (1/T) j" x(t)e~i0,,dt (4.3.1)

Jo

or

Re{T(yco)} = (1/T) j" x(t) cos tot dt Jo

lm{F(ja>)} = -(1/T) f x(t) sin cot dt Jo

Any method for approximating the integral using data obtained at discrete sampling points is called a discrete Fourier transform (DFT). The most common method for calculating the DFT is to approximate the integral as

follows. (Other methods are discussed in Section 4.5.’

x(t)exp(-jcot) dt ^ X xp exP( At/T (4.3.4)

p = 0

where xp is the pth sample of x(t), N the number of samples, and At the time interval between samples. Since N = Т/At, we write

F'(jco) = (l/N) £ xp exp( —jo)tp

P= o

where F'(j(o) is the DFT approximation to F(jco).

Since the DFT is only an approximation to the true Fourier transform, it is necessary to evaluate the consequences of discrete sampling. After a signal is sampled, it is not possible to detect a component with frequency above 1/2АГ cycles/time or n/At radians/time (see Fig. 4.4). This frequency is called

the Nyquist frequency. If the signal contains components above the Nyquist frequency, these higher-frequency components not only go undetected, they also contaminate the lower-frequency results. Since these higher-frequency results give apparent results for “names” (frequencies) different from their own, this effect is called aliasing.

The aliasing effect will be analyzed for the DFT approximation of Eq. (4.3.5) for a periodic signal

Fk = 0/ЛО Y xp exp( — M*p) (4.3.6)

P = 0

This may be written

CT

x(£) exp( — jcokt) d(t — tp) dt 0 ,

Fig. 4.5. The periodic delta function.

where d(t — tp) is the delta function. The delta function (see Fig. 4.5) may be written as a Fourier series (without approximation)

N — 1 oo

X S(t — tp)= X Dm exp(jcoj) (4.3.8)

p = 0 m = — oo

where u>m = 2mn/At. The evaluation of Dm gives

Dm = 1/At (4.3.9)

This series may be substituted into Eq. (4.3.7) to give

oo і» T

Fk’= X (1/Л x(t)exp[—j(a)k — wm)t] dt (4.3.10)

m = — oo * 0

The integral may be interpreted if we identify the following:

wk — wm = {2kn/T) — (2mn/At) = (2n/T)(k — mN) (4.3.11) and

CT

x(t)exp[-j{k — mN)(2n/T)t]dt = Fk_mN

0

Then we obtain

00

Fk = X F)c-miv

m = — oo

or

Fk — Fk + Fk_N 4- Fk + N + Fk_2N + Fk + 2N + ••• (4.3.14)

This shows that the DFT result at harmonic к is contaminated by higher harmonics.

It is useful to convert these results into the trigonometric form of the Fourier integral in order to interpret the results. Equations (2.5.9) and

(2.5.10) give

oo

Ak = A + X (AmN + k + AmN-k) (4.3.15)

m = 1

oo

Вк = В к + X (BmN + k ~ ^mN-k)

m = 1

where

Thus the weighting of the influence of the higher harmonics is as shown in Fig. 4.6. The first alias frequency is at N — к (m = 1). The frequency corre­sponding to harmonic N is 2л/At. This is recognized as twice the Nyquist frequency. Clearly, the DFT will be accurate only if the contributions from

к N-k N N + k 2N-k2N2N + k 3N-k3N3N + k

Harmonic

Fig. 4.6. Aliasing frequencies for analysis at harmonic k, for a sampling rate N.

the alias frequencies are small. There is nothing that can be done about aliasing after the signal is digitized. Thus the tester must make sure that there is insignificant energy in alias frequencies relative to the analysis fre­quencies. If the combination of sampling rate and frequency response of the system being tested accomplish this, then no additional action is required.

If not, then either the sampling rate must be increased or the high-frequency part of the signal must be suppressed by low-pass filtering prior to sampling. The low-pass filtering must be set to pass the frequencies of interest and to attenuate alias frequencies as shown in Fig. 4.7.

In some cases, the position-measuring device on a control rod will not be accurate enough for use in a frequency response test. If the mechanical characteristics of the rod and drive are essentially the same at zero power and at full power, then it may be possible to overcome the position-measure­ment problem.

The basis of the calibration technique is that the zero-power frequency response is known with negligible uncertainty. The measurement of the zero-power frequency response using a CRDM with unknown dynamic response may be interpreted as in Fig. 7.2. The system output 50 and the

input to the CRDM controller 5R are measured. The zero-power transfer function G is obtained from theory. Then the frequency response of the rod and drive system is obtained using

F = (50/5R)/G (7.4.1)

This dynamic control-rod calibration is subsequently used for the at-power

tests. In that case the desired at-power frequency response Gp is given by

Gp = (ёОт/F (7.4.2)

Of course, the experimenter must be careful to insure that the response characteristics of the CRDM do not change with power level. For example, if the temperatures of CRDM components change, then clearances between parts could change. This would cause changes in frictional effects that would alter the CRDM response characteristics.

Fourier transformations may be used to convert functions of time into functions of frequency. This procedure is used in analyzing data from a frequency response test. We consider periodic functions first, then non­periodic functions.

A periodic function that satisfies certain conditions! may be expressed as a Fourier series. The trigonometric form of the Fourier series is

A 00

f(t) = — Г2 + X cos comt + Bm sin comt (2.5.1)

Z m= 1

where com = Ъпк/Т is the angular frequency (radians/time) and T the period

of/(0-

The virtue of this sort of representation of the function is that certain orthogonality relations exist that facilitate the evaluation of the coefficients Ak and Bk. These are

t See Churchill (3).

tThese conditions, called the Dirichlet conditions, are that f(t) must (1) be single valued, (2) be finite, (3) have a finite number of maxima and minima, and (4) have at most a finite num­ber of finite discontinuities. It is clear that Dirichlet’s conditions are satisfied by signals of prac­tical importance.

The integration shown here has limits — T/2 to T/2 but any section of the function of duration T is suitable.

Cofficient Ak is evaluated by multiplying both sides of Eq. (2.5.1) by cos[(2kn/T)t] and integrating from — T/2 to T/2:

Therefore,

2 rT/2

= t /(0

і J — T/2

Example 2.5.1. Consider the square wave shown in Fig. 2.17. The coefficients Ak are given by

Ak = 0

The coefficients Bk are given by

Bk = —[1 — cos kn — cos kn + 1] kn

. 2n 1.6n 1 . 10я

sin—t + — sin—t + — sin—-t +
T З T 5 T

Note that only sin(o;tt) terms appear in the Fourier series representation of the square wave of Fig. 2.17. Any odd function [i. e.,/(t) = —/( —t)] has only sin(wkt) terms in its Fourier series. Even functions [i. e.,/(t) = /( —t)] have only cos(o>kt) terms. ■

It is also possible to write the Fourier series in a more compact exponential form:

oo

f(t)= £ C*expO'(ott) (2.5.7)

к = — oo

Where

Q — 2 (Л к — jBk) C-k = j(Ak + jBk)

Example 2.5.2. Compute the complex Fourier coefficients for the square wave considered previously:

1 f° 1 ГТІ2

Ck = — exp (~ja)kt)dt+ — exp (-jcokt)dt

‘ J-T/2 t J0

2 _ 2[exp(jcokT/2) + exp(—jcokT/2)]

2

Use the Euler formulas to give

C =_L_

k jcokT

Since cok = 2nk/T,

Ck = —r[l — cos kn] jnk

0 if к is even

Ck = 2

—- if к is odd

(jnk

We can compare this with the value of Ck predicted from the trigonometric form by Eqs. (2.5.9) and (2.5.10).

Ck = j(Ak — jBk) = ^(0 — 0) = 0 if к is even

= i(A — jBk) = £[0 — j(4/knj] if к is odd

Thus, the value of Ck for odd values of к is

Ck = 2/jkn

which agrees with the previous results. For negative k, we get C-k = i(A + jBk) = £(0 + 0) = 0 if к is even

= i(A + }Bk) = £[0 + j(4/kn)] if к is odd

The value obtained for — к is then

C-k = -2/jkn

which also agrees with the previous results.

We have seen that there is a definite set of Fourier coefficients associated with a given periodic function and that there is likewise a definite periodic function associated with a given set of Fourier coefficients. This unique relationship between a periodic signal and its Fourier coefficients means that the signal may be completely characterized by its Fourier coefficients (harmonics). Since the Fourier coefficients are related to a specific frequency a>k, harmonic analysis is also called frequency-domain analysis.

The set of Fourier coefficient amplitudes, taken as a function of frequency, is called the amplitude spectrum of the periodic signal. We have seen pre­viously that the Fourier coefficients of a square wave have zero values except at certain frequencies. The amplitude spectrum of the square wave is shown in Fig. 2.18 for the complex and the trigonometric forms. Note that

the complex spectrum extends from — oo to +oo while the trigonometric form extends from 0 to oo.

Modified square waves and their amplitude spectra are shown in Fig. 2.19. Note that the amplitude of the harmonics becomes smaller, and the number of harmonics in a given frequency range becomes greater, as T

increases. In general, as the period increases, the harmonic amplitude diminishes and the number of frequencies increases in the same ratio:

C(co)T = constant (2.5.11)

N/T = constant (2.5.12)

where C(a>) is the Fourier coefficient at frequency to, N the number of harmonics in a given frequency range, and T the fundamental period. The fundamental period is the shortest time T for which f(t) = f(t + T). The signal is also periodic with period 2T, ЗT,…, but an analysis using these longer periods will not change the harmonic amplitude or frequency spacing of signal harmonics.

A periodic function becomes nonperiodic if the fundamental period is allowed to become infinite. We have observed that the harmonic amplitude decreases and the number of harmonics increases as the period increases. Thus the Fourier coefficients approach zero as T -> oo. However, we noted
that the product of the Fourier coefficient and the fundamental period is a constant. Therefore, for the complex form, let us consider the function CkT as T -> oo and denote it by F{jto):

ЛТ/2

F{ja>} = lim CkT = lim f(t)exp{-j(okt)dt

T -* oo T-*ooJ—j/2

(2.5.13)

This is the definition of the Fourier transform F(jw). The inverse transform is obtained by a similar limiting process:

/(0 = lim £ (CkT) ехр(М0(1/Л

T~* 00 к = — oo

We observe that the increment in со is given by

сok+1 — cok = 2 n/T = Acu

We can substitute this into the equation for/(t) and use the fact that limj.,,^ is the same as 1ітДш_0 to obtain

J oo 1 Л00

/(t) = — lim £ (CkT) exp(M0 Act) = — F(jco)eia>t doo (2.5.14)

2,71 A(0~* <X) QQ 271 щ) —————— QQ

We saw earlier that the spectrum of a Fourier series has nonzero values only at discrete frequencies. The spacing of harmonics decreases to zero as the fundamental period goes to infinity. Thus the Fourier transform of a nonperiodic signal has a continuous spectrum. We also saw that the product CkT is a constant. Since CkT becomes F(jto) in the limiting process, we see that the amplitude of the Fourier transform is identical with the envelope of the curve through the amplitude points of CkT for the Fourier series.

The space-independent neutron kinetics equations are dn/dt = [{p — P)/A]n + £ ЛС;

І = 1

dCJdt = — ЛІСІ

where n is the neutron density (proportional to power level), p the reactivity, P the total delayed neutron fraction, A the neutron generation time, 2, the decay constant for ith delayed neutron precursor group, C, the concentration of the ith delayed neutron precursor, and the delayed neutron fraction for delayed neutron group i.

The transfer function obtained from these equations is

where n0 is the steady-state neutron density. The delayed-neutron parameters needed for theoretical calculations of the frequency response are shown in Table 5.1. Zero-power frequency response results for a reactor with 235U fuel are shown in Fig. 5.1.

An approximate transfer function based on a single effective group of delayed neutrons is often used. The transfer function is

Sn _ n0(s + X) dp As(s + X + p/A)

where X is the average decay constant (sO. l sec-*).

The n sequence is obtained by a simple modification of the PRBS. This consists of changing the sign of every other bit in the PRBS. The modifica­tion of a PRBS with Z bits (Z always odd) gives an n sequence with 2Z bits. For example, a seven-bit PRBS and the corresponding 14-bit n sequence are shown in Fig. 3.6. The power spectrum for a Z-bit n sequence is given by

7-bit PRBS

The power spectrum of the 14-bit n sequence is shown in Fig. 3.7. The n sequence may be generated by a PRBS generator with additional components added to accomplish the bit inversion.

The autocorrelation function of a general n sequence is shown in Fig. 3.8. From these results we can deduce the following:

1. The autocorrelation function has a large positive spike at 0, T, 2T,… and a large negative spike at T/2, ЗТ/2, 5T/2,… , with a string of smaller spikes between the large spikes. As the number of bits increases, the large

spikes become sharper and the small spikes become smaller. Thus the auto­correlation function looks more like a series of delta functions with alternating signs every half period as the number of bits increases.

2. The power spectrum becomes flatter as the number of bits increases, as with the PRBS.

3. The absolute magnitude of the harmonics decreases as the number of bits increases, as with the PRBS. Since the shape of the power spectrum is the same as for the PRBS, the same bandwidth relations apply (Eqs. 3.1.4-3.1.7).

4. The n sequence is antisymmetric (and thereby discriminates against nonlinear effects).

The most obvious way to perturb a reactor for a dynamics test is to change the reactivity by moving a control rod. In this chapter, the most common control-rod drive mechanisms (CRDMs) will be considered to assess their suitability for dynamics testing. The key characteristics are the speed of motion, the accuracy of the position-measuring system, and the worth of the control rod. There are five different main CRDM types: (1) magnetic jack, (2) roller nut, (3) rack and pinion, (4) locking piston, and (5) winch.

Other input signals also may be used in power reactors. Steam-flow perturbations introduced with the steam valve cause significant responses in a number of plant process variables in all types of reactors. Feedwater flow to the steam generator (or to the reactor vessel for a boiling water reactor) may also be modulated. Modern boiling-water reactors use forced recirculation, and the recirculation flow may be perturbed to excite the system. Also, special tests involving minor loops may be useful. An example is a test involving modulation of the pressurizer pressure set point in a pressurized-water reactor.

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

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.