Measurement of time-varying sinusoidal signals by a digital computer

COMPUTER PROGRAMS IN BIOMEDICINE 3 (1973) 105-111. NORTH-HOLLAND PUBLISHING COMPANY

MEASUREMENT OF T I M E - V A R Y I N G S I N U S O I D A L S I G N A L S BY A DIGITAL COMPUTER W.R.J. FUNNELL and Chas. A. LASZLO BioMedical Engineering Unit, McGill University, Montreal, Quebec, Canada

Various methods of measuring the frequency, amplitude and phase of sinusoidal signals are compared. A method is suggested which is specifically designed to exploit a general-purpose digital laboratory computer. The frequencymeasurement procedure permits accurate determinations over a wide range of frequencies. The principle of the method for measuring amplitude and phase (and also DC shift) is that one can calculate these three parameters given the values of three samples of the signal taken at known times. In this case the three samples are equally spaced in time, about one third of a cycle apart. By taking three samples from each of two signals, one can calculate the relative phases and amplitudes of the signals. The equations used are derived, and some details of an implementation for audio-frequency measurements on a PDP-12 are given. Amplitude

Phase

DC shift

Sinusoidal

1. Introduction In our experimental studies [1] we wished to use our computer (a DEC PDP-12) to measure continuously the time-varying amplitude and phase of sinusoidal signals over a wide frequency range (two decades). One approach to this problem would have been to buy or build specialized electronic hardware, such as a frequency meter and phase-sensitive demodulators, and to have the computer sample the relatively slowly varying meter and demodulator outputs. Such equipment is expensive, however, and also tends to be inaccurate if the signals being observed contain much low-frequency jitter. The latter consideration is especially important in physiological research, where the signals often jitter appreciably. For these reasons, we decided to exploit more fully the computer hardware already available by having it operate directly on the sinusoidal signals. Thus we do not have to pay for extra specialized hardware, and we can employ algorithms which minimize the effects of jitter when necessary. Specifically, using the clock and two of the analogueto-digital channels of the PDP-12, we measure the frequency of a reference signal and then measure the

Audio

Frequency

PDP-12

amplitude and phase of a second signal relative to the reference. In sect. 2 below we compare various ways of carrying out these tasks, and outline the methods we have used. In sect. 3 we describe the programmes used in more detail; while the particular programmes described here were designed for audio-frequency measurements (100 to 10,000 Hz), the principles can be used for sinusoidal signals of any frequency as long as the specifications of the computer (especially the clock) are adequate.

2. Theory 2.1. Frequency If the frequency of a sinusoidal signal is relatively low it can conveniently be measured by measuring the duration of a small number of cycles. However, as the frequency increases the duration becomes too small to be accurately measured, unless the number of cycles to be observed is increased as a function of the frequency. Similarly, if the frequency is relatively high it can conveniently be measured by counting the number of cycles that occur in some specified time

106

W.R.J. FUNNELL AND C.A. LASZLO

interval, but as the frequency decreases the number of cycles counted becomes too small unless the time interval is lengthened as a function of the frequency. If the frequency range of interest is wide (a couple of decades, for example), and if one does not know ahead of time whether the frequency will be high or low, it can be quite troublesome to apply either of these two methods. One way of combining the two approaches, so as to obtain a good measurement strategy independent of frequency, is to measure the time interval (1) required for an integral number of cycles, with the constraint that I must be greater than some minimum, Imin" (Imin is chosen to be much longer than the interval between clock pulses, to ensure accuracy.) Thus, the computer measures the duration of one cycle: if the duration is less than Imin (as it will be for high frequencies) the computer goes on to a second cycle. It continues measuring subsequent cycles until the total interval measured is greater than/min" The length of one cycle, and hence the frequency, is then obtained by dividing by the number of cycles that has been required. The mechanics of this method as we have implemented it on the PDP-12 are discussed below in sect. 3.2. Note that we have not worried about the effects of jitter on the accuracy of the frequency measurement. This is because in our application the frequency can be determined using a clean reference signal rather than a physiological signal.

2.2. Amplitude and phase To measure the amplitude of a sinusoidal signal one could simply sample as rapidly as possible for some time and then pick out the largest sampled voltage as representing the peak of the sine wave. Unfortunately this method is sensitive to high-frequency noise: a positive noise spike on top of the sine wave will make the peak look higher, but a negative spike will just be ignored. Thus, high-frequency noise will always tend to increase the amplitude estimate and cannot be averaged out. The method is also limited by the sampling rate: the shortest sampling interval of the PDP-12 corresponds to a phase change of more than 65 ° at 10 kHz (the highest frequency of concern in our particular application), so that one is not guaranteed of having any sample fall near the peak of

