A resource-efficient adaptive Fourier analyzer

We present a resource-efficient frequency adaptation method to complement the Fourier analyzer proposed by Péceli. The novel frequency adaptation scheme is based on the adaptive Fourier analyzer suggested by Nagy. The frequency adaptation method was elaborated with a view to realizing a detector connectivity check on an FPGA in a new beam loss monitoring (BLM) system, currently being developed for beam setup and machine protection of the particle accelerators at the European Organisation for Nuclear Research (CERN). The paper summarizes the Fourier analyzer to the extent relevant to this work and the basic principle of the related frequency adaptation methods. It then outlines the suggested new scheme, presents practical considerations for implementing it and underpins it with an example and the corresponding operational experience.


Introduction
Recursive estimation methods based on signal models allow the measurement of various parameters or the calculation of different transforms of signals. Their recursive nature renders them particularly suitable for real-time applications. These methods consider the signal under measurement as the output of a hypothetical dynamic system: the signal model. The procedure then consists of estimating the state vector of this system. Measurement systems based on a deterministic signal model are referred to as observers.
Among the observer-based approaches, Hostetter suggested a recursive algorithm to calculate the discrete Fourier transform of a signal [1]. Péceli extended this method to allow calculating arbitrary linear signal transforms recursively [2]. When Péceli's resonator-based observer structure is used as a spectral observer, it is referred to as Fourier analyzer (FA).
These algorithms rely on a linear system model with known parameters. In the case of spectral analysis, however, the frequencies of the harmonic components of the signal may be unknown. Nagy proposed a way of implementing frequency adaptation in the FA. This method, called adaptive Fourier analyzer (AFA), allows locking the resonator structure to the fundamental frequency of the signal in a PLL-like manner, following slow frequency drifts [3,4]. It has been used successfully in several industrial applications, such as a vector voltmeter [3,4], and active noise and distortion cancellation [5,6].   Several modifications to this procedure have been put forward. These frequency adaptation methods, outlined in section 2.2, all require calculating phase differences between complex numbers. Since the complex calculations required by the Fourier analyzer are best carried out in the algebraic form (z = a + jb), phase calculations become computationally intensive. In this paper, we describe a new frequency adaptation scheme avoiding phase calculations at the expense of slower adaptation of the fundamental harmonic frequency, while potentially reducing sensitivity to noise. Section 2 of this paper describes the Fourier analyzer and the basic principle of the frequency adaptation method along with an enumeration of the modifications that have been suggested. Section 3 introduces the new frequency adaptation scheme, while section 4 presents the application which motivated its elaboration, along with practical considerations for the implementation. Section 5 then concludes the paper.

The Fourier analyzer
Péceli's resonator-based observer estimates the state variables of the conceptual signal model [2]. The Fourier analyzer is obtained when using this structure as a spectral observer.
The state variables x i,n of the signal model shown in figure 1 correspond to the complex Fourier coefficients of the signal, therefore, the signal model can be interpreted as a complex multisine generator reconstructing the signal from its discrete Fourier transform (DFT), with each of its state variables (channels) representing a harmonic resonator of the corresponding frequency.
We are assuming the signal to be real-valued. In this case, the conceptual signal model is described by the following equations: Figure 2. Block diagram of the resonator-based observer. Notice how the structure of the observer replicates that of the conceptual signal model shown in figure 1. It is extended by the time-varying g i,n coefficients at the input of each integrator and the common feedback.
where y n is the output of the signal model, that is, the signal to observe. This signal consists of K harmonics, thus it has N = 2K + 1 complex Fourier coefficients, including DC, which compose the state vector x n . Since the signal is assumed to be real-valued, the Fourier coefficients will form complex conjugate pairs: x i,n = x * −i,n . f 0 is the relative frequency of the fundamental harmonic with respect to the sampling frequency: In case of equality, in the equations above, i = −K, . . . , K + 1 and N = 2K + 2.
y n is calculated in every time step by multiplying each x i,n by the corresponding element of the time-varying coupling vector c n to obtain the Fourier components of the signal, then summing them.
The above imply that the complex Fourier coefficients of the input signal can be estimated by an appropriately designed observer. Figure 2 shows the structure of the corresponding observer.
The system equations of the observer are as follows: x n+1 =x n + ge n =x n + g (y n −ŷ n ) , (2.6)

