IWAVE -- An Adaptive Filter Approach to Phase Lock and the Dynamic Characterisation of Pseudo-Harmonic Waves

We present a novel adaptive filtering approach to the dynamic characterisation of waves of varying frequency and amplitude embedded in arbitrary noise backgrounds. This method, known as IWAVE, possesses critical advantages over conventional techniques making it a useful new tool in the dynamic characterisation of a wide range of data containing embedded oscillating signals. After a review of existing techniques, we present the IWAVE algorithm, derive its key characteristics, and provide tests of its performance using simulated and real world data.


I. INTRODUCTION
The co-inventor of the MASER [1], Arthur Schawlow, is said to have advised his students: "Never measure anything but frequency!"There exist a great variety of techniques for measuring the frequencies of pseudo-harmonic waves, broadly addressing two different classes of problems.
In the first class, the average wave characteristics are estimated across a measurement interval, under the implicit assumption that the wave properties are essentially static, or when changes in these properties over the measurement interval are not of interest.This problem is addressed by a wide variety of methods, including Welch's method using discrete Fourier transforms (DFTs) [2,3], Pisarenko's method [4], MUSIC [5], and ESPRIT [6].
In the second class, the oscillator is constantly evolving, and we seek a time-evolving best estimate of oscillator parameters.This problem is most often solved using phase locked loops (PLLs) [7], or their hardware realisation, lock-in amplifiers [8].The myriad applications of PLLs in science and engineering include, for example, a recent proposal for digital PLLs as an alternative technology for the readout of the photodiode clusters used to stabilise the alignment of mirrors in future gravitational wave detectors [9].
In searches for almost continuous wave (CW) signals in gravitational wave detectors, the evolving oscillation of continuous wave (CW) signals is at a very low signalto-noise ratio (SNR) and so the frequency evolution of such signals in ground based gravitational wave interferometers cannot be inferred from the raw data as is necessary for successful tracking with a conventional PLL.Instead, coherent matched filtering techniques, such as the F-statistic used in CW searches [10,11], are able to achieve higher sensitivity to these weak signals at the price of substantially greater computational burden.Stack-slide-based semi-coherent algorithms expedite the computation to some extent, at the cost of sensitivity, by summing the signal power in multiple coherent segments after sliding the segments in the frequency domain to account for the signal phase evolution [12][13][14][15].More efficient semi-coherent methods, e.g., signal tracking algorithms based on hidden Markov models [16][17][18][19], have been developed to tackle the computational challenge as well as to allow for some uncertainties in the signal evolution model [20][21][22].
Several other techniques beyond the conventional PLL have been developed for a range of applications outside the field of gravitational wave research, though these algorithms are not adapted for the detection of the very weak CW signals expected from gravitational wave sources.One example is the second order generalised integrator (SOGI-PLL) [23].Developed for the problem of characterising the characteristics of AC voltages in power lines, SOGI was one of the first methods to successfully addresses the problem of efficient generation of a copy of an input sinusoid that is out of phase (the so-called quadrature or Q phase) with the input signal.Other more recent papers, for example [24][25][26][27], have developed the SOGI algorithm for a range of practical applications.A second example is the enhanced PLL (EPLL) [28,29], which solves the problem of tracking the amplitude of a harmonic wave in addition to its frequency.Further PLL developments are well summarised in [30].
The IWAVE technique described in this paper is a new type of PLL addressing the dynamic characterisation of evolving pseudo-sinusoidal signals.Unlike a conventional PLL, the adaptive element is a filter rather than an oscillator or counter.IWAVE has certain advantages over existing PLLs.Firstly, IWAVE produces benign output when the PLL is unlocked.In the case where the output of the PLL is being used to control something, for example, some parameter of a gravitational wave detector in a closed loop feedback system, then, in colloquial terms, IWAVE does no harm when it isn't working.IWAVE naturally tracks the amplitude as well as frequency using a single feedback loop, unlike EPLL which requires two loops internally for the amplitude and the phase.IWAVE is initialised using a small set of free parameters corresponding directly to physical oscillator propertiesjust an initial frequency and a single time constant; other PLL algorithms typically contain many control parameters that do not have a clear physical meaning.IWAVE also has the ability to characterise, simultaneously, multiple oscillations having almost-degenerate frequencies using the cross-subtraction method described in Section III C.This last advantage has led to detailed analysis of almost-frequency-degenerate violin modes of fused silica suspension wires in advanced LIGO [31].Comparisons of IWAVE with the SOGI and EPLL algorithms are given in Appendix A.
As we discuss in Section III B, IWAVE as currently implemented is not sensitive to signals where the ratio of the amplitude of the target wave to the root mean square (RMS) of the noise background is less than about 0.3.When applied to broadband gravitationalwave data, therefore, IWAVE cannot be expected to detect gravitational wave CW signals at the strengths anticipated for sources described in the literature.However, preprocessing steps to divide the data into narrower frequency bands, thereby significantly reducing the RMS of the noise, may yield promising search methods.Studies of these ideas are under investigation, and will form the subject of future papers.The potential advantage is that, by using the data itself to track the evolving frequency of the oscillation, IWAVE may be significantly less reliant on banks of templates for wave evolution, and hence may require significantly less computational resources than existing CW search methods.For now, however, IWAVE represents a simple, well-characterised and useful technique for analysing quite complex spaces of evolving oscillators at relatively low computational burden, which is already being applied to studies of important background oscillations in gravitational wave data.This, then, is a methods paper describing how IWAVE works, how it performs, and giving some examples of that performance on gravitational wave data.
The IWAVE algorithm has many potential applications beyond gravitational wave science.In the control of brushless electric motors, IWAVE could out-perform standard vector control methods at low rotation rates where the back-emf in the motor windings is weak [32].In radio communications, IWAVE might be used to separate closely spaced channels with frequency evolution and multipath splitting [33].In heart magnetometry, IWAVE could be used to estimate Fourier coefficients of the cardio-magnetic signal using an electrocardiograph signal to modify the phase evolution as the heartbeat rate evolves [34].In physics IWAVE could be applied to atomic force microscopy [35], or as part of an optical squeezing scheme for the readout of interferometers [36].The authors look forward with anticipation to seeing what other applications we have not noticed.
This paper consists of: a description of the IWAVE method in Section II; a discussion of the limits of applicability of IWAVE in Section III; and an overview of certain other PLL methods in Appendix A. Space constraints have led us to leave out many mathematical steps; a full treatment of the mathematics can be found at [37].A software library implementing IWAVE in C with wrappers into MATLAB and PYTHON/NUMPY is available on a public git repository here [38].

