Rapid parameter determination of discrete damped sinusoidal oscillations
Jim C. Visschers, Emma Wilson, Thomas Conneely, Andrey Mudrov, Lykourgos Bougas
RRapid parameter determination of discrete damped sinusoidal oscillations
Jim C. Visschers, Emma Wilson,
2, 3
Thomas Conneely, Andrey Mudrov, and Lykourgos Bougas ∗ Institut f¨ur Physik, Johannes Gutenberg Universit¨at-Mainz, 55128 Mainz, Germany Department of Mathematics, University of Leicester, Leicester, LE1 7RH, UK Photek Ltd., St Leonards on Sea, East Sussex, TN38 9NS, UK (Dated: October 23, 2020)We present different computational approaches for the rapid extraction of the signal parametersof discretely sampled damped sinusoidal signals. We compare time- and frequency-domain-basedcomputational approaches in terms of their accuracy and precision and computational time requiredin estimating the frequencies of such signals, and observe a general trade-off between precision andspeed. Our motivation is precise and rapid analysis of damped sinusoidal signals as these becomerelevant in view of the recent experimental developments in cavity-enhanced polarimetry and el-lipsometry, where the relevant time scales and frequencies are typically within the ∼ − µ sand ∼ −
100 MHz ranges, respectively. In such experimental efforts, single-shot analysis withhigh accuracy and precision becomes important when developing experiments that study dynamicaleffects and/or when developing portable instrumentations. Our results suggest that online, running -fashion, microsecond-resolved analysis of polarimetric/ellipsometric measurements with fractionaluncertainties at the 10 − levels, is possible, and using a proof-of-principle experimental demonstra-tion we show that using a frequency-based analysis approach we can monitor and analyze signals atkHz rates and accurately detect signal changes at microsecond time-scales. I. Introduction
Precise and rapid signal-parameter estimation is im-portant for both fundamental and applied research, andbecomes particularly crucial when observing and con-trolling fast processes in real time (e.g., chemical reac-tions), and in the development of portable instrumenta-tion where fast, real-time, data streaming and inspectionis essential.Several different research fields rely on the preciseand accurate extraction of the time constants and fre-quencies of damped sinusoidal signals. Prominent ex-amples include: nuclear magnetic resonance (NMR) [1],where information on the structure and the spin envi-ronment of a target molecule is extracted from precisedetermination of the frequency and decay constant of adamped sinusoidal signal; free-induction-decay (FID) op-tical magnetometry [2–7], where the magnetometric sen-sitivities depend on the precision of the measurementof the oscillating frequency; and pulsed/continuous-wavecavity ring-down polarimetry (CRDP) [8–15] and ellip-sometry (CRDE) [16–19], where polarization-dependentabsorption and refraction/reflection through/by an op-tical medium is extracted with high sensitivity throughthe precise measurement of the signal-decay time and itspolarization beat frequency.A distinction among the aforementioned examples canbe made according to their respective decay constantsand oscillating frequencies. In routine NMR, typical de-cay times are in the 10 − −
10 s range, while frequen-cies are in the 10 −
800 MHz range; especially, portableNMR instruments operate in the 10 −
30 MHz frequencyrange [20–23]. In FID optical magnetometry, typical de- ∗ [email protected] cay times are in the 10 − − − Hz range (see for instance,Refs. [2, 5]). In CRDP/CRDE demonstrations, however,decay times are typically in the 10 − − − s range, whilepolarization beat frequencies are in the 1 −
100 MHz range.For all these applications, significant data processingis typically required to determine the signal parameters,and, in general, experimental sensitivity is improved byaveraging over many measurement runs. As such, whendeveloping portable instruments one needs to appropri-ately select the instrument’s sampling and acquisitionrates, but also carefully consider the computational cost,i.e. the calculation time, to analyze each acquired signal.For applications where the relevant time-scales are rela-tively long (10 ms - 1 s), such as NMR or FID magnetom-etry, there are several options that can provide preciseresults sufficiently fast (with respect to “single-events”),such as, e.g., frequency counters (see Refs. [24, 25] andreferences therein). However, in applications where therelevant time-scales are much shorter than a few ms, as inthe case of CRDP/CRDE, acquisition and computationalspeeds ultimately define, respectively, the measurementand analysis repetition rates.It is instructive to consider here a general data acquisi-tion system, as depicted in Fig. 1: such a system collectsthe physical signal, in our case a damped sinusoidal sig-nal, through an analog front-end that typically performssome signal conditioning (e.g. analog filtering, signal am-plification). The acquired signal is subsequently digitizedusing an analog-to-digital converter and transferred to adata processing unit for parameter analysis. A partic-ular case example is the field-programmable gate array(FPGA), whose development has lead to the emergenceof stand alone multi-channel (e.g., four) high-precision(e.g., 16-bit) data acquisition systems with high sam-pling (e.g., 2 GS/s) and triggering rates (e.g., 1 MHz), a r X i v : . [ phy s i c s . i n s - d e t ] O c t FIG. 1.
Data acquisition block diagram - a damped sinusoidal, physical, signal acquired by an experimental apparatus isdigitized using an analog-to-digital (A/D) converter and then transferred to a data processing unit, which estimates the signalparameters and displays them on screen for real-time monitoring. which have nowadays become commercially available atcost-effective rates. In such systems, data are (typically)transferred via USB or PCIe interfaces and, as such, datatransfer rates as high as 5 GB/s are feasible (e.g. Tele-dyne SP Devices, ADQ series [26]). Using such systems,therefore, it is possible to perform (sub-)microsecond-resolved CRDP/CRDE measurements [19]. However, theprincipal limiting factors towards an online, real-time,processing system that needs to be operable in a running fashion with no dead-times, is the data-processing mod-ule of the overall data acquisition system (Fig. 1) andits limited memory storage capacity. Considering thatin CRDP/CRDE experiments demonstrated decay timeconstants are in the 10 − − − s range, it is importantto identify appropriate computational approaches thatcan be implemented in data acquisition systems, such asFPGA-based digitizers, to allow for online, real-time, sig-nal analysis at such fast time scales.Time- and frequency-based computational methodolo-gies for rapid parameter estimation have been developedwithin the context of CRD spectroscopy, and these havebeen evaluated and compared in terms of their speed andprecision [27–29]. Most notably, Fourier transform meth-ods have been implemented on FPGAs for fast analysisof exponentially decaying signals [30], with demonstratedanalysis rates as high as 4.4 kHz [31]. However, while sev-eral works discuss the performance of various time- andfrequency-domain analysis algorithms for damped sinu-soidal signals [32, 33], a direct comparison between theirattainable precision and computational speed is currentlymissing.In this work we compare three specific analysis meth-ods of discretely sampled damped sinusoidal signals interms of their speed, and attainable accuracy and pre-cision. These methods are: (a) a time-domain least-squares analysis based on a Levenberg-Marquardt algo-rithm [34]; (b) a frequency-domain analysis based on afast Fourier algorithm [28, 29, 35, 36] in combination witha quadratic interpolation of the frequency componentsof the resulting Fourier transform [17]; and (c) a time-domain analysis based on the Prony method [37]. Weevaluate their efficacy in terms of the signal’s parame-ters, and discuss how each of these affect the sensitivitylimits for each computational methodology. Finally, wepresent an experimental, proof-of-principle, demonstra-tion of the capabilities of such methods for the online analysis of CRDP signals. II. TheoryA. Damped sinusoidal signals
A damped sinusoidal signal can be characterized in termsof a model function as: y ( t ) = A · e − t/τ · cos (2 π · f · t + φ ) + y , (1)where t is the (discretely sampled) independent (time)variable of the signal, A is the amplitude of oscillation, τ is the characteristic decay time, f and φ the frequency andphase of the oscillation, respectively, and y is a globalsignal offset. Under realistic experimental conditions, allthe signal parameters will be time-dependent, and thepower spectral density of the signal will be proportionalto their respective noise contributions. Here, for sim-plicity, and to clarify the main results of our findings, weassume that τ and f are constant parameters and restrictthe investigation of noise contributions to the global off-set parameter, i.e. y ( t ), which we assume to be normallydistributed [ (cid:104) y ( t ) (cid:105) = 0, (cid:10) y ( t ) (cid:11) = σ y ].In Fig 2 we show an example of a discretely sampleddamped sinusoidal signal. For the analysis of such a sig-nal we consider four key parameters that affect the ex-pected precision and accuracy: a) the number of signaloscillations per typical decay time, f × τ ; b) the numberof samples per typical decay time, n × τ − (i.e. the sam-pling rate); c) the number of decay times measured ina measurement time window T m , T m × τ − ; and d) thesignal-to-noise ratio, defined as SNR = A × σ − y . B. Cram´er-Rao Lower Bound
The fundamental limit for the statistical uncertainty ofdetermining the oscillating frequency of a damped sinu-soidal signal (Eq. 1) is described by the Cram´er-Rao lowerbound (CRLB) [3, 38], which sets the lower limit on thevariance σ of any frequency estimator. The CRLB con-dition for the frequency extracted from a discrete dampedsinusoid is given by [6, 39–41], σ = 6(2 π ) SNR f BW T χ ( τ /T m ) , (2)where SNR is the signal-to-noise ratio of the signal; f BW is the sampling-rate-limited bandwidth of the measure-ment; T m is the measurement time window; and χ ( τ /T m ) FIG. 2. Example of a discretely sampled damped sinusoidas described by Eq. 1, for f × τ = 5, n × τ − = 200 /τ , andSNR = 2 . is a correction factor that takes into account the signaldecay, which is given by χ ( r ) = e /r − r cosh (2 /r ) − r ( r + 2) . (3)The factor χ ( τ /T m ) serves as a compensation factor inEq. 2 that penalizes measurement of the tails of the ex-ponential decay when the signal has effectively died out.Equation 2 remains valid under the condition that theperiod of the oscillation is much shorter than the decaytime of the signal and that a sufficient number of oscilla-tions occurs in it. Moreover, Eq. 2 dictates that any noisesources affecting the signal detection are contributing tothe fundamental CRLB limit through their effect on theSNR of signal.In Ref. [15], the authors demonstrate that the CRLBlimit is the appropriate estimator of the fundamentalsensitivity of frequency-based measurements within thecontext of CRDP, as the frequency measurements aredirectly translated into polarimetric results. However,one needs to carefully investigate whether different signalprocessing techniques can approach the CRLB, and if yes,under what conditions this is possible. Moreover, con-sidering our motivation is the development of a portableCRDP instrumentation operating with similar principlesas recent demonstrations of it [15, 19], we focus on in-vestigating and comparing different signal processing ap-proaches in terms of their speed and attainable accuracyand precision for damped sinusoidal signals with decaytimes in the range of 1 − µ s and frequencies in therange of 1 −
10 MHz.
C. Signal Analysis
1. Least-Squares Estimation of Nonlinear Parameters
For time-domain analysis we focus on an optimized leastsquares curve fitting approach based on the Levenberg-Marquardt algorithm (LMA) [34]. The algorithm mini- mizes the sum of the squared residuals, S = n (cid:88) i =1 [ y i − f ( t i , β )] , (4)where y i is the i th sample of the discretized recorded sig-nal y ( t ) (Eq. 1), t i is the i th time sample, and f ( x, β ) isthe non-linear fit function (Eq. 1) with β representing theguess fitting parameters for { A, τ, f , φ, y } . The LMA al-gorithm iteratively finds the optimal guess parameters β describing the recorded signal y .The LMA is, in itself, an efficient algorithm, but it re-lies heavily on the initial guess parameters of the iterativeprocess. However, we wish to identify the precision andspeed limitations of computational implementations of aleast-squares algorithm and, hence, we assume for ourcomputational investigations that the initial conditionsare well-defined and known in advance (with our exper-imental investigation we examine the dependence of theLMA algorithm on the initial guess parameters under re-alistic conditions; see Sec. IV C). Furthermore, the timerequired for the convergence of a fit using LMA is highlydependent on the platform used. In this work we chooseto work with a CPU-based code for the implementationof the LMA that employs a Python optimized package(SciPy) based on the MINPACK library [42].
2. Fast Fourier transform
For frequency-domain analysis we use a fast Fouriertransform (FFT) algorithm, as introduced by Cooley andTukey [35], to calculate the discrete Fourier transform(DFT) of the signal, F ( k ) = N − (cid:88) n =0 y n e − πiN n k , k = 0 , . . . , N − , (5)where y n is the n th sample of the discretized time-domainsignal y ( t ) (Eq. 1). We note here that the simplest andmost common implementations of the FFT algorithm in-troduced by Cooley and Tuckey assume that N is a powerof two.The Fourier transform of a monochromatic dampedsinusoidal signal corresponds to a single Fourier (fre-quency) component with a spectral width inversely pro-portional to the signal’s decay time. Our aim is to es-timate accurately and precisely the central value of thiscomponent, rapidly. One approach is to perform a least-squares curve fitting on the resulting FFT spectrum toobtain the central value of the frequency component andits width. However, the accuracy and precision of suchprocess depends strongly on the curve fit-model selectedand its initial guess parameters, but, importantly, thespeed of such an approach would be at least equal tothe overall time required to perform both the FFT andthe least-squares fitting. Furthermore, in order for suchan approach to be as precise as the direct time-domainanalysis approach using, e.g., LMA, one typically em-ploys additional data manipulation techniques (e.g. zero-padding, apodization).Here, we focus on algebraic approaches for the rapidextraction of the central value of the Fourier (frequency)component from the FFT spectrum. One such approachis to determine the center value of this component byconsidering the three closest neighbouring points to themaximum frequency value (peak): ( k i , F ( k i )) ≡ ( k i , b i ),( i = 1 , , max as:f max = k b ( b − b ) + k b ( b − b ) + k b ( b − b )2 [ k b ( b − b ) + k b ( b − b ) + k b ( b − b )] . (6)It is important to emphasize that the selection of theneighbouring points is crucial for the accuracy (not thespeed) of the frequency estimation using such an ap-proach. For high sampling rates, for instance, one canchoose - symmetrically, or even asymmetrically - pointsfurther away from the closest neighbouring points tothe peak, and preferably points lying near to the half-maximum of the Fourier component [this can be easilypre-set in the algorithm if the decay time and the sam-pling rate are (approximately) known in advance]. Suchan algebraic approach on analysing FFT spectra has al-ready been successfully implemented for rapid frequencyestimation in CRDP-based experiments (see Ref. [17]).By choosing such an approach, we ensure that the com-putational speed remains as close as possible to the speedrequired to employ a FFT algorithm.There exist several CPU-based codes available for FFTanalysis, but for an appropriate speed comparison be-tween the alternative signal processing methodologiespresented in this work, we use a DFT algorithm directlyfrom a Python-based scientific environment (NumPy; wenote here that we do not observe in our analysis any dif-ferences between different FFT libraries in Python suchas SciPy and NumPy).
3. Prony
The Prony method, is a time domain approach origi-nally designed for processing discrete time signals thatare superpositions of damped sinusoids. The Pronymethod is closely related to the Matrix Pencil method(both estimate the signal as a sum of complex expo-nentials) [43, 44], the latter being used in NMR analy-sis [45, 46]. However, Prony analysis takes a polynomialapproach in parameter (frequencies and damping factors)estimation whereas Matrix Pencil Method locates the sig-nal parameters by finding the eigenvalues to a matrixpencil.The application of the Prony method follows in threesteps: (a) an autoregressive model is built employing dis-crete measurements; (b) the roots of the characteristicpolynomial for the corresponding finite difference equa-tion are statistically estimated; and (c) estimates of theparameters of the signal are derived from the roots.For the special case of one samped sinusoid, a discretetime sampling of such a signal gives rise to an autore-gressive model of order 3 where the measurement y k at time k is expressed through 3 preceding measurements ina linear way: y k +3 + α y k +2 + α y k +1 + α y k = 0 , (7)where k varies from 0 to n + 2, and n + 5 is the numberof measurements (sampling points). The coefficients α i are determined by any of the linear systems y k +2 y k +1 y k y k +3 y k +2 y k +1 y k +4 y k +3 y k +2 · α α α = − y k +3 y k +4 y k +5 (8)with k = 0 , . . . , n .In the presence of noise, the 3 × α i can be found, e.g., by theleast square method minimising the loss function n (cid:88) k =0 ( y k +3 + α y k +2 + α y k +1 + α y k ) . These constitute the characteristic polynomial equation q ( z ) = z + α z + α z + α , (9)whose roots u , u ± incorporate the parameters of thesignal. In the special case under study, u = e − ∆ τ , u ± = e − ∆ τ e ± i π f∆ , where ∆ is the sampling time inter-val. Then the frequency and decay constant of the signalare found as τ = − ln( u ) / ∆ , (10)f = Im[ln( u + ) / ∆] . (11)The roots can be calculated, e.g., by the Cardano for-mulas, or by employing α = − u (in the limit of nonoise). The root u + can be distinguished from u − as theone with positive imaginary part if u + + u − >
0, andnegative otherwise.A practical realization of this scheme has to take intoaccount the role of the sampling rate n × τ − = . Evenin the absence of noise, there exist singular values at n × τ − = f πN (with positive integer N ) for which the matrixin Eq. 8 becomes degenerate (degeneracy occurs for half-integer N ), and the sampling rate has to be chosen tobe different from such singular values. Another thing isthat the coefficients α i depend on f via cos(f∆). Keepingthe sampling rate above f max ∆ confines f∆ in [0 , π ] anddetermines f from Eq. 11 uniquely.Furthermore, in the presence of noise, the accuracystill depends on n × τ − even if it exceeds f π . So, ifthe sampling rate is too high, the sampled points aretoo close to each other (note that their number is fixed),and small variations of the signal (a smooth function)from point to point are distorted by random jumps whichdeteriorate estimation. Therefore the frequency shouldbe bounded from below, say with f min . In practice, fora reasonable SNR the dependence of the result on thesampling rate is weak in a wide range of n × τ − values,and this observation can be used for estimation. III. MethodsA. Signal Simulation
To compare the three methods of analysis on their re-spective precision and accuracy in estimating the centralfrequency of damped sinusoidal signals (Fig. 2), we gen-erate and analyze sets of 500 such signals on a homemadePython CPU-code on a Windows 10 workstation [CPU:AMD Ryzen 7 2700, RAM: 16.0 GB 1330 MHz DDR4].All simulated signals have the following non-changing pa-rameter values: A = 1, (cid:104) y (cid:105) = 0, and φ = 0. We alsochoose the following baseline values for the key parame-ters of each simulated signal: f × τ = 5 , n × τ − = 1000 /τ , T m × τ − = 5, and SNR = 2 . We choose here a highbaseline value for the SNR to clearly examine whetherthe computational approaches can reach the fundamen-tal CRLB limit as a function of the other key signal pa-rameters. We proceed by varying each key parameterover several orders of magnitude while keeping the otherparameters at their baseline value, to explore the depen-dence of the precision and accuracy of each computa-tional approach on these parameters. B. Precision and accuracy
As a way to quantify the precision of each computationalapproach we use the standard deviation from the distri-bution of frequency values obtained through the analy-sis of the 500 simulated signals, i.e. σ f , to estimate thefractional uncertainty σ f / f (i.e. smaller fractional un-certainty corresponds to higher precision). Similarly, wedefine as the accuracy of a method asaccuracy = 1 N N (cid:88) i =1 | f i, est − f act | f act , (12)where f i, est is the frequency estimated by the analysismethod for a single signal, f act is the actual (input) fre-quency of the simulated signal (again here, N = 500).An analysis method is predicted to have no bias as longas the accuracy of its frequency estimation falls withinthe precision of the estimation. C. Computation time
We determine the speed of each computation methodby estimating the time required to analyze a single signalusing the internal timing functions of the Python soft-ware (e.g., function timeit ). IV. ResultsA. Precision and accuracy
In Fig. 3 we present results on the precision and accu-racy achieved when analyzing discrete simulated dampedsinusoidal signals using the least-squares, FFT, andProny computation analysis methods, as a function ofvarying signal conditions. In particular:
Frequency -
In Fig. 3 (a) we show a comparisonof the attainable precision and accuracy between thethree different approaches as a function of the frequency of oscillation. For all frequencies the least-squares ap-proach results in optimal accuracies compared to theother approaches, with the FFT approach being consis-tently less accurate (this is largely related to the peak-finding algebraic methodology we employ here). TheProny method becomes particularly inaccurate for lowfrequencies, which is related with the computational for-mulation of the Prony method that doesn’t allows us toinvestigate a large parameter space without approach-ing singular points in the analysis. In terms of preci-sion, both the least squares and FFT methods approachclosely the CRLB, the latter being approximately a fac-tor of three less precise than the former, while for bothmethodologies the precision is not influenced by the fre-quency value (the least squares method deviates from theCRLB at frequencies f × τ < .
2, as expected, since theobserved time window does not contain a full period ofoscillation). The Prony method yields results with pooraccuracy for f × τ < Sampling Rate -
In Fig. 3 (b) we show a similarcomparison as a function of the sampling rate, i.e. as afunction of the number of sample points per decay time n × τ − . We note again here that for these estimationswe choose a constant frequency of f × τ − = 5. As such,the Nyquist criterion limits the lowest sampling rate for asensible frequency estimate to n × τ − = 20. In terms ofaccuracy and precision the least squares method yieldsoptimal results, while the FFT method yields optimalprecision but relatively poor accuracy (at the 10 − level),both related to the peak-finding algorithm (these can beimproved by performing additional signal manipulation,such as zero padding, but this will significantly affect thecomputational speed). The Prony method provides accu-rate and precise results for low sampling rates, but thesedeteriorate for high sampling rates (see Sec. II C 3). Wealso observe that the least squares method is limited bythe CRLB over the entire range of sampling range weinvestigate, while the FFT method remains consistentlyless precise. Measurement window -
In Fig. 3 (c) we presentresults as a function of the measurement window, i.e. T m × τ − . The least-squares method yields optimalaccuracy ( ∼ − ) and precision ( ∼ − ) results for T m × τ − >
2, however, for short measurement windowsthe precision becomes poor [as predicted by the CRLBlimit, Eq. 2]. The FFT method reaches its optimum ac-curacy ( ∼ − ) and precision for T m × τ − ≈
4. Impor-tantly, we observe that both the least-squares and FFTmethods reach the CLRB limit for T m × τ − >
4, with anoptimal measurement window for precise signal analysisusing both methods to be ∼ τ . Similar conclusions havealready been reported in Ref. [29], suggesting that ∼ τ can be considered the optimum repetition rate for, e.g.,CRDP/CRDE experiments, as compared to the longeracquisition windows ( T m × τ − >
5) typically required intraditional CRD spectroscopy [47]. The Prony methodreaches similar accuracies as the FFT method but its
FIG. 3. Results of attainable fractional uncertainty, σ f / f, and accuracy using the FFT (black points), Prony (red points),and least-squares (blue points) analysis methods as a function of varying signal conditions: (a) signal frequency, f × τ ; (b)sampling rate, n × τ − ; (c) measurement time-window, T m × τ − ; and (d) signal-to-noise-ratio (SNR). In combination withthe performance of each method, the fundamental frequency estimation limit given by the Cram´er Rao lower bound (CRLB)is also shown (dashed gray line). The least-squares method yields results close to the CRLB limit in virtually all situations,with the FFT method yielding similar results, while the Prony yields poor results under most selected conditions. precision is two orders of magnitude larger than the pre-dicted CRLB limit. SNR -
The final key parameter we vary is the signal’sSNR, with the results seen in Fig. 3 (d). The least-squaresanalysis yields results close to the CRLB limit, while theprecision attained using FFT analysis method is approx-imately a factor of two ( ×
2) higher. Notwithstanding,we see that for an optimum measurement time-windowof 5 τ [Fig. 3 (c)] and a SNR ≈ both the least-squaresand FFT methods yield precisions at the ∼ − levels.In contrast, the Prony method does not provide reliablefrequency estimates for signals with SNR < B. Speed
In Fig. 4 we present results on the dependence of thefitting (computation) time for each method on the num-ber of data points in a single damped sinusoidal signal,which we also compare with the attainable precision foreach case. For these simulations, we use the results pre-sented in Fig. 3 to choose optimum values for the signal’skey parameters: f × τ = 5, SNR = 2 and T m × τ − = 5(such values are also realistically attainable in experi-ments; see discussions in Ref. [15]).Overall we observe a non-linear increase in the calcula-tion time as the number of samples is increased, with theleast-squares and Prony algorithms being more than anorder of magnitude slower than the FFT+peak-findingalgorithms. In addition, while the least-squares method FIG. 4. Dependence of precision and calculation time (greenpoints) on the number of samples for the: a) FFT (blackpoints), b) Prony (red points), and c) least-squares (bluepoints) evaluation methods. In all cases, we compare theobtained results with the expected Cram´er-Rao lower boundlimit (CRLB; solid gray line) for increasing sample size [ n ( × τ = 5, SNR = 2 and T m × τ − =5. Using the FFT algorithm methodology, for a signal with10 data points we obtain fractional uncertainties at the sub-10 − levels within a computational time of ∼ µ s, while forthe same conditions using a least-squares approach we obtainsimilar precisions but for a computational time of ∼ > ms computational times. reaches the CLRB limit, the FFT method yields frac-tional uncertainties approximately a factor of two largerthan the predicted CRLB limit, with the Prony methodpractically never reaching optimal precision levels. Mostimportantly, we observe that for a discrete signal with ∼ sample points, using the FFT+peak-finding algo-rithm one can achieve ppm sensitivities (10 − fractionaluncertainties) for computational times of ∼ µ s. Un-der the same conditions, the least-squares algorithm re-quires approximately ∼ − fractional uncertainties) requiring long ( > ms)computational times. C. Experiment
As a proof-of-principle demonstration for the capabil-ities of our analysis methodology we use a CRD-basedpolarimetric instrument we have developed in our lab-oratory for the attainment of experimental CRDP, i.e.damped sinusoidal, signals (see for details Ref. [15]).Briefly, the ring-down cavity of the instrument has atotal length of 0.60 m and consists of two concave mir-rors with radii of curvature of 1 m and specified reflec-tivity R ∼ λ = 408 nm) that we rapidly pulse to initiate ring-down events [with the use of an acousto-optic modula-tor (AOM; Gooch and Housego 3200-125)]. In our op-tical setup we can generate CRDP signals with ring-down times in the 0.3-1.5 µ s range (depending on theusage of intracavity optics), at repetition rates as highas 100 kHz, that we record and digitize using a 14-bitdigitizer (Teledyne, ADQ14DC-2X-PCIE, dual channelDC-coupled operation; sample rates of 2 GS/s per chan-nel), which has a maximum acquisition rate of 100 kHz(mainly limited by the data transfer rate), 14-bit resolu-tion per channel, and permits on-board channel subtrac-tion and signal averaging.For our demonstration, we use the (non-resonant)Faraday effect of a 6.35(1) mm thick, AR-coated SiO substrate (FiveNines Optics; AR coated by FiveNineOptics with specified R < . θ F [11, 15, 48] that result in CRDP signals with polariza-tion beat frequencies in the range of 1-3 MHz [the beatfrequency is proportional to the induced Faraday rotationas: f = θ F · FSR /π , where FSR= ( c/L rt ) is the cavity’sfree spectral range, with c the speed of light and L rt theround-trip cavity length].To demonstrate our ability to monitor and analyzeCRDP experimental signals in a running fashion at ∼ kHzrates using the FFT+peak-finding analysis algorithm [inaccord with the results shown in Fig. 4 (a)], we proceed asfollows: we initiate ring-down events at a rate of 100 kHzand continuously record the polarimetric ring-down sig-nals (i.e., the photo-detector signals) for ≈
220 ms, while
FIG. 5.
Rapid analysis of experimental signals: a) CRDP experimental signals showing polarization beat frequencies generatedvia the Faraday effect on a SiO substrate, as recorded by two orthogonal channels of a linear balanced polarimeter (rightside). By subtracting these and averaging 40 consecutive traces, at a repetition rate of 100 kHz, we acquire a CRDP trace(i.e. damped sinusoidal oscillation) with a SNR=40 within 400 µ s. (b) Direct observation of frequency shifts with respect tothe (bias) Faraday polarization beat frequency [f = 1 . ∼ µ s of calculation time [(c)]. The (black) line is the result of a nonlinear least-squares regression analysisused to fit a sigmoid function to the data, demonstrating that we can resolve sub- µ s dynamics in a running-fashion. Note thatthe CRLB estimation for such signals (with 10 data points, SNR=40, and τ ≈ . µ s) is ∼ . ∼ µ s per trace and remains constant during the analysis of the complete set of traces, while a least-squaresapproach requires computational times ranging from ∼ − channel subtraction allows for the generation of dampedsinusoidal signals (CRDP traces) and on-board signal av-eraging enables the average of 40 consecutive traces; eachtrace, therefore, requires an integration time of 400 µ s.In Fig. 5 (a) we show such an experimentally acquiredCRDP trace with a Faraday-rotation-related polariza-tion beat frequency of f = 1 . τ = 0 . µ s, and an SNR (cid:39)
40. Note that, giventhe digitizer’s sampling rate, the CRLB limit in estimat-ing the signal’s beat frequency is 1.3 kHz, i.e. σ f / f = 10 − (Eq. 2; Fig. 3). Furthermore, using an additional, home-made, solenoid [with a length of 2 . . substrate,we induce a rapid frequency shift on the recorded CRDPsignal by applying a (rapidly pulsed) external magneticfield [using a USB controlled metal-oxide-semiconductorfield-effect transistor-based switching circuit resulting inswitch-on times of < µ s].In Fig. 5 (b) & (c) we show our experimental results.We see that using an online FFT+peak-finding approach we can analyze CRDP traces at a constant rate of ∼ µ s time scales [evident from the analysis of our resultsusing a sigmoid function that yields a ≈ µ s rise-time;Fig. 5 (b)]. We emphasize here that the frequency fluc-tuations present in the recorded signals [Fig. 5 (b)] arethe result of experimental noise sources (see Ref. [15]).As a comparison, in Fig. 5 (c) we also demonstrate thata least-squares approach would require computationaltimes ranging from ∼ V. Discussion & Conclusion
In this work, we consider different time- and frequency-domain-based computational algorithms for the rapidestimation of the signal parameters of damped sinu-soidal signals. We analyze their accuracy and precisionin terms of key signal parameters and estimate thecomputational time required to obtain these usingstandard computational platforms (e.g. a desktopcomputer) and software (e.g. Python). Overall, wesee that a time-domain-based least-squares algorithmreaches the expected fundamental estimator limits interms of precision and accuracy, and requires ms-longcomputational times to obtain fractional uncertaintiesat the ppm levels (for signals with high SNR), while aFourier-based algorithm can achieve similar sensitivitiesat, at least, an order of magnitude faster computationaltimes, even when one employees standard computationalplatforms and algorithms. We also consider in compar-ison to the least-squares and FFT analysis methods analternative computational approach based on the Pronymethod, recently proposed for rapid analysis of dampedsinusoidal signals; we observe that an implementation ofthe Prony method for our parameter-range of interest,fails to provide comparable accuracies and precisions tothe least-squares and FFT methods, particularly withinsimilar computational time-scales. We then validateour results using an experimental CRDP setup and aFRGA-based acquisition system, to demonstrate theonline recording and analysis of damped sinusoidalsignals in a running fashion at ∼ kHz rates using anFFT+peak-finding alogrithm, and, in particular, wedemonstrate the ability to observe signal changes attime scales as fast as 10 µ s.Overall, based on our results, we recommend adoptingFFT analysis approaches for rapid parameter analysisof damped sinusoidal signals within the context ofCRDP/CRDE (and similar) techniques, which offers anoptimum combination of speed and performance.As a concluding remark we note that the exact com-putational speeds for the presented analysis methodsdepend significantly on the computational platform andsoftware used. In Ref. [37], the authors demonstratemore than an order of magnitude improvements incomputational times for the FFT and Prony methodsby implementing alternative computational pack- ages/softwares. As such, we expect that an optimizedFFT algorithm may very well outperform our results bymore than an order of magnitude, suggesting, therefore,that even with the use of standard data-acquisitionsystems, highly precise microsecond-resolved signalparameter estimation is possible. This becomes par-ticularly important for our application of interest,CRDP: using the results presented in Ref. [15] we canestimate that in CRDP experiments with typical decaytimes constants of a few µ s and oscillating frequenciesof > > , we see thatit is possible to perform online, in a running-fashion,microsecond-resolved experiments with µ rad polarimet-ric sensitivities (i.e. sub-nrad/ √ Hz sensitivities). Sucha possibility is of paramount importance for emergingapplications in chiral sensing and analysis, in surfacecatalysis, and indicates that real-time monitoring ofgas/liquid flows in GC/HPLC and of surface (chiral)dynamics is nowadays feasible. Finally, we anticipatethat an overall improvement in computational speedscould also result from implementing a network theoryapproach for rapid (spectroscopic/spectropolarimetric)signal analysis; for intance, network theory has beenrecently used for the classification and search of spectralfeatures of various molecules [49, 50]. However, rapidparameter extraction of time-sampled signals couldprove to be nontrivial with a network approach, and wewill investigate such a possibility in future works.
Acknowledgments
This work was supported by the European Com-mission Horizon 2020, project ULTRACHIRAL (GrantNo. FETOPEN-737071), and by the European Union’sSeventh Framework Programme for Research, Techno-logical Development, and Demonstration, under theproject ERA.Net RUS Plus, grant agreement no. 189(EPOCHSE). EW, TC, and AM were supported in partby Innovate UK KTP 010819. LB and JV are grateful toDmitry Budker for his constant help, support and com-ments. LB is grateful to Michael Everest for his help andsupport, and to Amelia Meath for fruitful discussions. [1] H. G¨unther,
NMR spectroscopy: basic principles, con-cepts and applications in chemistry (John Wiley & Sons,2013).[2] I. M. Savukov and M. V. Romalis, Phys. Rev. Lett. ,123001 (2005).[3] C. Gemmel, W. Heil, S. Karpuk, K. Lenz, C. Lud-wig, Y. Sobolev, K. Tullney, M. Burghoff, W. Kilian,S. Knappe-Gr¨uneberg, et al. , The European PhysicalJournal D , 303 (2010). [4] A. Nikiel, P. Bl¨umler, W. Heil, M. Hehn, S. Karpuk,A. Maul, E. Otten, L. M. Schreiber, and M. Terekhov,The European Physical Journal D , 330 (2014).[5] Z. D. Gruji´c, P. A. Koss, G. Bison, and A. Weis, TheEuropean Physical Journal D , 135 (2015).[6] D. Hunter, S. Piccolomo, J. D. Pritchard, N. L. Brockie,T. E. Dyer, and E. Riis, Phys. Rev. Applied , 014002(2018).[7] D. Hunter, R. Jim´enez-Mart´ınez, J. Herbsommer, S. Ra-maswamy, W. Li, and E. Riis, Opt. Express , 30523 (2018).[8] T. M¨uller, K. B. Wiberg, and P. H. Vaccaro, The Journalof Physical Chemistry A , 5959 (2000).[9] T. M¨uller, K. B. Wiberg, P. H. Vaccaro, J. R. Cheeseman,and M. J. Frisch, J. Opt. Soc. Am. B , 125 (2002).[10] D. Sofikitis, L. Bougas, G. E. Katsoprinakis, A. K. Spili-otis, B. Loppinet, and T. P. Rakitzis, Nature , 76(2014).[11] L. Bougas, D. Sofikitis, G. E. Katsoprinakis, A. K. Spili-otis, P. Tzallas, B. Loppinet, and T. P. Rakitzis, TheJournal of chemical physics , 09B603 1 (2015).[12] P. Dupr´e, Phys. Rev. A , 053817 (2015).[13] A. Spiliotis, M. Xygkis, E. Klironomou, E. Kardamaki,G. Boulogiannis, G. Katsoprinakis, D. Sofikitis, andT. Rakitzis, Chemical Physics Letters , 137345(2020).[14] A. K. Spiliotis, M. Xygkis, E. Klironomou, E. Kar-damaki, G. K. Boulogiannis, G. E. Katsoprinakis,D. Sofikitis, and T. P. Rakitzis, Laser Physics , 075602(2020).[15] J. C. Visschers, O. Tretiak, D. Budker, and L. Bougas,The Journal of Chemical Physics , 164202 (2020).[16] V. Papadakis, M. A. Everest, K. Stamataki, S. Tzortza-kis, B. Loppinet, and T. P. Rakitzis, in Instrumentation,Metrology, and Standards for Nanomanufacturing, Op-tics, and Semiconductors V , Vol. 8105, edited by M. T.Postek, International Society for Optics and Photonics(SPIE, 2011) pp. 104 – 112.[17] K. Stamataki, V. Papadakis, M. A. Everest, S. Tzortza-kis, B. Loppinet, and T. P. Rakitzis, Applied optics ,1086 (2013).[18] D. Sofikitis, K. Stamataki, M. A. Everest, V. Papadakis,J.-L. Stehle, B. Loppinet, and T. P. Rakitzis, Opt. Lett. , 1224 (2013).[19] D. Sofikitis, A. K. Spiliotis, K. Stamataki, G. E. Katso-prinakis, L. Bougas, P. C. Samartzis, B. Loppinet, T. P.Rakitzis, M. Surligas, and S. Papadakis, Appl. Opt. ,5861 (2015).[20] H. Lee, E. Sun, D. Ham, and R. Weissleder, NatureMedicine , 869 (2008).[21] J. Perlo, V. Demas, F. Casanova, C. A. Meriles,J. Reimer, A. Pines, and B. Bl¨umich, Science , 1279(2005).[22] K. Lei, H. Heidari, P. Mak, M. Law, F. Maloberti, andR. P. Martins, IEEE Journal of Solid-State Circuits ,284 (2017).[23] K.-M. Lei, D. Ha, Y.-Q. Song, R. M. Westervelt, R. Mar-tins, P.-I. Mak, and D. Ham, Analytical Chemistry ,2112 (2020).[24] R. Prigl, U. Haeberlen, K. Jungmann, G. zu Putlitz,and P. von Walter, Nuclear Instruments and Methods inPhysics Research Section A: Accelerators, Spectrometers,Detectors and Associated Equipment , 118 (1996).[25] H. Dong, H. Liu, J. Ge, Z. Yuan, and Z. Zhao, IEEETransactions on Instrumentation and Measurement , 2187 (2004).[28] M. Mazurenka, R. Wada, A. J. L. Shillings, T. J. A.Butler, J. M. Beames, and A. J. Orr-Ewing, Applied Physics B , 135 (2005).[29] M. A. Everest and D. B. Atkinson, Review of ScientificInstruments , 023108 (2008).[30] T. G. Spence, M. E. Calzada, H. M. Gardner, E. Leefe,H. B. Fontenot, L. Gilevicius, R. W. Hartsock, T. K.Boyson, and C. C. Harb, Opt. Express , 8804 (2012).[31] G. Bostrom, D. Atkinson, and A. Rice, Review of Sci-entific Instruments , 043106 (2015).[32] E. Aboutanios, IEEE Instrumentation & MeasurementMagazine , 8 (2011).[33] E. Aboutanios, IEEE Transactions on Signal Processing , 501 (2009).[34] J. J. Mor´e, in Numerical analysis (Springer, 1978) pp.105–116.[35] J. W. Cooley and J. W. Tukey, Mathematics of compu-tation , 297 (1965).[36] T. K. Boyson, T. G. Spence, M. E. Calzada, and C. C.Harb, Opt. Express , 8092 (2011).[37] E. Wilson, T. M. Conneely, A. Mudrov, and I. Tyukin,IFAC-PapersOnLine , 269 (2019).[38] Y.-X. Yao and S. M. Pandit, IEEE Transactions on signalprocessing , 878 (1995).[39] Ying-Xian Yao and S. M. Pandit, IEEE Transactions onSignal Processing , 878 (1995).[40] C. Gemmel, W. Heil, S. Karpuk, K. Lenz, C. Lud-wig, Y. Sobolev, K. Tullney, M. Burghoff, W. Kil-ian, S. Knappe-Gr¨uneberg, W. M¨uller, A. Schnabel,F. Seifert, L. Trahms, and S. Baeßler, The EuropeanPhysical Journal D , 303 (2010).[41] H.-C. Koch, G. Bison, Z. D. Gruji´c, W. Heil,M. Kasprzak, P. Knowles, A. Kraft, A. Pazgalev,A. Schnabel, J. Voigt, and A. Weis, The European Phys-ical Journal D , 202 (2015).[42] J. J. More, B. S. Garbow, and K. E. Hillstrom, (1980),10.2172/6997568.[43] Y. Hua and T. K. Sarkar, IEEE Transactions on Acous-tics, Speech, and Signal Processing , 814 (1990).[44] F. Sarrazin, A. Sharaiha, P. Pouliguen, J. Chauveau,S. Collardey, and P. Potier, in (2011) pp. 1–4.[45] Y.-Y. Lin, P. Hodgkinson, M. Ernst, and A. Pines, Jour-nal of Magnetic Resonance , 30 (1997).[46] S. Fricke, J. Seymour, M. Battistel, D. Freedberg,C. Eads, and M. Augustine, Journal of Magnetic Reso-nance , 106704 (2020).[47] H. Huang and K. K. Lehmann, The Journal of PhysicalChemistry A , 13399 (2013).[48] L. Bougas, G. E. Katsoprinakis, W. von Klitzing,J. Sapirstein, and T. P. Rakitzis, Phys. Rev. Lett. ,210801 (2012).[49] D. P. Zaleski and K. Prozument, The Jour-nal of Chemical Physics , 104106 (2018),https://doi.org/10.1063/1.5037715.[50] R. T´obi´as, T. Furtenbacher, I. Simk´o, A. G. Cs´asz´ar,M. L. Diouf, F. M. J. Cozijn, J. M. A. Staa, E. J.Salumbides, and W. Ubachs, Nature Communications11