7)
y n = c T nx n , (2.8) wherex n is the estimated state vector,ŷ n is the estimated signal and e n is the estimation error. The vector g n is the product of α N , a freely tunable parameter referred to as observer gain determining the dynamical behavior of the observer, and c n , a vector of complex roots of unity setting the poles and thereby the frequency response of the individual observer channels. If the vectors c n and g n are set according to (2.5) and (2.10), respectively, with α = 1 and f 1 = f s N , a deadbeat observer is obtained which settles in at most N steps and calculates the recursive discrete Fourier transform of the input signal [7].

JINST 11 P10014
The transfer function of channel i in the feedback loop from y n tox i,n · c i,n (from the input of the structure to the adder at the output, see figure 2) can be expressed as [8]: where r i = 1 N , ∀ i in case of the Fourier analyzer.

Frequency adaptation
If the fundamental frequencies of the signal and the Fourier analyzer are different, Nagy suggested a method to estimate the frequency difference from the state variables of the observer and to adapt the frequencies of the resonators to match the signal [3]. The input signal of the observer can be expressed in the following form, as reconstructed from its Fourier series, with a i corresponding to its complex Fourier coefficients: At present, a difference between the fundamental frequencies of the signal and the Fourier analyzer is assumed. They will be denoted by f 0 and f r relative to the sampling frequency, respectively. Accordingly, the signal doesn't necessarily have the same number of harmonics as the number of resonators in the observer.
If f 0 = f r , then after settling, the estimatorsx i are equal to the Fourier coefficients a i . Otherwise, the estimator of the fundamental harmonic in the nth time step can be expressed as follows [8]: wherex 1,i,n is the contribution from the ith harmonic component of the signal in the nth time step to the estimator of the fundamental harmonic. The transfer function H 1 e j2πi f 0 of the fundamental harmonic resonator from y n tox 1,n · c 1,n at the corresponding frequency was defined in (2.11). From (2.13), it follows that argx 1,1,n+1 − argx 1, , then x 1,n provides a good estimate ofx 1,1,n and the error in the frequency of the fundamental resonator relative to the sampling frequency can be expressed as: (2.14) Equation (2.14) can be used as a basis for adapting the resonator frequencies, and the principle of the frequency adaptation scheme stems from it. In the adaptive Fourier analyzer (AFA), Nagy suggested updating f r by 1 N · ( f 0 − f r ) as obtained through (2.14) in every time step. If the change in the frequency of the base harmonic warrants it, the number of resonator pairs K must also be -4 -updated so that K · f 1 < f s 2 ≤ (K + 1) · f 1 still holds, in order to ensure a uniform distribution of the resonator positions up to f s 2 [3]. Several modifications to this procedure have been put forward. Applying a damping factor to the frequency update has been suggested in order to improve noise rejection (robust AFA, RA). Simon and Péceli developed a block-wise frequency adaptation method (block AFA, BAFA) in order to evaluate its convergence properties [8]. With a view to decreasing sensitivity to noise, Ronk suggested extending these approaches by averaging the Fourier components used for the adaptation, thus obtaining the eBAFA from the BAFA [9], and the improved robust AFA (IRA) from the RA [10]. Dabóczi then proposed combining the IRA and the BAFA, yielding the extended improved robust AFA (eIRA) [11].