II. THE IWAVE METHOD
A. The Core Algorithm Before writing down the IWAVE core algorithm, we consider as a starting point the Z-transform [39] of a regularly sampled time series, x p , where p increases with time, and Ω = w − i∆ has real part w the reciprocal of one e-folding for the weighting of previous samples and imaginary part ∆ equal to the frequency of the Z-transform component in radians per sample.
Z n (Ω) obeys an iteration equation Unit amplitude phasor input, x n = e in∆ , to Equation 2 results in output Z n = e in∆ /(1 − e −w ), which has zero phase shift with respect to the input and a larger amplitude dependent on w.Scaling the input by a factor of 1 − e −w leads to an iteration algorithm which passes phasors at frequency ∆ with unit gain and zero phase shift y n = e −w e i∆ y n−1 + (1 − e −w )x n. ( This iteration algorithm is the core of IWAVE.As we shall see, it responds resonantly at frequency ∆.In the language of signal processing, Equation 3 is an infinite impulse response (IIR) filter because it generates its n th output using the current input x n , the previous output y n−1 and a two-input, two-output, multi-input, multioutput (MIMO) filter; y n and x n are, in general, complex variables.We will also need the real representation of the transfer function for IWAVE derived from Equation 3 by writing x n = x R n +ix I n , y n = y R n +iy I n , and y n−1 = z −1 y n , where z −1 is the sample delay operator.In these terms, Equation 3 can be re-written as where the elements of the transfer function matrix are In this form, the transfer function is seen to be second order in sample delay.We determine the response of the core algorithm to an input consisting of a phasor of arbitrary frequency Θ radians per sample by substituting x n = e inΘ into Equation 3. The response is also a phasor, y n = Ae i(nΘ+Φ) , where A (Θ) and Φ (Θ) are the frequency dependent magnitude and phase lag of the output phasor with respect to the input given by Thus, phasors of arbitrary frequency are eigenfunctions of the core algorithm with eigenvalues Ae iΦ .The resonant character of the eigenvalues at frequency ∆ can be seen in the denominator of Equation 5, which is identical to the resonant denominator in the SOGI filter [23].
< l a t e x i t s h a 1 _ b a s e 6 4 = " 0 4 j m  1, so that we are in a limit where the frequency is in the vicinity of a narrow resonance, the magnitude of the filter output can be approximated as so that the peak is approximately Lorentzian in shape with full width at half maximum (FWHM) of 2w radians per sample.Using the sampling rate, f s , in Hz we give other properties of narrow resonances occurring when w 1 in Table I, where τ s is the sampling period in seconds, τ is the response time in seconds and ∆ 0 is the resonant frequency in radians per sample.

B. Application of IWAVE to a Real Sinusoidal Input
We next discuss the usual case where the input data is a real oscillation at the IWAVE resonant frequency, x n = cos (n∆) .We decompose x n into two phasors, x n = x f + x b , where x f = e +in∆ /2 and x b = e −in∆ /2, each of which is an eigenfunction of the core algorithm.Substituting the frequencies ±∆ into Equation 6 and rearranging, we obtain the response where By inspection, the locus of y n is an ellipse in the complex plane having semi-major and semi minor axes 1 + a and 1−a. Figure 2 shows the result of driving the IWAVE core algorithm with an input x n = cos n∆ for n ≥ 0. The output starts at the origin, spiralling outward towards a limiting ellipse in the steady state.Notice that the argument of the output y n is always the same as the phase, n∆, of the input sinusoid.
. Symbols for IWAVE acting on two component phasor or real inputs.The parameters w and ∆ may be adjusted to reflect changes in the state of the drive.
Input sinusoids at frequencies other than ∆ also result in an elliptical steady state output, though with different inclination angles and eccentricities, and with smaller overall areas due to the falloff in the magnitude of the response for phasors having frequencies far from ∆.
A matrix transformation can be used to transform the elliptical locus into a circular one, with the real part of this circular locus being a sinusoid in-phase with the input wave, and the imaginary part having the same amplitude but lagging the input wave by 90 • .We refer to these in-phase and out-of-phase components as the D and Q phases, respectively.This transformation can be expressed as a sequence of three elementary operations on the vector whose elements are the real and imaginary parts of y n : a rotation through an angle of φ 2 about the origin; shears parallel to the real and imaginary axes by factors of 2 1+a and 2 1−a respectively; and finally a rotation through an angle of −φ 2 about the origin.These three matrices can be combined into a single transformation (10) The IWAVE core algorithm followed by this matrix transformation results in the output of both D phase and Q phase copies of the input drive.Thus, IWAVE is an example of what is referred to in signal processing parlance as an orthogonal state generator.Our convention follows other papers in the field in that the D (Q) phase output is in (out of) phase with the input at resonance.As we shall see, the Q phase quadrature can be used to generate an error signal to detect changes in frequency, allowing IWAVE to be used in place of a reference oscillator in a PLL.The sum in quadrature of the two phases, , is an estimate of the input signal amplitude.We will use the symbols in Figure 3 to denote the application of IWAVE either to complex phasor or real sinusoidal inputs.
We next consider the transfer functions from a real sinusoidal drive, an example of which is shown in Figure 4.Note that the dependence on frequency off-resonance is different for the D and Q outputs.The D transfer function rises linearly in frequency below the resonance, and falls linearly in frequency above it, having a phase lead of 90 • below the resonance and a phase lag of 90 • above it.The Q transfer function has a flat frequency   < l a t e x i t s h a 1 _ b a s e 6 4 = " R B N p k 5      response below the resonance but falls as f −2 above it, and is in phase with the drive below the resonance, but 180 • out of phase with the drive above it.
This behaviour is similar to that of a driven series RLC tank circuit, where the transfer functions from the input voltage to the voltage across the capacitor and resistor are similar to those between the input and the Q and D phase outputs, respectively.The relatively light suppression of low frequency off-resonance signals in the Q phase output can be important, particularly in cases where the quality factor Q of the circuit is set low by using a relatively large w coefficient.As in the resonant circuit, the ratio of the resonant to low frequency response is Q.
Changes in the amplitude of the incoming wave at the resonance result in a corresponding change in the quadrature sum, A n = D 2 n + Q 2 n , but with a response time τ = τ s /w leading to a single pole in the response to amplitude changes at s = −1/τ.Figure 5 shows the response at IWAVE's resonance frequency to an input having a sinusoidally modulated amplitude for a variety of response times.

C. An IWAVE-based Phase Locked Loop
In order to phase lock with IWAVE, we need a measure of departures in frequency from the frequency ∆ of the harmonic wave at the input.We achieve this by exploiting the response time τ of the IWAVE algorithm.Consider a harmonic wave initially at frequency ∆.IWAVE yields both a D n phase and a Q n phase copy of this wave at its output.The product of the out-of-phase copy and the input wave, A 2 cos (n∆) sin(n∆) = A 2 /2 sin (2n∆), is a pure harmonic signal at frequency 2∆.Now, consider an input signal where the wave develops anomalous phase and amplitude disturbances, δ and ε, respectively, so that x n = A(1 + ε) cos(n∆ + δ).For elapsed times significantly less than τ following the onset of these disturbances, the outputs D n = A cos n∆ and Q n = A sin n∆ are unaffected by them.Physically, where the oscillator has an unchanged frequency and amplitude, its complex plane representation precesses in a circle about the origin, and the current IWAVE filter output therefore provides a predictor of the point where such a static oscillator will land next sample.If the frequency or amplitude evolve, then the combinations n are estimates of the departure of the actual evolution of the complex coordinate of the oscillator state from this static model.These combinations can be written in matrix form as follows: This equation shows that E n and F n each consists of a static offset plus upper sidebands of frequency 2∆, linear in the phase and amplitude offsets, respectively.The upper sidebands, however, appear as components of a rotating phasor in the space spanned by the E n and F n signals.We remove them using the complex carrier version of IWAVE, subtracting its outputs from its inputs.We have determined experimentally that using twice the w factor in the 2∆ IWAVE filter yields an acceptable error signal for the detection of phase departures.
Figure 6 is a schematic for the full IWAVE-based phase and amplitude detector.We have only considered here the case where the input data is real and we are tracking harmonic waves.The case where the input data is two-component rotating phasors is simpler, as it can be shown that the combination contains a pure DC signal in phase offset with no upper sideband contamination.Also shown is the feedback path from the filtered error signal, δφ n , back to ∆ through an integrator, discussed below.
The response of IWAVE to modulation of the frequency of the carrier is shown in Figure 7. Knowing the analytic form of the frequency response, we can analyse the closed loop IWAVE PLL to determine an appropriate choice of feedback gain, G.A schematic for the feedback controller is shown in Figure 8.Any difference between the incoming wave frequency and the IWAVE filter central frequency results in an accumulating phase shift in the homodyne detector.This accumulation of phase is represented by the factor of 2π/s, where s = 2πif , with f being the signal's frequency.The response of IWAVE to phase, confirmed by the measurements underlying Figure 7, acts on the accumulated phase to produce the error signal.As shown in Figure 8, the feedback path from the error signal to a correction in the central frequency of IWAVE consists of an adjustable gain, G, an addition of the gain boosted error signal to the previous value of the phase shift per sample and a scale factor to convert from radians per sample to frequency in Hz.The closed loop gain is calculated by the usual consistency argument around the loop, from which we obtain the closed loop transfer function, This is the transfer function of a driven damped harmonic oscillator.We want a critically damped response since the control signal will be used to measure the oscillator frequency.For critical damping we require two coincident real poles, which is achieved if G = τ 2 s /(4τ 2 ).The closed loop transfer function then takes the simpler form The response to frequency or phase modulation is therefore flat below the knee frequency of 1/(4πτ ) Hz, where due to the 2 poles it has rolled off to −6 dB, and drops as 1/f 2 above that frequency.Larger values of G result in sharper turnover or a resonant peak, corresponding to the underdamped case.

III. LIMITS OF APPLICABILITY OF IWAVE
The time constant, τ , determines the responsiveness of IWAVE to changes in wave frequency and amplitude, as well as the noise bandwidth of the filter.Decreasing τ results in a faster response and better ability to stay locked on waves whose frequencies are changing, but also increases the bandwidth of IWAVE to input noise, resulting a noisier error signal.The optimal τ is small enough so that IWAVE stays locked when the frequency changes, but not so small that excessive background noise is admitted by the filter.Section III A is on frequency tracking, and Section III B discusses locking in the presence of additive noise and the character of the error signal.

A. Frequency Tracking
Consider a wave whose displacement at time t is h (t) = A cos (2πf (t) t), so that the frequency of the wave is changing.The ability of IWAVE to track the evolving wave depends on the value of the response time parameter, τ .Figure 9 shows the results of running IWAVE on a swept sinusoidal wave starting at 20Hz and increasing linearly in frequency at a variety of rates.At each sweep rate, a variety of values of τ were used.The sum of the squares of the deviation between the actual frequency and the frequency reconstructed by IWAVE over a time interval of 20 seconds, here referred to as χ 2 , was calculated as a measure of the accuracy of frequency reconstruction.
The value of χ 2 is seen to be a function of τ , with a pronounced minimum that is a function of the sweep rate.The χ 2 statistic also becomes noisy at both ends of the range of values of τ .Here we explain the smooth descent and ascent in χ 2 on either side of the minimum, and obtain a formula for the optimim value of τ .We also explain the onset of noise at either end of the range of τ .
Starting at the minimum, χ 2 rises with τ because of the response time of the closed loop transfer function of the servo, from Equation 14.This transfer function is that     of two RC lowpass filters in series, each having an exponentially decaying impulse response of timescale 2τ .The response time of the entire transfer function is therefore the time duration of the autocorrelation of this exponentially decaying function.This autocorrelation initially rises linearly in time after the impulse, before reaching a maximum at time 2τ and then decaying exponentially.The time between the impulse and the point where the exponential decay reaches 1/e of its maximum value is 6τ .This causes the frequency tracking of iwave to lag behind that of the input swept sine wave, leading to a frequency discrepancy at any given time of ∆f 1 = 6τ (df /dt).
Below the minimum of χ 2 (τ ), the rise in χ 2 with decreasing τ is explained by the frequency response of the core IWAVE algorithm, and is best understood by the analogy between the IWAVE algorithm and the characteristics of a damped harmonic oscillator discussed in Section II B. At higher values of damping, the frequency response is maximal at a frequency f below the natural frequency of the undamped oscillator, f 0 by an amount 0.01 0.10 1.00 R X e P O m 9 e O / e x 7 K 0 5 B U 9 p / A H 3 u c P v t q N V w = = < / l a t e x i t > 0.01 < l a t e x i t s h a 1 _ b a s e 6 4 = " E + P c o 0 A n q l K d g P U n g 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " t m 6 j S T y w b C + 0 F N 0 9 5 L 9 6 7 9 z F v X f H y m S P 4 A + / z B 9 P g j p g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " m Q < l a t e x i t s h a 1 _ b a s e 6 4 = " D P j s 7 T T q A v A 4 P n 1 N A E J 3 U r 1 q / q h 6 9 B I v g q S R F 1 G P B i 8 c K 9 g P a U D a b T b t 0 s 4 m 7 k 0 I J + R 1 e P C j i 1 R / j z X / j t s 1 B W x 8 M P N 6 b Y W a e n w i u 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 o w 3 4 J T l g g p g K 7 f g j c G v q l q U a 5 g = " > r F e l 6 1 r 1 m r m B H 7 A e v s E u 5 S N Z w = = < / l a t e x i t >

2.5
< l a t e x i t s h a 1 _ b a s e 6 4 = " 0 y r y j Figure 9. Measured χ 2 as defined in the Section III A from a simulation where an input swept sinusoid plus additive white gaussian noise was fed into IWAVE.A range of sweep rates from 0.1 to 2.5 Hz s −1 were used.For each injected wave, IWAVE was run with a range of τ between 0.01 s and 5 s.Also overlaid are three curves in bold.The curve labelled 'measured' is the minimum of χ 2 (τ ) for each df /dt.The curve labelled 'theory' is the calculated prediction for the minimum of χ 2 (τ ) from Equation 15.The curve labelled U.L. is the predicted upper limit on τ above which loss of lock is predicted as discussed in the final paragraph of Section III A, given by Equation 16, versus the measured χ 2 at that τ .The thick black line follows the values of τ where IWAVE loses lock, as demonstrated by the onset of noise in χ 2 , over simulations having a variety of df /dt values.
given by the relationship f 0 = f 2 + 1/ (2π 2 τ 2 ).This causes a systematic offset between the frequency, f , returned by IWAVE, and the frequency, f 0 , of the input signal.The square of this frequency difference is an additional contribution to the χ 2 statistic.By squaring and adding the frequency discrepancies arising from these two effects, and minimising with respect to τ , making the assumption that 2π 2 τ 2 (df /dt) 2 1, we arrive at a value of τ that minimises the χ 2 statistic for frequency tracking.
Figure 9 shows, in bold red the values of τ opt corresponding to the measured minimum τ and χ 2 across each of the χ curves, vs. the measured χ 2 at that minimum, along with, in bold blue, the value of τ opt , versus the value of χ 2 at τ opt .There is good agreement between the theoretical τ opt and the value determined from simulations, within 29% for the largest value of df /dt studied, and within 13% for the smallest one.At larger values of df /dt, τ opt is smaller, so there is more broadband noise in the error signal, and a more sophisticated optimisation on τ would lead to a larger optimal value than that predicted by Equation 15.
We next discuss the breakdown of IWAVE at high values of τ , where the χ 2 (τ ) curves become noisy, indicating loss of lock.The following argument leads to successful prediction of the value of τ where this breakdown occurs.Consider a wave whose displacement at time t is h(t) = A cos (2πf (t)t), so that the frequency of the wave is changing.Assume that IWAVE is locked at time t, so that the IWAVE output is m(t) and is equal to h(t).At time t + τ , the wave displacement has evolved to h(t + τ ) = A cos (2π(f + df )(t + τ )).In the time interval [t, t + τ ] the IWAVE output has not had time to respond to the frequency shift df , and therefore takes the form m(t + τ ) = A cos (2πf (t + τ )).The phase shift between the wave and the IWAVE output at time t + τ is ∆φ 2πτ df .The error signal for IWAVE is approximated by the integral h(t )m(t ) dt over time interval [t, t + τ ].If the phase shift between the incoming wave and the IWAVE output over this time interval is greater than π, the error signal will undergo a sign change causing loss of lock.Therefore, a condition for IWAVE to remain locked is that ∆φ < π, or τ ∆f < 1/2.Writing df = τ × df (t)/dt, we arrive at the upper bound that τ should obey for IWAVE to remain locked.This line is drawn on the χ 2 curves in Figure 9 in bold black, labelled U.L., versus the measured value of χ 2 at that value of τ in each simulation.
The limit tracks the onset of lock loss, as demonstrated by large excursions in χ 2 , well, across different trial values of df /dt.

B. Response to Noise
The error signal for the IWAVE PLL is derived from the product of the input data stream and the Q phase output of the IWAVE filter, minus the product DQ of the two iwave outputs.Where the wave frequency is static, this subtraction removes the upper sideband component at frequency 2f .When the frequency changes, this causes an additional transient upper sideband component, which is removed using a second IWAVE filter at frequency 2f having a response time half that of the primary IWAVE filter.Finally, we divide by the square of the wave amplitude, because both the amplitude of the incoming wave, and the amplitudes of both IWAVE outputs scale linearly with the wave amplitude.We need to divide this scale out otherwise the PLL loop gain will be dependent on the wave amplitude.This is a form of homodyne detector.Because the error signal incorporates the unfiltered wideband input, the error signal incorporates the broadband noise of the incoming data.The distribution of the error signal therefore reflects the spectral characteristics of the input data over the full Nyquist band.If, for example, the incoming data includes a time domain transient with a broad spectral distribution, this transient will be reflected in the time history of the error signal from IWAVE.If the out of band noise is sufficiently large in amplitude, then IWAVE will lose lock.
We have determined experimentally that at the optimal choice of response time, τ opt given in Equation 15, then IWAVE will stay locked when the ratio of the wave peak amplitude to the root mean square noise amplitude exceeds 0.3.Future work to improve the performance of IWAVE at lower signal-to-noise ratios could involve, for example, prefiltering the input data to focus on a narrower frequency band about the frequency of interest for the waves under study.The noise content of the error signal, which leads to noise also in the estimate of the IWAVE frequency, is greater at smaller values of τ where the bandwidth of the IWAVE filter resonance is larger.At sufficiently small values of τ , the incursion of noise leads again to loss of lock.This can be seen in Figure 9 for τ < 0.02.The optimal value of τ is affected by higher noise levels, with larger values than that given in Equation 15 becoming optimal.
The exact value of τ where lock loss occurs and the effects of noise on the optimal value of τ could be determined by a detailed stochastic differential equation analysis, and is beyond the scope of this paper.However, a simple scaling argument can be made.Noise in the error signal leads to a random component being added to the phase shift per sample.This means that the reconstructed frequency, which also forms the servo control signal, contains a component that undergoes a random walk, and hence grows with the square root of the number of samples.This means that you might expect the onset of loss of lock to occur at a τ which scales as one over the square of the signal to noise ratio, so that IWAVE works at τ above a lower limit that goes down by a factor of two when the signal to noise ratio is enhanced by a factor of √ 2. This was verified by injecting additive white Gaussian noise on top of the swept sine wave, and noting that the threshold for lock loss at low τ reproduced this predicted behaviour.
In terms of noise in the error signal, the most significant source is the product of noise in the input data, since this noise enters the error signal without bandpassing through the IWAVE filter.However, this broadband noise has an RMS amplitude independent of the amplitude of the tracked wave.If the tracked wave drops in amplitude, but the RMS of the broadband noise at the input does not, then the noise component of the error signal will scale inversely proportional to the amplitude of the line.This is exactly what is seen when IWAVE is run in practice on real data.The normalisation of the error signal with the amplitude squared is necessary to ensure that the feedback loop gain is independent of the wave amplitude.
The other function of the error signal is to provide an indication of whether or not the PLL is locked.With a non-stationary RMS, this is difficult.However, if we scale the error signal by multiplying it by the amplitude of the wave being tracked, this stabilises the RMS of the error signal noise at the level of the line amplitude.If we further divide by the long-term RMS of the input data, we then obtain a statistic that has a stable RMS of order unity when the servo is locked.Departures from lock manifest themselves as large transient spikes of amplitude greater than 10.This effect will be seen in the discussion of performance on gravitational wave data in Section IV.In particular, in Figure 11, the error signal plotted for each of the four harmonics in the study has been scaled in this way.

C. Use of multiple IWAVE filters in parallel
Multiple IWAVE filters can be applied to multiple harmonic waves in a single data stream.However, when those harmonics are frequencies spaced closely together, this causes crosstalk between the different filters.All harmonics present in the data will enter every IWAVE instance through the error signal.To mitigate this effect we employ a cross-subtraction scheme, as illustrated in Figure 10.
The schematic shows only two IWAVE instances, but the technique generalises to any number of filters.Assume that the two IWAVE filters are locked on line frequencies f 1 and f 2 .The D and Q phase outputs from each IWAVE instance are fed into a phase shifter which estimates the wave signal one sample in the future at the frequency of the wave tracked by the filter.This wave sample is stored until the next input data sample, when it is subtracted from the input data to the other filter.In this way, the input data to each IWAVE filter is purged of the harmonic being followed by the other filter.This technique has been determined to successfully lock multiplets of up to 20 filters, leaving the frequency estimates from each filter free of oscillations at the difference between frequencies in the multiplet, the effect seen in the absence of the cross subtraction technique.The technique works with IWAVE because the two quadrature outputs can be used to generate an output shifted through an arbitrary phase shift, and because the outputs are at the same amplitude as the input wave onto which IWAVE is locked.

IV. PERFORMANCE OF IWAVE ON REAL WORLD DATA
We present an example of the application of IWAVE to a set of harmonic noise components of gravitational wave data taken from the LIGO open data centre web site [40].The data was acquired by the LIGO Hanford interferometer on November 30 th 2016 and consists of 800 seconds of calibrated strain data from the Hanford interferometer during an 800 second lock stretch.The data was preprocessed with a fourth order Butterworth highpass filter at 30Hz, followed by four third order Chebyshev type 2 bandpass filters between 5Hz and 300Hz, applied in series.Finally, a fine adjustment to the spectrum was made using a single real pole at 10Hz.The resulting data is dominated by the 20-80Hz band and is approximately white in that range.It is not a requirement that the input data to IWAVE be whitened, but if a spectral feature is to be successfully tracked by IWAVE, its peak should rise above the noise floor in the surrounding background; whitening ensures that this is the case.Eight harmonic features were identified from a broad power spectrum and were tracked using eight parallel IWAVE instances.These originate from various instrumental sources present in the Hanford instrument at the time when the data was acquired.Violin mode harmonics have been studied in detail by Cumming et al. [31].
The results at four of the identified frequencies are shown in Figure 11.For each harmonic, the frequency, amplitude and scaled error signal (as described in Section III B) are displayed.Amplitudes are in dimensionless strain units, so 1 represents the 4km length of the LIGO detector arms.Though the lines are in some cases almost degenerate in frequency, there is no evidence of beats between the reconstructed frequencies.The error signal is roughly static with approximately unit RMS, though the distribution is non-Gaussian because of the modulation of the input noise by the sinusoidal IWAVE output.A non-statistical transient fluctuation in the modified error signal in the 36.7 Hz line at around 360 s corresponds to a jump in the IWAVE reconstructed frequency, indi- < l a t e x i t s h a 1 _ b a s e 6 4 = " j v X H l i T < l a t e x i t s h a 1 _ b a s e 6 4 = " S y n 9 g I L G U I P e y h R H q C y F X 9 N s 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " J h n W 0 a u w A H g z J + w r j + F q 0 5 J 5 s 5 h j 9 w P n 8 A j h u N O w = = < / l a t e x i t > (i) < l a t e x i t s h a 1 _ b a s e 6 4 = " b c 1 0 q V S E 9 7 4 E y 6 T 1 K B k i 0 V h K o i J y e x v 0 u c K m R F j S y h T 3 N 5 K 2 J A q y o x N p 2 B D 8 J Z f X i X N a s W 7 q l z e V 0 u 1 s y y O P J z A K Z T B g 2 u o w R 3 U o Q E M B v A M r / D m C O f F e X c + F q 0 5 J 5 s 5 h j 9 w P n 8 A k q q N P g = = < / l a t e x i t > (l) Figure 11.The results of IWAVE frequency tracking on four pseudo-harmonics present in data from the LIGO Hanford detector.Subfigures a-c, d-f, g-i and h-l pertain to harmonics at the nominal frequencies displayed at the top of each group of three sub-plots.At each frequency, the reconstructed frequency, amplitude and scaled error are plotted.Refer the discussion in Section III B for a description of error signal scaling.Some of the harmonics have frequencies that exhibit significant time evoultion, others are more static.In a few cases, IWAVE can be seen to have lost and re-acquired lock.
In the case of the 36.7 Hz wave, at about 350 s, a lock loss can be seen in frequency, and also via a spike in the scaled error; similarly at about 50 s in the 37.3 Hz wave, and at about 680 s in the 46.1 Hz wave.
cating that IWAVE momentarily lost lock when either the frequency of this sinusoid shifted rapidly or IWAVE jumped between two almost degenerate harmonics.A second loss of lock can be seen in the 37.3 Hz data at around 30 s accompanying a sudden drop in the amplitude of the harmonic.The eight IWAVE instances all successfully tracked their target harmonics, with the reconstructed error signal providing a useful performance indicator.

V. SUMMARY AND FUTURE WORK
We have described IWAVE, a novel orthogonal state generator, resonant filter and phase locked loop for the dynamic tracking of harmonic waves.The method has a single input parameter, its response time.The algorithm has a low computational load, so that many harmonics can be tracked in real time using a single CPU core.The ability to track multiple closely-spaced harmonics means that the method lends itself well to applications where there are dense 'forests' of harmonics, such as in LIGO violin mode clusters, and communications applications.IWAVE has been applied to LIGO strain data and used to study the character of violin modes in ultra low loss fused silica suspensions.There are many possible applications of the IWAVE method.It is complementary to existing PLL algorithms that we have described in the Appendix.We have supplied software implementations of IWAVE in C with MATLAB and PYTHON wrappers to encourage the community to find other applications and uses [38].
The authors can see several directions in which IWAVE could be improved.The error signal is susceptible to broadband noise contamination, and a narrower band alternative would be of benefit in applications with very weak signals, though narrowbanding will reduce the responsiveness of the method to frequency changes.There are also applications where feedback is not important, for example, the use of IWAVE for novel resonators that is promising for resonant detectors of weak signals in physics, such as those of dark matter axions [41].We look forward to seeing what the community finds to do with our harmonic tracking algorithm.

VI. ACKNOWLEDGMENTS
The members of the Sheffield gravitational wave research group wish to acknowledge the support of the Science and Technology Facilities Council under grants ST/V005693/1.ST/V001752/1, ST/V001744/1, ST/V001019/1 and ST/R000336/1.One of us, (I.J.H.) was supported by the Hollows Scientific Foundation.L.S. acknowledges the support of the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), Project No. CE170100004, the United States National Science Foundation, and the LIGO Laboratory.M.F.acknowledges the support of the Fonds de la Recherche Scientifique-FNRS, Belgium, under grant No. 4.4501.19.This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org/), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector.Additional support for Advanced LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the United States National Sci-ence Foundation, and operates under cooperative agreement PHY-1764464.Advanced LIGO was built under Grant No. PHY-0823459.Additional support for Advanced LIGO was provided by the Australian Research Council.Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain.
detection, for example.The SOGI orthogonal state generator does not have the same transfer functions as the IWAVE one and this is not surprising given the different methods by which they are obtained, although the denominators of the two transfer functions are identical, as they both represent a resonant response to a harmonic drive.

