Sensing Frequency Drifts: A Lookup Table Approach

Real-time frequency analysis is a task commonly implemented in last generation signal processing devices. Even if the reliability of these systems is widely recognized, detecting the presence of drifts in the frequency of analog signals is still an open problem. Besides being simple in principle, it needs a suitably long-term memory in order to determine not only in which direction the frequency drifted but also along which specific path. In this paper, we propose a novel strategy based on the concept of jump resonance in nonlinear forced oscillators to determine a digital signal processing paradigm hinged on a lookup table able to provide reliable and real-time information on complex paths of frequency drifts.

suitable control actions, whose efficiency is linked on the drift 23 amplitude and direction [3]. 24 The concept of frequency drift can be further generalized time-series [4] or pandemic related quantities [5]. 30 The associate editor coordinating the review of this manuscript and approving it for publication was Derek Abbott .
The necessity to oversee frequency drifts, however, goes 31 beyond the merely technological field. Let us consider the 32 monitoring of the heartbeat or other biomedical quantities 33 whose frequency is a key indicator to discriminate patho-34 logical states. As concerns the human heartbeat rate, in fact, 35 a fundamental diagnostic information to prevent diseases is 36 gained monitoring its variability [6], which may also occur 37 through complex paths of increasing/decreasing drifts. 38 Detecting frequency drifts is a non-trivial problem from 39 the computational point of view, since it needs to determine 40 the main component of the frequency spectrum of the input 41 signal and to store a significant number of samples to infer the 42 trend. Modern digital signal processing (DSP) units allows 43 to quickly determine the frequency properties of analog sig-44 nals [7], but this feature alone is not sufficient to efficiently 45 infer frequency drifts and complex paths of frequency drifts. 46 In order to provide a simple but effective solution to the 47 problem of real-time detection of frequency drifts either 48 increasing or decreasing, a the theoretical concepts adopted 49 to design a device specifically intended to this aim [8] will 50 be coupled to a fast and efficient algorithm based on the 51 construction of a Lookup-Table (LUT). drift, thus creating an hysteresis in the frequency response. 112 Moreover, the amplitude of the input signal plays a crucial 113 role, since it must be rescaled to a precise value in order to 114 elicit the jump resonance behavior, otherwise, as observed 115 in [21], chaotic oscillations may arise. 116 The device described in [8] is based on the digital 117 implementation of a discrete-time non-autonomous nonlin-118 ear system with cubic polynomial nonlinearity, whose fre-119 quency response displays the jump behavior. It has been 120 shown its efficiency also when the input analog signal con-121 tains more harmonics, but the device allows to discriminate 122 drifts increasing or decreasing from the nominal value of 123 the frequency. The approach presented here is intended to 124 improve the efficiency of the frequency drift detection, thus 125 discriminating among more complex patterns of increas-126 ing/decreasing drifts and consists in the implementation of a 127 LUT based on the peculiar behavior of jump resonance sen-128 sibly reducing both the computational needs and the memory 129 storage. A specifically designed windowed approach allows 130 to detect not only the single frequency drift, but also nested 131 paths of deviations, either increasing or decreasing, thus fully 132 implementing the hysteretic behavior of multijump resonance 133 systems. Moreover, the jump resonance paradigm to design 134 the LUT guarantees the possibility to cope with signals in 135 which direction and amplitude of frequency drifts are not 136 known a-priori. 137 The paper is organized as follows. In Sect. II the basic 138 elements and the design guidelines of jump/multijump res-139 onance systems are outlined, in Sect. III the windowed 140 approach to design and implement the LUT is presented. The 141 experimental setup obtained implementing the LUT in hard-142 ware is discussed in Sect. IV, while the experimental cam-143 paign showing the capabilities of the LUT to detect nested 144 complex path of frequency drifting is presented in Sect. V, 145 together with an applicative example aimed at showing the 146 detection of frequency drifts in the mains frequency of Con-147 tinental Europe power grid. A comparative analysis which 148 allows to determine the main advantages of the proposed 149 approach with respect to different solutions is reported in 150 Sect. VI. Conclusions are drawn in Sect. VII.