The zero crossing adaptation scheme
The frequency adaptation methods referred to earlier all rely on calculating phase differences between differently spaced samples of the estimator of the fundamental harmonic, then approximating the frequency difference through numerical manipulation.
We propose estimating the frequency mismatch by tracking when the phase of the fundamental harmonic, argx 1,n , crosses zero. The phase difference between two consecutive zero crossings is 2π, thus if N z zero crossings are detected spanning a time interval t m , an approximation of the frequency error relative to the sampling frequency can be calculated as follows based on (2.14): where the sign of the frequency difference can be determined based on the direction of the zero crossings in phase, as described in section 4.5. The new value of f r can then be used for adapting the frequency of the resonators. This method of calculating the frequency mismatch presents the following advantageous characteristics: • Detecting zero crossings in phase is very efficient in terms of FPGA resources for most practically feasible implementations. For complex numbers represented in polar form (z = A e jϕ ), finding changes in the sign bit of the phase is sufficient. In case the algebraic form (z = a + jb) is used, zero crossings in phase can be identified by considering the signs of the real and imaginary parts simultaneously.
• In case the observed signal and thus the resonator phase are subject to noise, this approach may allow reducing the resulting influence on the estimator of the frequency, as illustrated in section 4.7.

Detector connectivity checks with the Fourier analyzer 4.1 Motivation
Beam loss monitoring (BLM) is an essential method for the protection of particle accelerators. It consists in measuring the showers of secondary particles generated by the primaries escaping from the particle beams under acceleration, impacting on an aperture limit of the machine. BLM systems -5 -are the backbone of the strategy for beam setup and machine protection of all particle accelerators at the European Organisation for Nuclear Research (CERN), one of the world's leading particle physics research facilities. The measurement data acquired by the BLM system are processed in real time and if required for protecting the machines from excessive energy deposition by the beam, a decision maker process running continuously is able to trigger the safe extraction of circulating beams and the inhibition of further injections. Specifically, the rapid response of the BLM system is critical for protecting the machines from short and intense particle losses [12]. Moreover, subsets of the measurement data are displayed in the control room to assist the operators and archived for machine tuning and calibration.
Continuous supervision of the connectivity and correct operation of the BLM detectors, and the functionality of the entire BLM system is therefore essential, yet to our knowledge, no particle accelerator in the world has this feature at present. The current best solution is in operation at CERN's flagship machine, the Large Hadron Collider (LHC), where a connectivity check requiring the accelerator to be offline is enforced on each detector channel every 24 hours. This procedure relies on modulating the bias high voltage of the ionization chambers used as detectors and detecting the corresponding modulation generated in their output current [13,14].
Within the framework of a major accelerator upgrade project at CERN, a new BLM system for the preaccelerators of the LHC is in an advanced state of development. This complex of accelerators is referred to as LHC injectors. In the new BLM system, it is foreseen to use ionization chambers similar to those in the LHC system in most cases, however, other detector types such as diamond detectors, secondary emission monitors and Cherenkov detectors will be installed in some specific locations [15].
Developing a new process continuously supervising the functionality of the BLM system for the injectors is mandatory. Contrary to the LHC, where the cycle of the machine contains long periods without beam while the new fill is being prepared, the operation of the injectors is essentially uninterrupted. This means that any such process needs to run in parallel to the data acquisition without interfering with it.
In the digitized output signal of the ionization chambers, the switching frequency of about 30 kHz of the high voltage power supply (HVPS) and its harmonics have been identified as parasitic components. Detecting these harmonics might allow realizing a non-destructive connectivity check, a highly appealing way of implementing detector supervision without external intervention by simply post-processing the digitized output signal.