EPLL
The EPLL, although apparently seeming to differ from quadrature signal generation-based PLLs, such as SOGI-PLL, is closely related thereto [28].A conventional PLL is embedded within a second feedback servo that uses the integral of the quadrature phase of the PLL output, multiplied by the unintegrated quadrature phase to synthesise a 2ω signal, which can then be subtracted from the input data, compensating for the 2ω component of the phase detector output inside the PLL.Further, the common amplitudes are also equal to half the amplitude of the sinusoid in the input data.In the EPLL structure, the input data are normalised to the PLL by this amplitude, thereby removing the amplitude dependence of the phase detector [42].
Unlike IWAVE, there is a necessity for two servos to be locked at once, and several numerical parameters that must be adjusted, including parameters necessary to implement the s-plane design digitally.EPLL is susceptible to interference from other waves present in the input and to DC offsets.The latter two issues are addressed by prefiltering the input to the EPLL.EPLL has been widely implemented because of its comparative simplicity, ability to track wave amplitude and the availability of code implementing the method [43].

Figure 1 .
Figure 1.Magnitude (a) and phase (b) of the response of IWAVE to phasor input as a function of phasor frequency in radians per sample, for different values of w and ∆ = 1.257.A smaller w results in a sharper resonant peak.

Figure 1
Figure1shows A and Φ as a function of Θ for ∆ = 1.257 and four values of w: 1, 0.1, 0.01, 0.001.Starting from Equation6and writing δ = (∆−Θ) 1, and w 1, so that we are in a limit where the frequency is in the vicinity of a narrow resonance, the magnitude of the filter output can be approximated as t e x i t s h a 1 _ b a s e 6 4 = " n S E D m V 6 z 7 Z F C b 0