the sine wave. Another way of measuring amplitude, if the frequency of the signal is known, is to detect a zerocrossing, then sample after a delay equal to one quarter of the period. This method is relatively unaffected by the timing limitations of the clock since an error of, say, 10° corresponds to an error of less than 0.2 dB, due to the flatness of the peak of a sine wave. However, the method requires the clock's Schmitt trigger to be set somewhere near the centre of the sine wave, and would be adversely affected by low-frequency jitter in the signal. The most direct way of measuring the relative phase of two sinusoidal signals of known frequency would be to measure the time between a zero-crossing of one signal and a zero-crossing of the other. Just as for the amplitude measurement, however, this requires an accurate adjustment of the thresholds of the Schmitt triggers used to detect the zero-crossings, and is sensitive to any jitter in the signals. The method is also limited by the counting rate of the clock: for the standard PDP-12 clock the shortest interpulse interval corresponds to a phase change of more than 9 ° at 10 kHz. The disadvantages of these various methods may be summarized as (1) the necessity of carefully adjusting a Schmitt-trigger threshold, (2) sensitivity to signal jitter, (3) the biased detection of high-frequency noise, and (4) a low upper-frequency limit due to the rate of the clock or to the analogue-to-digital conversion time. A method which avoids these difficulties is to take a number of samples of the signal and to compute the signal parameters from the sample values. To obtain the three parameters, amplitude, phase and DC shift, requires three samples. For convenience these samples can be evenly spaced in time. If the signal were timeinvariant, and if the sampling procedure did not involve quantization errors, the spacing between the three samples could be anything except a multiple of 180 ° . However, if the signal may be time-varying the samples should be kept fairly close together to permit the measurement of rapid signal changes. On the other hand, if the samples are too closely spaced their values will be almost equal (especially if they happen to fall at the peak of the sine wave), magnifying the effects of the quantization errors. As an intuitively reasonable compromise, we have chosen to take the

MEASUREMENT OF TIME-VARYING SINUSOIDAL SIGNALS three samples about one third of a period apart. To derive the equations used to calculate the signal parameters from the sample values, we define the signal to be measured as A sin(21rft) + B , as shown in fig. 1. ¢ and 0 are defined as shown in the figure, ¢ being the location of the second sample with respect to the beginning of the sine wave, and 0 being the interval between samples. LetYo, Yl and Y2 be the values of the three samples. Then YO = A s i n ( Q - 0 ) + B Yl = A sin ~b+ B Y2 = A sin(~b+0) + B . Solving this system of equations yields

~b= tan-1

A =

1 - cos 0

171 2 sin ~b(cos 0 - 1)

172 2 c o s ~ sin0

B =Yl - A s i n ¢ , where 171 =Yo +Y2 - 2Yl and 172 =Yo - Y 2 ' Note that there are two alternative expressions for

107

A, depending on 771 and 172 respectively. The exclusive use of either expression may sometimes lead to loss of significance in the answer due to the quantization level of the analogue-to-digital conversion: in a particular case, for example, 171 may turn out to be very small and thus to have few significant bits. One can guarantee acceptable accuracy by using whichever of 1/1 and 172 is larger in a particular case. In the programme this choice is made on the basis of a comparison of the exponents of the two quantities. The programme checks for two special cases before calculating the phase. (1) If ??2 is exactly zero, which invalidates the equation given above for ~b, one need only set ~ = 7r/2. (2) To avoid misleading phase calculations when the signal is very small or absent, the amplitude is set to zero and the phase is not computed if the sum of the absolute values of the three samples is less than some threshold. The equations derived above could be simplified somewhat if one could assume 0 to be exactly one third of a cycle. However, 0 must be an integral multiple of the clock interpulse interval, so that at high signal frequencies it may be significantly different from a third of a cycle. The requirement that the general forms of the equations be used is made less serious by the fact that 0 changes only when the signal frequency changes, so that sin 0 and cos 0 need not be recalculated as long as the frequency remains the same.

t

A

J

E3

,

t

t

t

Fig. 1. Measurement of amplitude, phase and DC shift of a sinusoidal signal. The three signal parameters (A, ~ and B) and the three sample values (.Vo,Yl and Y2) are shown. The vertical arrows on the time axis represent the instants of sampling.

108

W.R.J. FUNNELL AND C.A. LASZLO

The method as described so far yields the phase of the signal with respect to the instant of the second sample, which is not a particularly valuable piece of information. One will usually wish to measure the relative phases of two signals which have the same frequency. This can easily be done by taking six evenly spaced samples, the first three from one signal and the last three from the other signal, as shown in fig. 2. One calculates the phase for each signal independently, giving q~l and ¢2- Then the phase of the second signal with respect to the first is given by