The suggested solution
In the BLM system under development for the injectors, each front-end card acquires and digitizes the output current of eight detectors at 500 ksps, using an innovative measurement technique featuring a fully differential integrator [16]. The data is then forwarded to the corresponding backend card over a high-bandwidth optical link for further processing. In order to ensure high throughput and flexibility, the treatment of the digitized data on these cards is based on reconfigurable FPGA devices. This offers a convenient way to implement additional data processing. For the on-line connectivity check, we chose the Altera Cyclone IV FPGA of the back-end module.
-6 -A highly customized version of Péceli's Fourier analyzer [2] was chosen for detecting the HVPS harmonics. It's inherently stable due to the unit negative feedback in the observer, it offers flexibility in locating the resonator poles [17] and a straightforward frequency adaptation method has been proposed to complement it [3], making it possible to track the slowly varying frequency of the HVPS contributions in the signal. The same structure could also be used to detect a low-frequency modulation if desired.

Practical considerations for the Fourier analyzer
We are implementing the resonator-based observer on an FPGA. This means that additions and multiplications have a moderate hardware cost, while divisions and angle calculations are more expensive and thus are better avoided in order to save resources.
According to (2.6)-(2.10) describing the operation of the resonator-based observer, the state update requires additions and multiplications of complex numbers. Using the polar form of complex numbers would be computationally advantageous for multiplication, however, conversions to and from the algebraic form would be required for the additions. Nevertheless, complex multiplications can be realized by a series of multiplications and additions in algebraic form, avoiding the additional FPGA resource use that would result from the conversions between the two forms. Therefore, we implemented our project using this form. We also exploited the fact that the input signal of the observer is real-valued: this implies that the resonators form complex conjugate pairs, thus realizing one of the pair is sufficient and the summed output of the two can be obtained as 2 · Re {·}.
We considered several different solutions for the arithmetic used in the processing, since the values are expected to have a high dynamic range. In MATLAB, we implemented and tested a block floating point architecture based on fixed point arithmetic with a common but variable scaling factor 2 n (n ∈ Z) for all fixed point numbers. However, we discontinued this approach due to its complexity. We first realized the resonator structure on an FPGA with IEEE-754 compliant 32-bit floating point arithmetic using IP cores provided by the manufacturer. In order to reduce FPGA resource use, we then converted this realization to fixed point arithmetic. Our tests showed that at least 26 fractional bits are required for satisfactory performance with data samples acquired from the prototype system when using fixed point arithmetic.
The real and imaginary parts of each state variable and coefficient are treated separately. According to (2.4) and (2.5), the c n coefficients are calculated from their initial values by exponentiation. This calculation can become numerically unstable. Since the c i,n coefficients are unit complex numbers, we circumvented this problem by calculating their phase instead, which can be done by simple addition. Then, the real and imaginary parts of the coefficient required by the processing are obtained in every iteration using memory-based lookup tables storing values of the cosine and sine functions. If required, the precision of this method can be improved by interpolating the values obtained from the lookup tables, but we deemed this unnecessary in our implementation.
We implemented the Fourier analyzer on an FPGA with one single processing channel, which calculates the real and imaginary parts of the new state variable for each resonator sequentially and then the error of the estimator. This processing channel consists of two multipliers and two adders, along with memory blocks storing the resonator values, the phase of the c n coefficients and the lookup tables. The operations related to the state update are parallelizable, which provides a

Tailoring the Fourier analyzer to the application
The output current of the detectors is acquired at a sampling rate of 500 ksps, and since the goal is detecting a signal around 30 kHz and its harmonics, there is little margin in decreasing the sampling rate. In order to ensure continuous, on-line data processing, the observer operates at 500 kHz. If the observer were used for calculating a recursive DFT, with a resolution on the order of 10-100 Hz ensuring proper selectivity, this would imply a prohibitively high number of resonators in the entire Fourier analyzer structure realized on the FPGA we have at hand in our system: in case of a serial implementation, the processing wouldn't be able to cope with the input data flow, and in case of a parallelized implementation, the FPGA resource cost would be too high for the design to fit. We then considered a structure with an extremely coarse frequency resolution, with the DFT points positioned at the HVPS frequency and its harmonics. At the expense of slower convergence, the α N parameter in (2.10) can be used to control the bandwidth of the resonator channels. As shown in figure 3, this can ensure an appropriate frequency selectivity, without impacting on the stability of the processing, while reducing noise in the acquisition [18]. We then turned some of the individual resonator positions into resonator triplets. This means that for some of the lowest order harmonic resonators, we established two additional resonator positions on either side of each of the already existing ones, at frequencies relatively close to but separated from that of the initial resonator. Our aim was to provide estimates of the noise floor along with the measurement of the actual signal component. Figure 4 illustrates how the triplet resonator positions are set with respect to the initial ones, while figure 5 shows how adding the triplets changes the transfer function of a channel. Note that if the frequencies are too close to each other so the bandwidths of adjacent resonator channels overlap, this approach may cause resonances.