2 Figure 2 .
Figure 2. Output of the IWAVE core algorithm for an input xn = cos n∆ for ∆ = 0.3068, 0 ≤ n < 27 and w = 0.3.The points are the individual outputs yn and the solid line is the limiting ellipse discussed in Section II B.
t e x i t s h a 1 _ b a s e 6 4 = " 3 U Y f T n y n 3 q w r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k R w b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t o 5 T x b D r 5 J 2 v e Z d 1 i 7 u 6 p V G N Y + j C C d w C l X w 4 A o a c A t N a A G D I T z D K 7 w 5 w n l x 3 p 2 P R W v B y W e O 4 Q + c z x + E r I 0 4 < / l a t e x i t > (b) r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k R w b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t o 5 T x b D 4 h S p 4 c A U N u I U m t I D B E J 7 h F d 4 c 4 b w 4 7 8 7 H o r X g 5 D P H 8 A f O 5 w + D J 4 0 3 < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " S S P 4 i 4 B 6 l M e y n Q 3Y M C g J t M Z R 4 I s = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k p S R D 0 W v H i s a G u h D W W z n b R L N 5 u w u x F K 6 E / w 4 k E Rr / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k R w b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t o 5 T x b D 4 h S p 4 c A U N u I U m t I D B E J 7 h F d 4 c 4 b w 4 7 8 7 H o r X g 5 D P H 8 A f O 5 w + G M Y 0 5 < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " T r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k R w b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t o 5 T x b D