possibility is the use of the floating-point-arithmetic hardware which is available for the PDP-12 [3]. (The computation time quoted above was obtained using table look-up methods for calculating sines, cosines and arctangents, to speed up the procedure.)

3. Methods

3.1. Introduction The clock being used is the DEC KWl2A [4]. It consists of a variable-rate pulse generator which increments a 12-bit counter with each pulse; the generator is run at a rate of 105 sec -1 in this programme, thus producing an overflow of the counter every 41 msec. The programme can detect the occurrence of a counter overflow. It can also detect an "event", which is defined as the crossing of the threshold of the Schmitt trigger by a signal at the clock's input. The reference signal is connected to one of the clock's inputs, and also to channel 8 of the 10-bit analogueto-digital converter. The second, "object", signal is connected to channel 9 of the converter. The readings taken from the clock and from the analogue-to-digital converter are 12-bit and 10-bit

Aq~ = ~2 - q~l - 30 + 2rr. The speed of the phase-and-amplitude-measurement method presented above is limited by the computation time. As it is implemented at the moment it requires about 30 msec for one calculation of amplitude and phase if the frequency is already known. (To this must be added the time required for the actual sampling.) At the moment all of the computations are done using a floating-point-arithmetic interpreter [2] for ease of programming, and this slows things down considerably. If extra speed were required one could avoid the overhead of the interpreter by using more specialized arithmetic subroutines. Another

~

Reference

I

t

i--.--.J

t

t

t

t

t

~t

Fig. 2. Measurement of relative phase. The phases of the individual signals (~1 and ~2) correspond to ~ in fig. 1. The relative phase is A~. As in fig. 1, the vertical arrows on the time axis represent the instants of sampling.

MEASUREMENT OF TIME-VARYINGSINUSOIDALSIGNALS integers, respectively. They are converted to threeword floating-point format for most of the calculations.

3.2. Frequency The flow diagram for the frequency-measuring subroutine, MSRFRQ, is shown in fig. 3. The computer zeros the counter and starts the clock running. The first event (crossing of the Schmitt-trigger threshold by the reference signal) to occur is detected, and the value of the counter at the time of the event is stored as t 1 . The computer then waits for a second event, which will occur at the same point in the next cycle of the sine wave. If an overflow has not yet occurred the computer goes back to wait for another event, after incrementing n, the number of cycles that

have occurred. If an overflow has occurred, the value of the clock counter at the time of the last event is stored as t 2. This process of waiting for an overflow ensures that the total interval measured will be at least 20 msec, and usually about 40 msec (depending on the signal frequency). At this point the programme checks whether the overflow may in fact have occurred just after the last event (that is, between the machine instruction that detects events and the one that detects overflows) rather than before it. If this is the case, the measurement is still acceptable but the variable t2h must be set to zero. (It was initially set to 1 at the beginning of the subroutine.) This variable is used as a thirteenth bit, to supplement the 12-bit counter value t 2. In other words, the time of occurrence of the last event is taken to be

T 2 = t2h. 212 + t 2 •

Set t~=1~ n=O

1

Finally, the programme converts tl, T 2 and n to floating-point format and calculates the signal frequency according to the formula

Zero counter. Start clock.

1

Await event. When it occurs, tl~counter.

f = ~ T2-t 1 '

J

w h e r e f c is the clock rate (105 sec-1).

q

no+

3.3. Single measurement of amplitude and phase

n--n÷ 1. Await event.

t t ~ c o u n t e r at time of last event.

109

I

If overflow occurred I after last event, I set t=~=O I

1

Ca Iculate frequency I Fig. 3. Flow diagram of the frequency-measurement subroutine, MSRFRQ.

Fig. 4 shows the sequence of operations performed for a single measurement of the amplitude and phase, relative to a reference signal (number 1), of an object signal (number 2) of unknown frequency. (Omitted are the details of the address computations used for storing the sampled values, intermediate results, and the final amplitude and phase.) After the frequency has been measured by the subroutine described in sect. 3.2, subroutine CALC 1 calculates m, 0, and three functions of 0 that are used by CALC2 and CALC3.0 has been defined in sect. 2.2; m is the number of clock pulses between samples. The next step is to take three samples from each channel; this is done by SAM2. CALC2 and CALC3 then calculate the phase and amplitude, respectively, of the first (reference) channel, using the equations derived in sect. 2.2; they are then called again to do

110

W.R.J. FUNNELL AND C.A. LASZLO

MSRFRQ" Measure frequency

CALC1: Calculate m, e0 and various functions of e

J

I

CALC1 SAM2 CALC2(1) CALC3(1)

J