152
Let us consider a nonlinear system in which the linear and 153 nonlinear parts can be separated and organized in a feedback 154 scheme, as shown in Fig. 1 where the linear part is represented 155 feedback systems is referred to as Lur'e form [22]. 158 We focus now on the case in which the linear part is a 159 generic second-order transfer function of the type When subjected to an external sinusoidal driving input 169 u(t) = R sin (ω(t)t + ψ) with time-varying frequency ω(t), 170 the output of such system is a nonlinear oscillation whose 171 amplitude at timet depends not only on the instantaneous 172 frequency ω(t) but also on the values assumed by ω(t) for 173 t <t. An example of multijump frequency response is 174 reported in Fig. 2, it can be observed that for a given range 175 of ω, the frequency response is a multi-valued function.

176
The multijump resonance curve in Fig. 2  the frequency continues to grow, the response will undergo a 185 jump down at ω(t) = ω 5 , conversely if ω(t) decreases a jump 186 up at ω(t) = ω 2 can be retrieved in the output amplitude.

187
In this latter case, the upmost branch can be visited either 188 for increasing values of ω(t), leading to a jump down at 189 ω(t) = ω 6 , or for decreasing values of ω to a jump down at 190 ω 1 . Finally, if the frequency of the input signal is swept down 191 from a value larger than ω 6 , a single jump up at ω(t) = ω 3 is 192 observed in the output of the system.

193
It appears evident that such kind of resonance, and its 194 hysteretic behavior, allows to discriminate among complex 195 paths of increments/decrements of the frequency of the given 196 input signal by simply inspecting the jumps occurring in the 197 amplitude of the output signal. Our approach is based into 198 transferring such behavior, suitably rescaling the frequency 199 ranges, into a simple and effective LUT.

200
To design the linear and nonlinear parts of the feedback 201 oscillator, so that jumps around a specific value of the 202 frequency, it is useful to adopt an approach based on the 203 describing function method [23]. The method is based on 204 two hypotheses: (i) the system variables can be approxi-205 mated as the sum of a finite number of sinusoidal terms, 206 and (ii) a low-pass filtering effect occurs in the loop. Pro-207 vided that these hypotheses are verified, the nonlinear term 208 can be approximated by a function of the amplitude of its 209 input N (U ). Notwithstanding the effective approximation 210 introduced by this method, the describing function approach 211 proved to be effective also in presence of higher harmonics, 212 including chaotic oscillations [24]. Assuming ψ = 0, and indicating with (ω) and (ω) the 224 real and imaginary part of G −1 (jω), the modulus of Eq. (4) 225 can be written as: and thus rearranging terms that, defining now X = U 2 , it can be written as Since jump resonance occurs when the frequency response 236 is a multi-valued function of ω, in order to determine the 237 conditions on the parameters to obtain a frequency response 238 VOLUME 10, 2022 shaped as in Fig. 2, it necessary that for a range of ω Eq.  2 have a remainder equal to zero. Therefore 251 it follows: where the quotient is a third order polynomial on which the  modelling. The extraction of the values is simply obtained by 292 using a straight line between two points and stepping them 293 by a reasonable amount in order to minimise MCU memory 294 footprint. However, in the proposed approach, an important 295 concept which is referred afterwards as indexer function can 296 be used to model more complex scenarios, where a large 297 number of steps is required or a more complex pattern has 298 to be implemented, instead of simple memory arrays.

299
Each of the piecewise approximation parts is afterwards 300 referred as frame which has two boundaries, i.e. the maxi-301 mum and minimum x-coordinates for the segment.