Figure 4 .
Figure 4.An example of the transfer functions between a real sinusoidal input and the D and Q phase outputs.Here the sampling rate was 256 Hz, the resonant frequency was 1 Hz and the filter quality factor Q was set to 12.3.Notice that the low frequency attenuation in the Q phase output is about 0.08, which is 1/Q.Subfigures a and c (b and d) are the magnitude and phase of the D phase (Q phase) output.

Figure 5 .
Figure 5.An example of the response of IWAVE to amplitude modulated signals.The carrier frequency was 60 Hz and the modulation depth was 10%.Modulation frequencies in the range 33 mHz ≤ fAM < Hz were used and results are shown for five different values of τ .The 3dB point for turnover between flat response and proportionality to 1/fAM is f 3dB AM = 1/ (2πτ ) in each case.Subfigures (a) and (b) are the magnitudes and phases of the responses, respectively.

Figure 6 .
Figure 6.A schematic of the IWAVE PLL.The error signal after upper sideband filtering is δφn.The switch closes the feedback loop to adjust ∆ and 2∆ for the two IWAVE instances.Use of the complex input IWAVE for attenuation of the 2∆ component in the error signal reduces the computational load compared with the use of a real input IWAVE on the En signal alone.
r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k R w b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t o 5 T x b D 7 Y y o G e l l b y b + 5 3 V T E 1 7 7 G Z d J a l C y x a I w F c T E Z P Y 3 G X C F z I i J J Z Q p b m 8 l b E Q V Z c a m U 7 I h e M s v r 5 J 2 v e Z d 1 i 7 u 6 p V G N Y + j C C d w C l X w 4 A o a c A t N a A G D I T z D K 7 w 5 w n l x 3 p 2 P R W v B y W e O 4 Q + c z x + E r I 0 4 < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " R B N p k 5 I E