Practical considerations for the frequency adaptation
Since the source of the parasitic signal we are aiming to detect is common for all eight channels processed by the same FPGA, we decided to implement frequency adaptation on one channel only and use the same g n and c n coefficients for all other channels as well. For practical use, (3.1) needs some more elaboration. Technically, the phase function also crosses zero when it wraps from π to −π (or vice versa). These events, henceforth referred to as π crossings, must not be considered as zero crossings. Excluding them, however, is very straightforward when using the algebraic form: we only register an event if Im x 1,n changes sign while Re x 1,n > 0. In addition, the formula only gives the absolute value of the frequency difference. To obtain its sign, one also needs to consider the sign of the zero crossings registered: a positive zero crossing in phase is when Im x 1,n changes from negative to positive, which corresponds to f 0 > f r , and the opposite for negative zero crossings. Noise in the input signal may cause spurious zero crossing events as demonstrated in figure 6. These events can be filtered out by requiring a π crossing, that is, a sign change of Im x 1,n while Re x 1,n < 0, to occur between any two accepted zero crossings.
All these considerations assume that the resonator value used for frequency adaptation rotates sufficiently slowly for the adequate detection of zero and π crossings. A necessary condition for this is that the vector under consideration should never make a transition between non-adjacent quadrants of the complex plane. This requires the frequency mismatch to be inferior to f s 4 , a condition automatically fulfilled in most cases. Since we are tracking a slowly varying frequency with a relatively good initial estimate, we don't expect failure cases related to this condition to occur.
The main motivation behind the development of the zero crossing-based frequency adaptation method was to avoid the division required for calculating the phase of the resonators, thus spare the FPGA resources that would be necessary for implementing a divider. However, according to (3.1), this method also requires a division at every adaptation step for calculating the estimator of the frequency difference. Nevertheless, by adapting the frequency whenever there is a suitable pair of zero crossings, the numerator of the fraction will always be 1. In this case, the division can be avoided by constructing a lookup table: for different values of t m , storing the corresponding increment of the phase step size ∆ϕ. The precision of these values could be improved by interpolation between the points of the lookup table if required. When adapting higher order harmonics, the resulting phase estimator values will need to be multiplied by the order of the harmonic. If desired, the phase estimators thus obtained can be averaged over several zero crossings.

Test results
The algorithm with frequency adaptation was tested on a development system operating in the laboratory and a prototype system installed at the Proton Synchrotron Booster (PSB) accelerator.
The frequency of the HVPS contributions varies with time. The main influence on the frequency appears to be temperature, as shown in figure 7. The data were acquired in the summer, when the temperature inside the building and, consequently, inside the rack, tends to follow the exterior -10 -2016 JINST 11 P10014 Resonator frequency HVPS frequency Figure 8. Fundamental frequencies of the resonator structure and of the HVPS contributions. The resonator frequency was sampled at each frequency adaptation step, while samples of the HVPS frequency were taken every 10 minutes. On the latter curve, each marker corresponds to a data sample.
temperature. We see a corresponding cooldown overnight, echoed in the evolution of the HVPS frequency. Figure 8 demonstrates the efficiency of the frequency adaptation in the laboratory system. The adaptation process is capable of following the slow drift of the frequency satisfactorily.
In our implementation, we use α N = 1 50000 , yielding a resonator bandwidth of approximately 5 Hz, which seems to be a good compromise between reduced noise in the acquisition and sufficient bandwidth for signal detection. Figure 9 shows the typical output signal of a resonator properly tuned to the frequency of the HVPS signal. The input signal of the observer is also shown for comparison. The peaks in the signal correspond to actual beam loss events. These events have