302
Notice that the intervals are not strictly contiguous and 303 there are overlapping parts. This is the key to the jump 304 mechanism implementation and to the hysteresis behavior, 305 that is translated into a current frame pointer. The current 306 frame pointer, which models the state of the LUT, is the dis-307 criminant for the jump management logic. A query (lookup) 308 to the table is matched on the boundaries first, to discern 309 if a jump is needed. In the occurrence of a jump, the logic 310 is recursively applied to each of the frames, by moving the 311 current frame pointer. Once the queried value falls in between 312 the boundaries, the algorithm falls back to standard LUT 313 • jump to the previous frame and re-run the query • jump to the next frame and re-run the query The query re-run step implicitly means moving the current 338 frame pointer. Using this approach, the hysteresis behaviour 339 is implemented. When a new frequency has been computed, 340 a query is done using the LUT manager entity. The manager 341 has a state where the current frame of the LUT is stored 342 and the frame has the two boundaries (lower and upper).

343
If the queried value is inside the current frame, stored in 344 the manager, the value is directly taken as a regular LUT. In the following sections the aforementioned terms are 361 explained in depth, in order to clarify the implementation 362 logic and strategy.

363
A. THE FRAME

364
The approach foundation is the frame, which is a LUT with 365 ancillary data. A frame is composed by: Key to the understanding of the whole strategy is the 397 current frame pointer. This pointer is read when the LUT is 398 accessed, to match search item with the frame boundaries, 399 and written when a lookup refers an off-boundary request and 400 VOLUME 10, 2022 a new frame must be searched into. The search-and-reassign strategy is recursively applied across frames, consolidating 402 into the last, value-resolving call. configurations (e.g., pits or loops).

412
Once the instantiation process has been fulfilled, the LUT 413 can be queried through the manager entity lookup method.

414
The manager algorithm can be synthesized into the following

451
The adopted solution is based on a mean value detection plus 452 a zero-crossing strategy. The first one ensures that eventual 453 FIGURE 7. Experimental setup: a standard PC equipped with a GPIB interface is used to drive the Agilent 33250A generator imposing a given frequency drift path; the signal is conditioned using standard OP-AMP configurations to fit the MCU amplitude input spefications; the LUT implemented on the STM32F446RET6 microcontroller provides an analog output indicating the occurrence of jumps which is visualized on an oscilloscope.
DC components lead to erratic results while the latter simply 454 counts the times the signal crosses the mean value during 455 an interval slice. To enhance the mean-crossing strategy and 456 reject false crossing flags, signal is low-pass filtered inside 457 the same detection loop.

459
The previously described system was fully implemented on 460 a hardware platform, to assess its validity and performances, 461 as schematically represented in Fig. 7. A series of choices 462 were done during the implementation transposition, which 463 are discussed below.  To test and evaluate the algorithm implementation correct-479 ness and integrity, the MCU board was fed using an Agilent 480 33250A function/arbitrary waveform generator, an oscillo-481 scope was also connected to watch the DAC (output). To drive 482 a frequency sweep to stimulate all the possible hysteresis 483 paths, the 33250A generator was driven using custom devel-484 oped software (WinForms C GUI) and GPIB which, in turn 485 used UART debug output from the MCU to watch the current 486 detected frequency and current frame pointer extracted from 487 the running algorithm. The complete system was thoroughly tested by using the 496 whole toolchain described. Especially, using the specifically 497 designed GUI and setup, every frequency path was fed to 498 the system, by ensuring the right path on the curve was 499 followed. Recorded data was compared with evidence on the 500 instruments displays to assess the whole setup validity.