Figure 8 .
Figure 8.An s-plane model of the IWAVE PLL.
K b 5 7 y X r x 3 7 2 P R W v D y m W P 4 A + / z B 8 9 R j p U = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " j K M 2 6 V s 9 F f E s W F 8 g H n O o N F + N T q 4 = " > A A A B 7 X i c b V B N S w M x E J 3 1 s 9 a v q k c v w S J 4 K p s q 6 r H g x W M F + w H t W r J p t o 3 N J k u S F c r S / + D F g y J e / T / e / D e m 7 R 6 0 9 c H A 4 7 0 Z Z u a F i e D G + v 6 3 t 7 K 6 t r 6 x W d g q b e a d T y O I p w A q d w D g F c Q g N u o A k t I P A I z / A K b 5 7 y X r x 3 7 2 P R W v D y m W P 4 A + / z B 9 J b j p c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U I l 5 B C o 8 E 5 a X X 8 j M B m z 9 W 9 4 I U Q w 6 G P R u o b y m R P 4 A / T 5 A 9 V l j p k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Q M 3 y c T I S a M k h j i e D J r X V C B W M 9 5 4 = " > A A A B 7 X i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e w W U Y 8 F L x 4 r 2 A 9 o 1 5 J N s 2 1 s N l m S r F C W / g c v H h T x 6 v / x 5 r 8 x b f e g r Q 8