Time [sec] Resonator amplitude [bits]
Resonator amplitude Input current value Figure 9. The amplitude of the fundamental resonator tuned to the HVPS frequency, superimposed onto the unprocessed input current values. The peaks in the data correspond to beam loss events, having an influence on the output of the resonator, which recovers quickly nevertheless. a clear impact on the output of the resonator, however, it recovers quickly and the influence of these peaks on the measurement can be eliminated by introducing a blind period after detecting a high peak. The resulting effect on frequency adaptation has been negligible in our acquisitions. In general, the level detected by the resonator properly tuned to the HVPS peak in the presence of HVPS contributions tends to be at least 20 dB higher than the signal in the other two resonators of the triplet, which provides a satisfactory margin for detection.

Merits of the zero crossing frequency adaptation scheme
We implemented the processing described above on an Altera Cyclone IV FPGA device present in the new BLM system. The system clock frequency was 100 MHz.
When using a fully serialized realization in order to minimize resource use, we implemented a functional observer with zero crossing based frequency adaptation using about 1300 logic cells and 21 9-kilobit memory blocks out of the device total of 149760 and 720, respectively. We estimated the resource use of a corresponding "original" adaptive Fourier analyzer (AFA) implementation including a divider at about 2450-2800 logic cells, with 12 9-kilobit memory blocks. At the price of increased memory use, the zero crossing scheme promises a saving in logic cells on the order of 50%.
At our system clock frequency, the typical latency of the divider required for the calculations of the "original" AFA is about 15 clock cycles. If a realization computing all resonators in parallel were implemented, depending on the number of resonators, our observer could process one new input sample in about 15-20 clock cycles. Part of the latency of the divider could be spared, but the "original" AFA would still take approximately 25-30 clock cycles to process a sample. Thus, with our scheme based on zero crossings, we can expect a gain of up to 50% in terms of maximum input sample rate, at the price of increased FPGA hardware use.
Since the estimator of the frequency is only updated at zero crossings of the resonator phase, it may be less susceptible to apparent frequency changes resulting from the influence of noise in the signal. Figure 10 shows an example with a noisy input signal where the zero crossing adaptation Frequency [kHz] Figure 10. Example: frequency estimator yielded by the zero crossing-based frequency adaptation method in comparison to some other adaptation schemes with a noisy input signal. The input signal is a sinusoid with linearly varying frequency from 29.95 kHz to 30.05 kHz contaminated with white Gaussian noise, SNR = 20 dB. In the resonator structure, α = 1. f 1 = 30 kHz, f s = 500 kHz, thus N = 17. In this configuration, the fluctuations exhibited by the frequency estimator obtained with the "original" AFA are about an order of magnitude higher than with the block AFA, therefore the former curve is omitted for clarity. For the robust AFA, we used a damping factor of η = 32. The block AFA was executed with a measurement period P = 333.
scheme yields a considerably smoother frequency estimator than some of the other frequency adaptation methods.

Conclusions
We presented a principle to complement an existing Fourier analyzer implementation with a resource-efficient frequency adaptation method, practical considerations for realizing it and a use case with operational experience. The method promises implementational advantages or reduced noise in the frequency estimator at the price of a less reactive frequency adaptation. It can be implemented efficiently in applications where the computational complexity is critical, especially on synthesized hardware. The performance delivered by the scheme appears suitable for the application we tested it with.