502
Using the previously described setup, a testbench was devel-503 oped to assess and validate the strategy. Essentially, the evalu-  to the specific green slope. In the right column ( Fig. 8(b),   Notice also that as per LUT design, the only allowed path to 538 reach the uppermost values must follow this path.

539
The third experiment follows a path which highlights the 540 possibility to form loops by falling from the third window 541 back to the first. This happens thanks to the fact that no 542 constraints in the LUT enforce the previous/next window 543 relationships. With a similar mechanism, pits can be created, 544 allowing a state-machine like behaviour. This occurs when 545 the frequency of the input signal increases above the first 546 jump value and then decreases producing a jump up around 547 28 kHz and a further jump down at 25 kHz.

548
The last path shows a reverse path following, from the 549 highest to the lowest frequencies values mapped in the table. 550 The reverse sweep is again highlighting the loose relationship 551 between LUT windows. A single jump occurs in this case 552 around 30 kHz.

553
Clearly, more complex paths of frequency drifts can occur 554 in real cases, however they can be reconstructed simply look-555 ing at the series of jumps in the output generated by the LUT. In order to test the LUT approach outlined in this paper, 559 we focused on detecting the occurrence of frequency drifts 560 in an historical dataset reporting the mains frequency of the 561 Continental Europe power grid taken with a sampling time of 562 one minute in the week 9th−15th of September 2012 [27]. 563 The LUT has been realized designing the multivalued fre-564 quency response so that the nominal value of the power grid 565 frequency, i.e. f n = 50 Hz, lies within the jump resonance 566 range and the jump frequencies are located at: f 1 = 49.98 Hz, 567 f 2 = 49.99 Hz, f 3 = 50.01 Hz, f 4 = 50.02 Hz, f 5 = 50.04 Hz, 568 f 6 = 50.06 Hz. This is obtained choosing ω 0 = 270 rad/s 569 and ξ = 0.01. However, due to the modularity of the LUT 570 implementation, note that the jump frequencies can be fixed 571 by suitably rescaling the frequency ranges in the previously 572 adopted LUT.

573
The output of the LUT, reported in Fig. 9 together with 574 the trend of the drift of the mains frequency f around the 575 nominal value f n , allows to determine the specific path of 576 frequency drift by inspecting its jumps and the levels among 577 which they occur. The output of the LUT, in fact, shows a 578 series of jumps which unveils specific path of drifts in the 579 mains frequency. In order to clarify this point, let us refer to 580 the zoom reported in Fig. 10 where two jumps are highlighted. 581 The first jump occurs between sample 7200 and sample 7201, 582 corresponding to the jump down at f 1 from the upper branch 583 of the multivalued frequency response. This indicates a drift 584 from a frequency above the nominal value (f = 50.02 Hz) to 585 a frequency below it (f = 49.98). Moreover, the second jump 586 still indicates that the frequency is decreasing but, as a jump 587 up is observed at f 3 this indicates that the frequency remains 588 above the nominal values (drifting from f = 50.02 Hz to 589 f = 50.01 Hz).

591
The dynamic behaviour of the LUT approach described in 592 this paper allows the creation of loops and initial/final states, 593 thus allowing the creation of complex state machines which 594 are based on the input values. 595 VOLUME 10, 2022 FIGURE 8. Experimental results: Test 1: (a) frequency of the input signal sweeping from 100 Hz to 90 kHz and (b) oscilloscope trace of the LUT output; Test 2: (c) frequency of the input signal sweeping from 100 Hz to 33 kHz, then from 33 kHz to 28 kHz and then from 28 kHz to 90 kHz and (d) oscilloscope trace of the LUT output; Test 3: (e) frequency of the input signal sweeping from 100 Hz to 33 kHz, then from 33 kHz to 100 Hz and (f) oscilloscope trace of the LUT output; Test 4: (g) frequency of the input signal sweeping from 100 kHz to 100 Hz and (h) oscilloscope trace of the LUT output.
The proposed LUT approach has several advantages 596 over other approaches. Leaving out the frequency detec-597 tion portion of the system, which can be implemented 598 using several strategies, the hysteresis path behaviour 599 can be obtained using one of the below described 600 algorithms.    As concern the computational cost of the proposed solu- prescaler gives the ability to span from (virtually) 629 0 to the maximum ADC frequency (2.4Msps in the 630 STM32F443RE). The implemented system is working 631 on a single ADC channel, thus the sample rate also 632 depends on the number and the operating mode of 633 channels that are acquired (e.g. 7.2 Msps in triple 634 interleaved mode). 635 2) The number of bits of the ADC, which translates in 636 clock cycles to obtain a sample. In the specific setup, 637 the ADC requires 15/13/11/9 cycles for 12/10/8/6 bits, 638 respectively.