l e 4 Q 0 p 9 I 2 <
L e 0 c e i t Y D y m W P 4 A / T 5 A y S F j s w = < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " 5 o Z 8 R + u C z 4 1 A n 4 p e m n i f I 0

4 < l a t e x i t s h a 1 _ b a s e 6 4 =
e n N d l a 8 5 Z z R y D H 3 D e P g H L 9 o 1 y < / l a t e x i t > 0." x L X o 8 2 y J 5 o

1 7 d 9 a 6 5 < l a t e x i t s h a 1 _
r l Y a / j K M I R 3 A M p + D B O T T g C p r Q A g p D e I A n e H a E 8 + i 8 O K + L 0 o K z 7 D m E H 3 D e P g G Z L 4 1 P < / l a t e x i t > 0.b a s e 6 4 = " D

n 9 g f f 4 A
b m 6 S f A = = < / l a t e x i t > df dt < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 n E y 6 n B 3 B Z 9 2 V I U m S h y Y U S b P x / w = " > A A A B 6 n i c d V D J S g N B E K 2 J W 4 x b 1 K O X x i B 4 G n q S m I m 3 g B e P E c 0 C y R B 6 O j 1 J k 5 6 F 7 h 4 h D P k E L x 4 U 8 e o X e f N v 7 C y C i j 4 o e L x X R V U 9 P x F c a Y w / r N z a + s b m V n 6 7 s L O 7 t 3 9 Q P D x q q A A A B 6 n i c d V D L S s N A F L 2 p r 1 p f V Z d u B o v g K i T p w 7 o r u H F Z 0 T 6 g D W U y n b R D J 5 M w M x F K 6 C e 4 c a G I W 7 / I n X / j 9 C G o 6 I E L h 3 P u 5 d 5 7 g o Q z p R 3 n w 8 q t r W 9 s b u W 3 C z u 7 e / s H x c O j t o p T S W i L x D y W 3 Q A r y p m g L c 0 0 p 9 1

8 <
y s i Y y w x 0 S a d g g n h 6 1 P 0 P 2 l 7 t l u z 3 Z t K q e G t 4 s j D C Z z C O b h w A Q 2 4 h i a 0 g M A I H u A J n i 1 u P V o v 1 u u y N W e t Z o 7 h B 6 y 3 T 8 f 4 j W 8 = < / l a t e x i t > 0.l a t e x i t s h a 1 _ b a s e 6 4 = " e k k 2 s E 2P N B Z K c j q z e K G G e 3 Z Y 0 b 0 = " > A A A B 6 n i c d V D L S s N A F L 2 p r 1 p f V Z d u B o v g K i R t a J J d w Y 3 L i v Y B b S i T 6 b Q d O n k w M x F K 6 Ce 4 c a G I W 7 / I n X / j 9 C G o 6 I E L h 3 P u 5 d 5 7 w p Q z q S z r w y h s b G 5 t 7 x R 3 S 3 v 7 B 4 d H 5 e O T t k w y Q W i L J D w R 3 R B L y l l M W 4 o p T r u p o D g K O e 2 E 0 6 u F 3 7 m n Q r I k v l O z l A Y R H s d s x A h W W r q 1 T H 9 Q r l i m 6 9 d c z 0 G W W a 8 6 j u d q Y j u e 7 9 e Q b