SAM2: Sample 3 paints each from 2 channels

I

CALC2(1): Calculate

I

phase of 1st signal

MEA$2: MSRFRQ

J

CALC3(1): Calculate magnitude of 1st signal

J CALC2(2):Calculate phase of 2nd Si~lnal

-1

MEAS3:SAM2

CALC2(1) CALC2(2) CALC3(2)

Fix & store magnitude and phase Display once

n.o~

I CALC3(2):Calculate

magnitude of 2nd signal

J CALC4: Co,cu,ate J relative phase & magnitude J onvert to deg & dB, fix and store Fig. 4. Flow diagram of the programme for measuring frequency, amplitude and phase once. Each box represents a subroutine. The first one, MSRFRQ, is the one shown in fig. 3. the same thing for the second (object) channel. Finally, CALC4 calculates the amplitude and phase of the second signal relative to the first. (Note that the calculation of the DC shift, B, has not been implemented in the programme as described here. Its addition would be straightforward, but we do not need it for our present application.)

3. 4. Repeated measurements of amplitude and phase Fig. 5 shows the flow diagram for taking repeated measurements of amplitude and phase of a signal of timed frequency. The reference signal is assumed to have constant amplitude. (Again, details of address

Divide by A1. Convert data to dB.

(end) Fig. 5. Flow diagram of the programme for continuously measuring amplitude and phase of a signal of fixed frequency. The various subroutines indicated (MSRFRQ, CALC1, and so on) are those described in fig. 4. manipulation have been omitted.) The first part of the programme, MEAS2, measures the frequency, phase and amplitude of the reference signal, as was described in the preceding section. The routine MEAS3 is then executed repeatedly for as long as one wishes to measure. Each time through it calls SAM2 to take the samples; CALC2 to calculate the phase of the reference signal; then CALC2 and CALC3 to calculate the phase and amplitude of the second signal. It then fixes and stores the magnitude of the second signal; subtracts the phase of the reference signal from that of the second signal; and ffLxes and stores this relative phase. The data are stored in a

MEASUREMENT OF TIME-VARYINGSINUSOIDALSIGNALS fLxed-length buffer; when the buffer is full the programme starts Idling it again from the beginning. The contents of the buffer are displayed once each time through MEAS3. When one has stopped measuring, the programme can go back over the data that have been collected, in order to divide the amplitudes by the amplitude of the reference signal and to convert the amplitudes to decibels. These operations could be done within MEAS3, at the expense of slowing down the measuring loop.

3.5. Programme specifications This programme is written in the DEC LAP6-DIALMS assembler language of the PDP-12; both LINCmode and 8-mode instructions are used. The memory allocation for the programme in an 8000-word machine is shown in fig. 6. The measuring programme proper takes up about 34008 words. In addition it makes direct use of the floating-point interpreter, the

Address ~ctal) 00000 01000 02000

Measuring programme

111

question-and-answer routine, the tables of sines and arctangents, and the data-storage space. The other items shown in fig. 6 (the file handler, monitor, and I/O routines) are parts of a system of programmes to which the measuring programme has access. This system permits the display and editing of the experimental data, and also the formation of a data library. No write-up of the measuring programme is available apart from the present article, although a listing is available which contains a fair number of comments. It is expected that the actual measurement methods will be of more general interest than the details of their implementation within the particular software system that we are using.

4. Conclusion A method has been presented for measuring the frequency, amplitude, phase and DC shift of sinusoidal signals. The method is intended to obviate the necessity of specialized hardware by taking full advantage of the capabilities of the small general-purpose digital laboratory computer. The particular programmes described here are designed to handle audio-frequency signals, but the principles are applicable to sinusoidal signals of any frequency.

03000 040OO

Acknowledgements

050OO Interpreter for 060OO floating-point arithmetic (DEC-08-YQ4B, modified) 07OO0

This work has been supported by the Medical Research Council of Canada.

10000 File handler (DEC-12-FZFA)

References

11000 Table of sines 12000 Question &answer routine (DECUS-12-56) 13000 Monitor 14000 15000 Data storage 16000 Table of arctangents 17000 LAP6-DIAL-MS I/O routines Fig. 6. Diagram of core allocation in an 8000-word machine.

[1] W.R.J. Funnell, The acoustical impedance of the guineapig middle ear and the effects of the middle-ear muscles (Master's thesis, McGill University, Montreal, 1972). [2] Floating-point system (DEC, Maynard, Massachusetts, DEC-08-YQYB-D, 1969). [3] FPP-12 (DEC, Maynard, Massachusetts, DEC-12-GQZAD, 1971). [4] System reference manual (DEC, Maynard, Massachusetts, DEC-12-SRZA-D, 1970).