639
3) The DMA buffer size, where the frequency detection 640 routine is taking place. Considering that the frequency 641 detection is essentially a threshold passing count to 642 have a reliable detection, a sufficient number of traver-643 sal points must be found in the buffer. When working 644 at higher ADC sample rates, problems may arise with 645 low frequency signals, as the needed buffer size must 646 be increased to become able to see a sufficient number 647 of traversals.

648
4) The DMA buffer size also directly influences the drift 649 detection speed. Given that the detection algorithm is 650 executed over the whole buffer to compute the fre-651 quency, an increase in its size directly translates in an 652 increase of the time between two successive frequency 653 computations.

654
The choice of these parameters determines, there-655 fore, the effective bandwidth. The clock prescaler puts a 656 limit on the maximum detectable frequency: the lower the 657 prescaler, the higher the sample rate and, consequentially, the 658 higher the maximum detectable frequency. However, having 659 a higher sample rate requires a larger DMA buffer size to 660 contain one or more input waveform cycles, in order to 661 properly detect the frequency. This means that with higher 662 sampling rates, the low frequency detection capability is 663 capped by the DMA buffer size. 664 Moreover, the varying the number of ADC bits allow a 665 slight clock cycle improvement, thus an higher sample rate, 666 at the cost of the loss of resolution in the signal. This can be 667 considered a mild issue when only the frequency detection 668 algorithm must be executed on the signal.

669
Summarizing, higher sample rates, obtained by means 670 of lower prescaler values or decreased ADC bits resolu-671 tion, translate into higher maximum detectable frequency 672 but requires a DMA buffer size increase to allow also low 673 frequency detection; larger DMA buffer sizes allow increased 674 span of detectable frequencies but decreases the detection 675 speed.

676
Carefully tuning both the parameters, the frequency span 677 and resolution can be adjusted, to suit the specific detection 678 needs. In the experimental analysis presented above, we fixed 679 the ADC clock at 11.250 MHz, with a resolution of 12 bits, 680 and DMA buffer size of 4096 samples. This leads to admissi-681 ble input signals with frequency in the range between 10 Hz 682 and 300 kHz.    In fact, it arises also in mechanics and civil engineering, 737 where vibrations and their frequency must suitably moni-738 tored. Moreover, a recently emerging area of interest is also 739 in the context of the pandemic monitoring. As a specula-740 tive example, in the Covid-19 pandemics, vaccination rate 741 assumed a fundamental role in controlling the dimension of 742 the outbreak. Sensing the variation of the vaccination rate, 743 and the effect that this produces, can be considered as a 744 frequency drift detection problem. By monitoring the vacci-745 nation rate, it may be possible to qualitatively infer the impact 746 of the outbreak of new SARS-CoV-2 variants.

747
The paradigm introduced in this paper offers an easy and 748 MCU friendly way of implementing a device for frequency 749 drift detection based on a nonlinear dynamical feedback 750 system not designed to exclusively monitor electrical quan-751 tities. This key feature allows to obtain a device which is 752 highly reconfigurable, as opportunely shaping the frequency 753 response of the nonlinear system, it is possible to implement 754 complex LUTs with hysteresis. Moreover, reconfiguring the 755 nonlinear part of the oscillator, increasing the order of the 756 polynomial nonlinearity easily ensures the possibility to get 757 the desired selectivity for the drift detection, allowing to 758 discriminate over complex paths of increase and decrease 759 of the frequency without the a-priori knowledge of the drift 760 dynamics. It should be noted that the only a-priori knowledge 761 needed to construct the LUT is the nominal frequency of 762 the signal under monitoring. Slightly tuning the parameters 763 of the nonlinear system ensures the capability to adapt the 764 LUT to the specific characteristics of the frequency drifts, 765 thus optimizing its capabilities.

766
The frequency range over which it is possible to design 767 the desired nested detection hysteresis is limited by several 768 factors which, however, make the solution suitable up to 769 3 MHz signals for the considered STM32F446RET6 micro-770 controller. Its ease in terms of algorithm complexity, renders 771 its implementation immediate even on power-limited devices. 772 The possibility to implement the LUT on low-cost devices is 773 a fundamental advantage of the proposed scenario since puts 774 this new strategy in the class of cheap DSP equipment.