9 <
0 m 7 a t p 1 0 7 5 x K o 3 q O o 4 i n M E 5 X I I N L j T g G p r Q A g J j e I A n e D a 4 8 W i 8 G K + r 1 o K x n j m F H z D e P g H x 6 I 2 M < / l a t e x i t > 0.l a t e x i t s h a 1 _ b a s e 6 4 = "

5 <
2 5 p z V z D H 4 A e f t E 7 / I j W o = < / l a t e x i t > 1.l a t e x i t s h a 1 _ b a s e 6 4 = " / i 3 B Y h z x j h v 7 f 6 N I D W j W S g t E X o M = " > A A A B 6 n i c d V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 9 j E m n o s e P F Y 0 d p C G 8 p m u 2 m X b j 7 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N m 7 a C i j 4 Y e L w 3 w 8 y 8 I B V c a Y w / r N L K 6 t r 6 R n m z s r W 9 s 7 t X 3 T + 4 U 0 k m K W v T R C S y G x D F B I 9 Z W 3 M t W D e V j E S B Y J 1 g c l n 4 n X s m F U / i W z 1 N m R + R U c x D T o k 2 0 o 1 r 4 0 G 1 h u 2 6 W / e w i 7 D t e Y 5 3 f m

Figure 10 .
Figure 10.A schematic showing the cross subtraction technique for two parallel IWAVE filters.The feedback loop components are omitted for clarity; each IWAVE filter is separately instrumented as shown in Figure 6.
t e x i t s h a 1 _ b a s e 6 4 = " 7 f 9 H l h J 3 l G / w u n 83 5 V h 1 o N Y D 8 H 8 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B q M Q L 2 E 3 i H o M e P E Y 0 T w g W c L s p D c Z M j u 7 z M w K I e Q T v H h Q x K t f 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G t z O / 9 Y R K 8 1 g + m n G C f k Q H k o e c U W O l h z K 9 6 B V L b s W d g 6 w S L y M l y F D v F b + 6 / Z i l E U r D B N W 6 4 7 m J 8 S d U G c 4 E T g v d V G N C 2 Y g O s G O p p B F q f z I / d U r O r d I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b + 5 3 V S E 9 7 4 E y 6 T 1 K B k i 0 V h K o i J y e x v 0 u c K m R F j S y h T 3 N 5 K 2 J A q y o x N p 2 B D 8 J Z f X i X N a s W 7 q l z e V 0 u 1 s y y O P J z A K Z T B g 2 u o w R 3 U o Q E M B v A M r / D m C Of F e X c + F q 0 5 J 5 s 5 h j 9 w P n 8 A g f O N M w = = < / l a t e x i t > H N E c 6 L 8 + 5 8 L F p z T j Z z D H / g f P 4 A g 3 i N N A = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " B X M e b m 6 g e 9 z I 9 e Z 2 F M a E o g J S d R M r O r d I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b + 5 3 V S E 9 7 4 E y 6 T 1 K Bk i 0 V h K o i J y e x v 0 u c K m R F j S y h T 3 N 5 K 2 J A q y o x N p 2 B D 8 J Z f X i X N a s W 7 q l z e V 0 u 1 s y y O P J z A K Z T B g 2 u o w R 3 U o Q E M B v A M r / D m C O f F e X c+ F q 0 5 J 5 s 5 h j 9 w P n 8 A h P 2 N N Q = = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " B p s U m u y k 0 H W m e g 0 P D A I l Q E 5 7 A P w e y t i I y w w s T Y d E o 2 B G / 5 5 V X S r t e 8 q 9 r l f b 3 S O M v j K M I J n E I V P L i G B t x B E 1 p A I I J n e I U 3 h z s v z r v z s W g t O P n M M f y B 8 / k D i x G N O Q = = < / l a t e x i t > (g) < l a t e x i t s h a 1 _ b a s e 6 4 = " V 7 m i d 9 g G G 3 r P S H U 3 g d p 7 Y 4 7 l r O r d I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b + 5 3 V S E 9 7 4 E y 6 T 1 K B A 2 p p g x t O g U b g r f 8 8 l / S r F a 8 i 8 r 5 b b V U O 8 n i y M M R H E M Z P L i E G t x A H R r A Y A B P 8 A K v j n S e n T f n f d G a c 7 K Z Q / g F 5 + M b j 6 C N P A = = < / l a t e x i t > (j) < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 0 H/ m k N 7 B y N d U I W 5 x d F O 2 q M K e a Y = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B q M Q L 2 E 3 i H o M e P E Y 0 T w g W c L s Z J I M m Z 1 d Z n q F s O Q T v H h Q x K t f 5 M 2 / c Z L s Q R M L G o q q b r q 7 g l g K g 6 7 7 7 e T W 1 j c 2 t / L b h Z 3 d v f 2 D 4 u F R 0 0 S J Z r z B I h n p d k A N l 0 L x B g q U v B 1 r T s N A 8 l Y w v p 3 5 r S e u j Y j U I 0 5 i 7 o d 0 q M R A M I p W e i i P L 3 r F k l t x 5 y C r x M t I C T L U e 8 W v b j 9 i S c g V M k m N 6 X h u j H 5 K N Q o m + b T Q T Q y P K R v T I e 9 Y q m j I j Z / O T 5 2 S c 6 v 0 y S D S t h S S u f p 7 I q W h M Z M w s J 0 h x Z F Z 9 m b i f 1 4 n w c G N n w o V J 8 g V W y w a J J J g R G Z / k 7 7 Q n K G c W E K Z F v Z W w k Z U U 4 Y 2 n Y I N w V t + e Z U 0 q x X v q n J 5 X y 3 V z r I 4 8 n A C p 1 A G D 6 6 h B n d Q h w Y w G M I z v M K b I 5 0 X 5 9 3 5 W L T m n G z m G P 7 A + f w B k S W N P Q = = < / l a t e x i t > (k) < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 C r a + o 9 y h u M + s A H T W I B Z N 3 D 5 3 P 8 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B q M Q L 2 E 3 i H o M e P E Y 0 T w g W c L s p D c Z M j u 7 z M w K I e Q T v H h Q x K t f 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G t z O / 9 Y R K 8 1 g + m n G C f k Q H k o e c U W O l h 7 K 4 6 B V L b s W d g 6 w S L y M l y F D v F b + 6 / Z i l E Ur D B N W 6 4 7 m J 8 S d U G c 4 E T g v d V G N C 2 Y g O s G O p p B F q f z I / d U r O r d I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b + 5 3