
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

39313497
70392
10.1038/s41598-024-70392-9
Article
An analytic, efficient and optimal readout algorithm for compact interferometers based on deep frequency modulation
Eckhardt Tobias tobias.eckhardt@physik.uni-hamburg.de

Gerberding Oliver
https://ror.org/00g30e956 grid.9026.d 0000 0001 2287 2617 Institut für Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany
23 9 2024
23 9 2024
2024
14 219888 8 2024
16 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Compact laser interferometers with large dynamic range are one of the core emerging tools to improve low frequency performance in gravitational wave detectors by providing local displacement sensing with sub 1 pmHz-0.5 precision. Strong sinusoidal frequency modulations are used in such laser interferometers to create heterodyne-like photodetector signals from which the phase and other parameters, such as the absolute distance, can be extracted. The nested sinusoidal function in such signals is a challenge for the real-time parameter estimation in low-noise applications. In this article, we present an algorithm to calculate exact signal parameters in a non-iterative way from such interferometric signals. The algorithm makes use of a recurrence relation between Bessel functions to enable a direct extraction of modulation parameters from the signal. Additionally, the algorithm is capable of dealing with high phase dynamics where the Doppler-shift of the signal becomes relevant and can limit the range and precision of the parameter estimation, if not accounted for. Simulations show that the algorithm is computationally efficient, can be well parallelised and the phase estimation is close to optimal precision given by the Cramer–Rao lower bound of the signal parameters.

Subject terms

Techniques and instrumentation
Mathematics and computing
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)EXC 2121 -- 390833306 German Federal Ministry of Education and Research05A20GU5 05A23GU5 Gerberding Oliver Universität Hamburg (1037)Open Access funding enabled and organized by Projekt DEAL.

issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Displacement sensing is a core metrological task that can be well addressed using laser interferometers. Compact interferometers that operate over more than a single fringe of the laser wavelength and measure in the mHz–Hz–kHz regime with sub 1 pmHz-0.5 precision are studied for fundamental and applied research. This prominently includes the displacement readout of ultra-precise inertial sensors for space-based geodesy missions and gravitational wave experiments1–4 and local inertial and displacement sensing in current and future ground-based gravitational wave detectors like LIGO, Virgo, KAGRA and the Einstein Telescope5–8. E.g. in LIGO, dozens of very compact and precise displacement sensors are required to isolate the mirrors from the ground-motion. More precise sensors that probe the 6 degrees of freedom can improve the mirror alignment, better suppress the mirror motion at mechanical suspension resonances and improve the active seismic isolation. Higher readout precision will generally result in lower mirror control noise, which is currently dominating detectors at low frequencies. Having compact sensors enables their deployment within the complex suspension systems9.

A class of interferometers studied for such applications use a scheme that we refer to here as Deep-Frequency Modulation Interferometry (DFMI)10. In DFMI a laser is strongly modulated in frequency and sent into several unequal arm-length interferometers, achieving both compact form and high precision displacement readout1,11. The interferometric signals in DFMI contain a nested sinusoidal function, which hides the desired signal parameters. To provide a linear estimate of these parameters the first realisations of DFMI have relied on an iterative fit algorithm that was developed by Heinzel et al. for a related technique called Deep-Phase Modulation Interferometry (DPMI)12. DPMI can be described as a modern variant of techniques like Phase Shifting Interferometry13. It uses a strong phase modulator in one arm of an interferometer (or multiple interferometers) to generate heterodyne like signals which are effectively the same as for DFMI. The iterative fit used for DFMI and DPMI, as well as other readout algorithms studied for DFMI14 and similar techniques13,15–17, use either approximations (which limit the precision) or computationally expensive, iterative fitting routines, giving them a disadvantage for a high precision test mass readout or the realisation of many channels and sensors. This paper presents an analytic algorithm to realise this parameter estimation in an efficient and optimal way for both DFMI and DPMI.

The algorithm we present aims to improve upon the readout by making use of a recurrence relation between the (non-linear) Bessel functions that appear when describing the signal in Fourier domain. It allows us to calculate all appearing parameters analytically without the need for approximations while achieving close to optimal precision (given by the Cramer–Rao lower bound of the estimates18). This is necessary to achieve typical precision levels around or below 100fm/Hz in displacement1, corresponding to a phase noise of ≈4·10-7rad/Hz for a laser wavelength of 1550nm. Additionally, we extend this algorithm to operate even with larger linear dynamics and we provide simulation results to verify and quantify the algorithms performance.Table 1 Nomenclature of the signal parameters.

A	Amplitude	P0	Mean laser power	
B	Offset	ω0	Mean laser frequency	
m	Modulation index	Δω	Modulation depth	
ψ	Modulation phase	fS	Sampling frequency	
φ	(Global) phase	fR	Readout frequency	
ωm	Modulation frequency	t	Time	
ΔL	Arm-length difference	τ	Propagation-time	
			(of the light)	

Fig. 1 Sketch of a DFMI setup with laser frequency ωDFM=ω0+Δωsin(ωmt+ψ) and its unequal arm-lengths. The setup measures the arm-length difference ΔL:=L2-L1 once encoded in the phase φ (referred to as microscopic distance)and once in the modulation index m macroscopic arm-length difference ΔL:=L2-L1.

The signal

Figure 1 shows the core topology of DFMI. The general form of the measured signal is given by1 s(t)=:B+Acosm·sinωmt+ψ+φ.

For the coefficients we use the nomenclature used by10 as written in Table 1.

The main parameter of interest is the phase φ, which encodes the microscopic displacement of the test mass, as well as other relevant information, like relative beam tilts if a quadrant photodiode is used for differential wavefront sensing19–21.

Using the Jacobi-Anger identity, we write the Fourier series of the signal (1) as2 s(t)=B+A∑n∈ZJn(m)cos(n(ωmt+ψ)+φ)

with Jn(m) as Bessel function of the first kind. Using the symmetry of the Bessel functions J-n(x)=(-1)nJn(x)forn∈Z we can also rewrite the series with a summation index n∈N as3 s(t)=B+AJ0(m)cosφ+2A∑n=1∞[J2n(m)·cosφ·cos(2n(ωmt+ψ))-J2n-1(m)·sinφ·sin((2n-1)(ωmt+ψ))]

Throughout this article we refer to the individual parts of the Fourier series (2) and (3) as ’signal harmonics’ or harmonic frequencies. Figure 2 shows a time-series and the power spectral density of an example signal with these harmonics clearly visible as ’peaks’ in the spectrum.Fig. 2 Time-series and power spectral density of an example ideal signal as defined in (1) with modulation index m=7. The frequency peaks are the individual parts of the Fourier series of the signal (3). The dotted line (here given by Jn(7)) is the enveloping function for the peaks. In this simulation, no additional noise besides the intrinsic digitization noise is present leading to a noise floor around 10-32. In a real experimental setup, this noise floor will be higher due to other additional noise sources which will ultimately limit how many harmonics can be resolved and used for the parameter estimation.

For DFMI and similar techniques, the modulation index scales with the differential arm-length ΔL. E.g. the modulation index is given by the product m:=Δω·δτ of the modulation depth Δω and propagation time-difference δτ (related to the arm-length difference via ΔL=c·δτ, with c as speed of light)10. Hence, m encodes a macroscopic length in DFMI that is complementary to the microscopic length in φ:=ω0·δτ (which is periodic and limited to φ∈[0,2π)).

From (3) we see that even and odd harmonics behave like different quadratures of the phase φ. These ’quadratures’ contain Bessel functions which are not periodic and scale differently with changing path-length differences ΔL. The Bessel functions are not invertible, meaning that there is no direct way to calculate the modulation index m (and the corresponding macroscopic path-length difference) for these measured harmonics. By use of a recurrence relation for the Bessel functions one can still calculate m from multiple harmonics to get the macroscopic distance and subsequently divide them out to reveal simpler trigonometric relations.

Most studied readout algorithms10,12,13,15–17 make use of the Fourier-series description ((2) or (3)), where they employ either a fast Fourier transform or a dedicated demodulation22 to determine the amplitude of a set of the signal harmonics and then apply further computations. These readout schemes realise the parameter estimation by approximating the signal13,15,16, the use of a numeric, non-linear fit algorithm12 and/or the use of a fixed modulation index11,13,15,17 (corresponding to specific absolute arm-length difference in DFMI).

Analytic algorithm

The algorithm presented here operates on a set of signal harmonics which are demodulated individually. Then it uses a recurrence relation between the Bessel functions to provide an analytical estimate of the modulation index m. Knowing m we calculate the Bessel functions and extract them to isolate the interferometric quadratures. And with these quadratures, the phase estimate can be calculated. The algorithm can thereby calculate all parameters of interest without any initial assumption and in a non-iterative fashion for arbitrary input parameters. A sketch of the algorithms as flowchart can be seen in Fig. S.1 in the Supplementary Information with the letters following the names of the paragraphs in this section.

We assume the signal in (1) to be measured on a photo-diode, sampled and digitized with a sampling frequency of fS, leading to a measured time series s(t1),⋯,s(tN) over a time period T which is chosen to be an integer multiple of the modulation time 2π/ωm.

Initial coefficients from the measured signal

We assume that the modulation frequency ωm is well known, if not it could be extracted from the FFT of the time-series, since the measured signal is a sum of ’harmonics’ frequencies n·ωm with n∈Z.

Definition of In and Qn

The first step in our algorithm (similar to10,12) is to calculate the coefficients of the Fourier series of the signal. By I-Q-demodulation of the harmonics of the signal, we define the In and Qn coefficients as:4 In:=1T∫0Tdts(t)·cos(nωmt)·W(t)Qn:=1T∫0Tdts(t)·sin(nωmt)·W(t).

The resulting coefficients of this demodulation are listed in Table 2. Besides the measured signal s(t) and the sin and cos factor of the demodulation, we apply a window function W(t) in (4) which acts as additional low-pass filter to suppress the leakage of the other harmonics of the discreetly sampled data.Table 2 Table of the In and Qn coefficients obtained from the I-Q-demodulation of the different DFM harmonics.

Low-pass filter of	=:	n even	n odd	
s(t)·cos(nωmt)	In	A·Jn(m)·cosφ·cosnψ	-A·Jn(m)·sinφ·sinnψ	
s(t)·sin(nωmt)	Qn	-A·Jn(m)·cosφ·sinnψ	-A·Jn(m)·sinφ·cosnψ	

Filtering of the demodulated signal (Regarding the window function W(t))

Choosing a measurement time T of multiples of 2π/ωm (corresponding to a rectangular window function Wrect(t)=rect(t/T-1/2)), already provides a major suppression for the other harmonics. The Fourier transform of such a window function5 W~rect(ω)∝sinc(ωT/2π)

has zeros at ω=2π/T·k with k as integer and choosing T=2π/ωm ensures that all the other harmonics (with k≠0) are suppressed as shown in Supplementary Fig. S.2 in the Supplementary Information. For real signals, additional filtering may be necessary as the measurement time T or the modulation frequency ωm might not exactly satisfy this condition, as described by Schwarze et al.22. In our reference implementation, we use a squared Hanning window6 W(t)=43T2sin2(πt/T)

(also shown in Fig. S.2) which provides a falloff to higher frequencies in the order of ∝1/f5. The integration operation in (4) acts as additional low-pass filter (which factors in as another 1/f in Fourier space).

Calculating the modulation phase ψ

The calculation of the modulation phase ψ is done identical to the algorithm by Heinzel et al.12. E.g. by calculating the arctan of the ratio In/Qn for the measured n’th harmonics and unwrapping the resulting values one obtains ψ. For completeness we write this as7 arctan(-Qn/In)n evenarctan(In/Qn)n odd.=nψmodπ/2fornψ≥0nψmod-π/2fornψ<0

The resulting values are multiples of ψ modulo 2π. To get ψ, we add π at jumps of the calculated list of (n·ψmodπ) and average the gradient.

Definition of the cn coefficient

Next we eliminate the ψ factor by calculating8 cn:=+In·cosnψ-Qn·sinnψn even-In·sinnψ-Qn·cosnψn odd=A·Jn(m)·cosφn evenA·Jn(m)·sinφn odd

using the previously calculated estimate of ψ.

Calculating the modulation index m

The crucial step in our algorithm is the estimation of the modulation index m. We make use of a recurrence relation (9) for Bessel-functions of the first kind that can be rearranged to yield m.9 Jn-1(m)+Jn+1(m)=2nmJn(m)⇒m=2nJn(m)Jn-1(m)+Jn+1(m)form∈Z.

Since the harmonic amplitudes / the cn coefficients contain sinφ and cosφ factors, one can not directly plug them in (9). Hence we apply the recursive relation twice and solve for m to get the relation10 m=4n(n-1)(n+1)Jn(m)2nJn(m)+(n+1)Jn-2(m)+(n-1)Jn+2(m).

Plugging the cn coefficients obtained from (8) into this relation in place of the Jn, the A and φ dependant factors cancel out and only the Jn(m) factors remain such that11 mn:=4n(n-1)(n+1)cn2ncn+(n+1)cn-2+(n-1)cn+2=mfor∀n.

Next, we average the calculated mn values from every set of 3 harmonics (n-2,n,n+2) to get the m estimate.

For noise dominated harmonics (leading to highly erroneous cn), the values in the square-root of (11) is also highly erroneous and can even become negative. In this case we continue with the absolute value (to calculate any real value from the square root) and effectively discard the erroneous mn later by using a weight that takes the harmonics individual signal-to-noise ratio into account. Without any prior knowledge we would use e.g. the power of the harmonics as weight when averaging over the mn values. To achieve the highest possible precision it is however necessary to use a specific weighting function as specified in Sect. 3.7.

Calculating the interferometric phase φ

Having calculated an estimate for m before, we calculate the estimate of φ straight forward by 3.5.1 Eliminating the m dependency from the cn coefficients via 12 cnJn(mestimate)=A·cos(φ)for n evenA·sin(φ)for n odd

3.5.2 Calculating a weighted average over the even and odd peaks (as explained in 3.7) leading to one value for each quadrature, A·sin(φ) and A·cos(φ) .

3.5.3 Similar to other interferometric phase readout procedures, we calculate the final estimate for φ from the two quadratures via 13 φestimate=arctanA·sin(φ)A·cos(φ)

It is common practice for interferometric measurements to use a atan2 routine (which compares the sign of the quadratures to extend the output range to φ∈(-π,π]). The phase is then calculated from the quadratures via 14 φestimate=atan2(sin(φ),cos(φ))

Additionally, for continuous interferometric measurements, the phase is usually unwrapped using a ’phase-tracking’ algorithm. This further helps to resolve ambiguities at exactly multiples of π/2.

Calculating the amplitude A and offset B

For a clean signal one can estimate the amplitude A and the offset B from the minimal and maximal values. A more elaborate method to estimate A is to extract an estimate for each harmonic and then to average using the optimal weights. Amplitude and offset are often required with less precision and are used to measure effects due to beam tilts that decrease the optical contrast and the constant optical power impinging on each photodiode or a segment of the same.

Averaging over individual peak values

In the algorithm, (almost) every harmonic yields an estimate for m, φ and nψ. To mitigate noise at specific frequencies and to take the individual SNR of the harmonics into account, we found that a crucial part of achieving optimal performance is to perform a weighted average over these values.

A simple weight is to directly use the signal power of the individual harmonics |s~(nωm)|2, which is proportional to ∝N·cn2. For the calculation of the mn which involves three harmonics (with n-2, n and n+2), we use the product of the signal power of these three harmonics as weight for the n’th harmonic. To achieve optimal performance, we use however a different weight given by15 wn=2|Jn(m)|m1+m21+12n(n2-1)-1·|s~(nωm)|2

which is the result of the error-propagation of equation (10) as derived in in Sect. 4.1 and an additional |s~(nωm)| factor. This additional |s~(nωm)|2 factor accounts for the individual signal-to-noise ratio of the harmonics. At certain phase values (e.g. φ≈0 or multiples of π/2), all even (or odd) harmonics contain very little signal energy. In such a case the remaining other odd (or even) harmonics contain most of the signal energy and need to be weighted higher.

A problem of the ideal weights (15) is that they themselves depend on the modulation index m. If no prior information is available; we run the algorithm twice, with the first iteration using only the harmonics signal power as weight; and further iterations use the resulting m value from previous iterations to calculate the weighting function, using the same data for every iteration. In an experimental setup with continuous measurements, we use the values from the previous run to calculate the weights so that the averaging of the mn coefficients only happens once.

Algorithm performance and error analysis

The use case of our algorithm in DFMI has (in general) 5 free parameters (B,A,m,φ,ψ) with ωm being fixed. Our algorithm should be able to deal with arbitrary real values for any of these parameters. In an experimental setup, the parameter space is however restricted to certain ranges. E.g. the sampling and modulation frequency limit the number of resolvable harmonics and any present noise will degrade the readout results.

A more extensive study on the maximally achievable readout performance for various additive white noise contributions can be found in18. The Cramer–Rao lower bound (CRLB) derived there will serve as optimal goal for our algorithms performance analysis. To test the readout algorithm, we simulate a signal with the same parameters as Gerberding et al.10 (including the additive white noise in the order of σ=2·10-7) and let the algorithm calculate the parameters. From the equations derived in Eckhardt & Gerberding18 (as shown in (18) as amplitude spectral density) we can calculate what the minimal achievable noise floor for the individual signal parameter estimates is (the CRLB).16 CRLBm(f)=8σfSNA3+4J0(m)·cos(φ)+J2(2m)·cos(2φ)

17 CRLBψ(f)=4σfSNAm2-m·J1(2m)·cos(2φ)

18 CRLBφ(f)=8σfSNA1-J0(2m)·cos(2φ)

Figure 3 shows the result of this simulation for varying distance ΔL between 1cm and 6m, leading to different modulation indices m and phases φ respectively. For almost the entire range, the phase φ reaches close (smaller than ×2) to the CRLB. At the ’edges’ other effects start to influence the precision. For too small m, not enough harmonics are visible above the noise floor and the algorithm cannot properly run since it needs at least 6 harmonics (3 even and 3 odd) for its calculations. If there are just enough harmonics visible to run the algorithm, a single harmonic with a bad signal-to-noise ratio degrade the performance as they are not averaged out by other higher order harmonics. For too large m, the bandwidth of the measurement (set by the sampling frequency) is filled with harmonics and higher-order harmonics start to spoil via the alias-effect the lower-order harmonics. Increasing the sampling frequency in simulations resolves this issue; in an experimental setup an anti-alias filter (as they are commonly applied before ADCs) similarly mitigates this effect.Fig. 3 Absolute error of the new readout algorithm for varying arm-length difference ΔL for a simulated DFMI setup. The x-axis is shown once in absolute arm length difference ΔL (lower axis) and once in modulation index m=Δω/c0·ΔL as it appears as signal parameter. The other signal parameters are the phase φ=2πΔL/λ and B=1, A=1, ψ=0.1rad, Δf=9GHz, λ=1064nm, fS=2MHz, fm=1kHz, fR=1/T=100Hz and some added Gaussian noise in the order of σ≈2·10-7V/Hz. For each datapoint in the plot we ran the algorithm 100 times with the same signal but with varying noise and calculated the variance. The initial distance value (with m≈2) corresponds to a signal where only a 3 harmonics are visible in the frequency spectrum while the last value (m≈1131) corresponds to the measurement band being filled with harmonics up to the Nyquist frequency (e.g. 1131·fm⪆fs/2). Since no anti-aliasing filter was applied the noise increase above (m≈700) is dominated by aliased harmonics.

Error approximation for m

When calculating the m parameter as outlined in Sect. 3.4 but without using the weighted averaging, the resulting error of the mn estimates can be many orders of magnitude above the CRLB as seen in Fig. 4.Fig. 4 Absolute error of calculated mn values for a simulated signal with the same parameters as in Fig. 3, once with m≈94 (distance ≈0.5 meter) on the left side and once with m≈2.8 on the right side. The upper error prediction (yellow dashed lines) correlates moderately well with the measured error (blue dots). (The error prediction is given by (21) with σ=2·10-7 as before and an additional scaling factor of fm/fS.

Due to the linearity of the equations, any white noise in the signal will lead to the same white noise level for the In, Qn and cn coefficients. For the mn coefficients calculated from equation (10), the linear error propagation, when adding the errors (δcn-2,δcn,δcn+2), yields:19 mn(Jn-2(m)+δcn-2,Jn(m)+δcn,Jn+2(m)+δcn+2)≈m+∑k∈{n-2,n,n+2}∂mn∂Jk(m)·δck+O(δc2)

20 =m+m8Jn(m)[4δcn-m2n(n2-1)((n-1)δcn+2+2nδcn+(n+1)δcn-2)]

21 ≲m+m2|Jn(m)|[1+m21+12n(n2-1)]·δc

When approximating the upper error in (21) we simplify by assuming that δck∈[-δc,δc], which allows us to calculate a simple upper bound for the error. Figure 4 also shows the upper error as yellow dashed line which shows good agreement with the measured errors of the given sample. We find that this linear error estimation already yields a good estimate for the precision of the mn values. Which is why we use the inverse of this error estimator as weight for the individual harmonics as mentioned in (3.7).

In (21) we see that the error of our m estimate scales with m and m3 (ignoring the |Jn(m)| scaling for now) suggesting that small values for m result in smaller errors. The plots in Fig. 4 shows the error of the mn coefficients once for a large m value (left plot) and once for a small value (right plot). In case of a small m, the signal energy is distributed between fewer harmonics compared and the individual harmonics can have a smaller relative error close to the CRLB. The average is however roughly the same for both cases, indicating that there is no general favorable parameter region for m.

Error approximation for φ

Figure 5 is the result of running the readout algorithm with the same parameters as in Fig. 3 but giving the algorithm the exact m and ψ values before calculating the φ parameter. The resulting error is almost the same as in Fig. 3 where the exact parameters were not given before. The increasing error at the edges of the plot can be similarly explained by insufficient harmonics for proper averaging for small m and noise due to aliasing of higher order harmonics for large m. In between, the phase error remains around 2× the CRLB.Fig. 5 Absolute error of calculated φ estimates for varying distances with the same parameters as in Fig. 3 but with the exact m parameter given before calculating φ.

Since the CRLB is the best achievable limit for the given (additive) white noise; we believe the remaining error comes from non-ideal filtering of the higher order harmonics (with power above the white noise) which can leak and add further “noise” to the other In,Qn and cn coefficients. This noise depends on the implementation of the filtering and a stronger filter might improve the results.

Signal dynamics

In this section we consider the effects of a linear phase term δω·t and its effect on the readout algorithm. Such terms can appear i.e. when the mean frequency ω0 of a DFMI laser is not constant, but drifts over time, or when the target motion is so dynamic that the signal is Doppler shifted by the frequency δω. This breaks with the static signal assumption used in some of the previously studied readout algorithms.

We write the dynamic signal as22 smoving(t):=B+Acosm·sinωmt+ψ+φ+δωt

23 ≈B+A∑n∈ZJn(m)cos(nωm+δω)t+nψ+φ.

In a DFMI setup, a target moving with speed v would cause a Doppler shift of the signal of δω=2π/λ0·v. Similarly, the modulation index m also changes to m↦m+Δωv/c·t with v being the target speed and c as speed of light. For small enough speeds v this additional time dependant term is negligible. For large signal dynamics it can become relevant for the estimation of m and influence the parameter estimation of the other coefficients. When simulating dynamic DFM signals for this publication, we always used the full m+Δωv/c·t expression for the signal.

For δω=0 the parts in the sum of (23) with positive and negative indices have the same frequencies, leading to overlapping signals for every nωm harmonic. For δω≠0 this is no longer true as the spectrum appears shifted in the direction of δω. For small enough δω, the harmonics in the PSD appear to split into one peak belonging to the positive indices/frequencies (nωm+δω) and one to the negative indices/frequencies (-nωm+δω) (which are folded onto the positive region when calculating the PSD) as shown in Fig. 6Fig. 6 Signal with ’small’ (δω<ωm/2) frequency shift. The positive and negative frequency parts of the signal harmonics become visible and do no longer overlap.

For the algorithm to work in the presence of a large Doppler shift, an estimator for δω is needed to (a) demodulate the DFM harmonics at their correct frequencies (nωm+δω) and (b) to account for the additional phase shift δωT that is added to φ over the course of the measurement period T. So instead of a single demodulation for each harmonic we implement two that track the splitting tones. The calculation of the parameters and coefficients used in the algorithm change then slightly compared to the previously shown calculation and is shown in Sect. 5.1. In the following we only consider the case where Doppler-shifts are smaller than half of the modulation frequency δω≤ωm/2. Our solution to deal with frequency shifts δω≤ωm/2 in continuous measurements is to initially run the algorithm, without any corrections, and then calculate the gradient of the resulting phase values and use it to generate a linear prediction of the frequency shift δω.

Dynamic readout algorithm (for large δω)

For dynamic signals with δω≠0, the algorithm can be modified to account for the individual tones (nωm+δω and nωm-δω) which were treated as single harmonic before. (Supplementary Fig. S.4 shows a sketch of the demodulation scheme that calculates the additional coefficients for a single (n’th) harmonic). Demodulating at exactly (nωm+δω) and (nωm-δω) yields slightly different I and Q coefficients (↦In,±andQn,±) as defined in (24) and (written in Table S.1 and S.2 in the Supplementary Information).24 In,±:=1T∫0Tdts(t)·cos((nωm±δω)t)·W(t)Qn,±:=1T∫0Tdts(t)·sin((nωm±δω)t)·W(t)

In the limit of δω↦0, the Qn,± and In,± coefficients converge to the previously introduced In and Qn coefficients from Table 2.

Using these coefficients and the additional derived coefficients displayed in Table S.2 (in the Supplementary Information), the algorithms calculation of the individual parameters changes slightly: 5.1.1 for ψ, instead of calculating the Q/I, we calculate (Qn,++Qn,-)/(In,++In,-) and continue with the result as before

5.1.2 for m, we calculate the ψ-free cn,+ and cn,- coefficients as 25 cn,±:=-(Qn,+±Qn,-)cos(nψ)-(In,+±In,-)sin(nψ)forneven-(Qn,+±Qn,-)sin(nψ)+(In,+±In,-)cos(nψ)fornodd

and plug the cn,+ coefficients into equation (11) to continue as before

5.1.3 To obtain the φ estimate, we first eliminate the remaining m (and n) dependency in the cn,± coefficients by calculating 26 d±,n:=2cn,±AJn(m),

which yields for even and odd n: 27 d+,neven=cosφ+cos(δωT+φ)sinc(δωT)d-,neven=sinφ-sin(δωT+φ)sinc(δωT)d+,nodd=sinφ+sin(δωT+φ)sinc(δωT)d-,nodd=cosφ-cos(δωT+φ)sinc(δωT)

Next, we average these coefficients over all even and odd harmonics using our custom weights, leading to exactly 4 averaged coefficients d+,even,d+,odd,d-,even,d-,odd. From these 4 coefficients we calculate the phase via: 28 φ:=arctand+,odd+d-,evend+,even+d-,odd

Precision in case of a (linear) frequency shift δω

Figure 7 shows the simulated performance of the readout algorithm for different relative frequency shifts δω between two consecutive measurements for DFM signals as given by29 sDFM, moving(t)=B+Acosm(t)·sinωmt+ψ+φ+δωt

30 =B+A∑n∈ZJn(m+δm·t)·cos(nωm+δω)t+nψ+φ

with m(t):=m+Δωvc·t=:m+δm·t.

For Doppler shifts below ≈0.6Hz the target is moving so slow that the algorithm consistently reaches the CRLB (green dashed line). For larger Doppler shifts >50Hz the extended algorithm (Sect. 5.1) is still able to achieve close to optimal precision (orange line) while not accounting for the Doppler shift will yield gradually worse results (blue line).

For the ’bump’ between 0.6 and 50 Hz, where the ’static’ and ’dynamic’ algorithm yield the same result but do not reach the CRLB, we found three contributing factors of similar size leading to the higher noise.

Firstly, insufficient filtering during demodulation. Since the harmonics are shifted by δω, they no longer ’sit’ exactly at the zeros of our filter (as described in Sect. 3.1.2) and start to leak into the calculation of the In and Qn coefficients of their neighboring harmonics.

Secondly; the filtering of the splitted tones for the dynamic algorithm. During demodulation of the splitted tones, small enough Doppler shifts are not filtered out sufficiently (like the higher order harmonics) and perturb the calculated In,± and Qn,± coefficients. This is why after the initial decrease in precision (after ≈1Hz in Fig. 7), the precision increases again for the dynamic algorithm as the splitted tones are separated more clearly and filtered out during demodulation of the other tone respectively.

The third effect comes from the DFM signal harmonics containing an additional time dependency in the m(t) parameter. We were yet unable to calculate a closed expression for the In and Qn coefficients when writing the signal with the exact time dependence including m(t)=m0+Δωvt/c. Even with perfectly filtered harmonics, the demodulated harmonic itself has an additional phase error proportional to δm that we disregarded in our algorithm so far.Fig. 7 Algorithm performance for modeled DFM signal (B=0,A=1,ψ=π/4,fs=524kHz,fm=1kHz, fR=16Hz, λ=1550nm and initial distance L0=20cm) with additive white noise (σ=2·10-7) and varying relative frequency shifts δω=2πv/λ. For each point in the plot, we simulated a target moving with constant speed 100 times (but always starting at the same position L0) and averaged the result. The blue line shows the readout algorithm result in case this Doppler shift is ignored and the algorithm runs as described in Sect. 3. The orange line shows the readout algorithm result in case the Doppler shift is known before and algorithm can account for it as explained in Sect. 5.1. The dashed green line is the CRLB showing the highest reachable precision. The vertical purple line marks a Doppler shift of fm/2 (here at 500Hz). At this point two neighboring harmonics (e.g. nωm+δω and (n+1)ωm-δω) overlap and the demodulation of a single isolated harmonic is not possible, leading to the algorithm failing at exactly multiples of fm/2.

Example readout for an oscillating mass

As more realistic example of a measurement, Figure 8 shows the results of the readout algorithm for an oscillating target. Its speed varies over time and the movement contains not only linear but also higher order terms in t. From the time-series (Fig. 8, left plot) we can clearly see that at the points of maximum slope (when the linear speed term is largest and all others small), the ’static’ algorithm version has the largest error, while at the turning points, the deviation from the exact phase is similarly for both ’static’ and ’dynamic extension’ of the algorithm.Fig. 8 Algorithm results for a simulated moving target, oscillating with a frequency of 8 Hz and an amplitude of 2.5 wavelengths (a total path of 5 wavelengths). The DFM signal (fm≈1kHz) was sampled with fS=2MHz and a readout frequency of fR≈8kHz and otherwise the same parameters and additive noise (σ=2·10-7) as in Fig. 7 . The black line marks the exact phase of the target. The blue line is the algorithm output as written in Sect. 3 (without the dynamic extension) and the orange line is the result when using the dynamic extension detailed in Sect. 5.1. The left plot shows one period of the resulting time-series and the right plot shows the amplitude spectral density calculated with an LPSD algorithm23. The dashed red line marks the frequency of the movement of the oscillating target.

It should be noted that the difference in the ASD plot (Fig. 8, right plot) between the static and the dynamic algorithm version is only around a factor of 40 while Fig. 7 shows a difference of around 5 orders of magnitude for the maximum speed of ≈125μm/s (a Doppler shift of δω=2π·125Hz) of the simulated target. The difference between Figs. 7 and 8 is that Fig. 7 is the result of a purely linear movement. This is where the difference between static and dynamic algorithm is greatest. The simulated signal for Fig. 8 has periods of very slow movement, here the static algorithm is just as precise as the dynamic one. Additionally, the movement of Fig. 8 is non-linear and relatively fast and strong compared to the modulation frequency, which causes errors in both the ’static’ and ’dynamic’ algorithm version that are not accounted for, leading again to a similar error in both versions. More optimal implementations of the dynamic algorithm that take into account not only linear speed might mitigate these errors further.

Conclusion and outlook

For the interferometric phase φ our new algorithm achieves close to optimal performance, within a factor of 2 or better of the CRLB. For the modulation index m, our analysis of the linear error propagation shows that it can have an error up to ≈102×CRLB for large m values. The lowest errors for estimating m of ⪅10×CRLB are achieved for m<20.

Besides the basic algorithm in Sect. 3, we provide an extended version detailed in Sect. 5.1 to deal with very dynamic signals (with large Doppler shifts). By using adaptive/dynamic filters it might even be possible to use the algorithm for even larger dynamics than shown in this paper. The precision achieved with our reference implementation does not use any initial values. Especially for larger Doppler shifts (δω) it is possible to use the algorithm e.g. with a close initial estimate for δω, run the algorithm and then calculate a more accurate and potentially larger δω from the result. As long as the frequency difference between this initial estimate and the true value is small enough, even much larger absolute Doppler shifts can be accounted for. For future improvements we aim to add a scheme to improve the algorithms precision over time by, for example, deploying a Kalman filter for the Doppler-shift estimation.

Real DFMI setups contain additional noise sources like 1/f laser frequency noise. The influence of the most relevant of such noise sources have been analysed shortly in the Supplementary Material. In the case of non-white additive noise or the presence of parasitic tones one may adjust the weights of the algorithm to optimise the available SNR. Stray light and ghost beams lead to additional unwanted DFM-like signals overlapping with the ideal signal and potentially cause additional errors. An actual interferometer application also has additional requirements which will need to be met. This includes the ability to provide real-time estimates to run stabilisation control loops. E.g. to lock the average laser frequency to a reference interferometer1 and being able to cope with the presence of parasitic ghost beams24. The latter will require an extension of the algorithm in the presence of more than one beat note. This can be supported by utilizing additional information provided by the readout of two complementary photodiodes at the two outputs of an interferometer, so called balanced detection. In the case of balanced detection or when using quadrant photodiodes to also implement a beam-tilt readout via differential wavefront sensing, the readout of multiple channels can be further optimised, because they will share many of the parameters to be estimated, similar to the optimisation of a phase-locked-loop for heterodyne interferometry by Heinzel et al.25. The core scheme of the algorithm, the estimation of the modulation depth, will also be possible when multiple overlapping signals are present, but it will require either more a-priori parameter knowledge or additional estimation steps.

Lastly, our current implementation is purely sequential and does not use any parallelization of the individual algorithm steps. Unlike an iterative fitting of the signal, the algorithm allows for parallelization e.g. by performing the calculation steps for the used harmonics in parallel. Further speed can be gained by implementing the computational expensive but repetitive calculations of the I and Q coefficients into an FPGA, as they involve a large number of additions and multiplications with constant factors. Using much higher modulation and demodulation frequencies, or achieving higher rates of phase estimation and readout frequency fR to deal with more dynamic signals also requires a less computationally expensive and therefore faster readout algorithm, as presented here.

The algorithm we present will be crucial in the further deployment and development of compact interferometric sensors using DFMI. It relies only on simple analytical calculations, has a deterministic processing time and operates close to optimal noise performance. The analytic nature of the algorithm also make it easily compatible with other readout schemes to provide either initial parameter estimates or to get a low noise parameter estimation, depending on the specific use case.

Supplementary Information

Supplementary Information.

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-70392-9.

Acknowledgements

This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—EXC 2121 “Quantum Universe”—390833306, and by the German Federal Ministry of Education and Research (BMBF, Projects 05A20GU5 & 05A23GU5). We also liked to thank Alexander Franzen for the graphics of his “ComponentLibrary”. We acknowledge financial support from the Open Access Publication Fund Of Universität Hamburg.

Author contributions

T.E. and O.G. contributed equally to this research paper.

Funding

Open Access funding enabled and organized by Projekt DEAL. We acknowledge financial support from the Open Access Publication Fund of Universität Hamburg.

Data availability

The datasets generated and analyzed for this publication are available from the authors on reasonable request at tobias.eckhardt@physik.uni-hamburg.de.

Competing interests

The authors declare no competing interests.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Tobias Eckhardt and Oliver Gerberding.
==== Refs
References

1. Isleif K-S Heinzel G Mehmet M Gerberding O Compact multifringe interferometry with subpicometer precision Phys. Rev. Appl. 2019 12 034025 10.1103/PhysRevApplied.12.034025
Isleif, K.-S., Heinzel, G., Mehmet, M. & Gerberding, O. Compact multifringe interferometry with subpicometer precision. Phys. Rev. Appl. 12, 034025. 10.1103/PhysRevApplied.12.034025 (2019).
2. Yang Y Single-element dual-interferometer for precision inertial sensing Sensors 2020 20 4986 10.3390/s20174986 32899128
Yang, Y. et al. Single-element dual-interferometer for precision inertial sensing. Sensors 20, 4986. 10.3390/s20174986 (2020).32899128
3. Schuldt T Picometer and nanoradian optical heterodyne interferometry for translation and tilt metrology of the LISA gravitational reference sensor Class. Quantum Gravity 2009 26 085008 10.1088/0264-9381/26/8/085008
Schuldt, T. et al. Picometer and nanoradian optical heterodyne interferometry for translation and tilt metrology of the LISA gravitational reference sensor. Class. Quantum Gravity 26, 085008 (2009).
4. van Heijningen, J. V. et al. The payload of the lunar gravitational-wave antenna. J. Appl. Phys. 133, 244501. 10.1063/5.0144687.
5. The LIGO Scientific Collaboration et al. Advanced LIGO. Class. Quantum Gravity. 32, 074001. 10.1088/0264-9381/32/7/074001.
6. Acernese F Advanced virgo: A second-generation interferometric gravitational wave detector Class. Quantum Gravity 2014 32 024001 10.1088/0264-9381/32/2/024001
Acernese, F. et al. Advanced virgo: A second-generation interferometric gravitational wave detector. Class. Quantum Gravity 32, 024001. 10.1088/0264-9381/32/2/024001 (2014).
7. Akutsu T Overview of KAGRA: Detector design and construction history Prog. Theor. Exp. Phys. 2021 2021 05A101 10.1093/ptep/ptaa125
Akutsu, T. et al. Overview of KAGRA: Detector design and construction history. Prog. Theor. Exp. Phys. 2021, 05A101. 10.1093/ptep/ptaa125 (2021).
8. Design report update 2020 for the einstein telescope (2020).
9. van Dongen, J. et al. Reducing controls noise in gravitational wave detectors with interferometric local damping of suspended optics. AIP Publishing 13.
10. Gerberding O Deep frequency modulation interferometry Opt. Express 2015 23 14753 14762 10.1364/OE.23.014753 26072834
Gerberding, O. Deep frequency modulation interferometry. Opt. Express 23, 14753–14762. 10.1364/OE.23.014753 (2015).26072834
11. Smetana, J. et al. Compact michelson interferometers with subpicometer sensitivity. PRApplied 5.
12. Heinzel G Deep phase modulation interferometry Opt. Express 2010 18 19076 19086 10.1364/OE.18.019076 20940802
Heinzel, G. et al. Deep phase modulation interferometry. Opt. Express 18, 19076–19086 (2010).20940802
13. de Groot PJ Vibration in phase-shifting interferometry J. Opt. Soc. Am. A 1995 12 354 10.1364/JOSAA.12.000354
de Groot, P. J. Vibration in phase-shifting interferometry. J. Opt. Soc. Am. A 12, 354. 10.1364/JOSAA.12.000354 (1995).
14. Schwarze, T. S. Phase extraction for laser interferometry in space: Phase readout schemes and optical testing. 10.15488/4233. Tex.ids: schwarze2018a.
15. Sasaki O Okazaki H Sakai M Sinusoidal phase modulating interferometer using the integrating-bucket method Appl. Opt. 1987 26 1089 10.1364/AO.26.001089 20454274
Sasaki, O., Okazaki, H. & Sakai, M. Sinusoidal phase modulating interferometer using the integrating-bucket method. Appl. Opt. 26, 1089. 10.1364/AO.26.001089 (1987).20454274
16. Sudarshanam VS Claus RO Generic j 1...j 4 method of optical phase detection J. Mod. Opt. 1993 40 483 492 10.1080/09500349314550481
Sudarshanam, V. S. & Claus, R. O. Generic j 1...j 4 method of optical phase detection. J. Mod. Opt. 40, 483–492. 10.1080/09500349314550481 (1993).
17. Krauhausen, M., Priem, R., Claßen, R., Prellinger, G. & Pollinger, F. High-resolution absolute range sensors based on the combination of frequency modulation and laser triangulation for heavy industry application.
18. Eckhardt T Gerberding O Noise limitations in multi-fringe readout of laser interferometers and resonators Metrology 2022 2 98 113 10.3390/metrology2010007
Eckhardt, T. & Gerberding, O. Noise limitations in multi-fringe readout of laser interferometers and resonators. Metrology 2, 98–113. 10.3390/metrology2010007 (2022).
19. Wanner G Methods for simulating the readout of lengths and angles in laser interferometers with gaussian beams Opt. Communi. 2012 285 4831 4839 10.1016/j.optcom.2012.07.123
Wanner, G. et al. Methods for simulating the readout of lengths and angles in laser interferometers with gaussian beams. Opt. Communi. 285, 4831–4839. 10.1016/j.optcom.2012.07.123 (2012).
20. Morrison E Meers BJ Robertson DI Ward H Automatic alignment of optical interferometers Appl. Opt. 1994 33 5041 5049 10.1364/AO.33.005041 20935885
Morrison, E., Meers, B. J., Robertson, D. I. & Ward, H. Automatic alignment of optical interferometers. Appl. Opt. 33, 5041–5049. 10.1364/AO.33.005041 (1994).20935885
21. Morrison E Meers BJ Robertson DI Ward Henry Experimental demonstration of an automatic alignment system for optical interferometers Appl. Opt. 1994 33 5037 5040 10.1364/AO.33.005037 20935884
Morrison, E., Meers, B. J., Robertson, D. I. & Ward, Henry. Experimental demonstration of an automatic alignment system for optical interferometers. Appl. Opt. 33, 5037–5040. 10.1364/AO.33.005037 (1994).20935884
22. Schwarze TS Gerberding O Cervantes FG Heinzel G Danzmann K Advanced phasemeter for deep phase modulation interferometry Opt. Express 2014 22 18214 18223 10.1364/OE.22.018214 25089440
Schwarze, T. S., Gerberding, O., Cervantes, F. G., Heinzel, G. & Danzmann, K. Advanced phasemeter for deep phase modulation interferometry. Opt. Express 22, 18214–18223. 10.1364/OE.22.018214 (2014).25089440
23. lpsd algorithm at https://pypi.org/project/lpsd/.
24. Gerberding O Isleif K-S Ghost beam suppression in deep frequency modulation interferometry for compact on-axis optical heads Sensors 2021 21 1708 10.3390/s21051708 33801264
Gerberding, O. & Isleif, K.-S. Ghost beam suppression in deep frequency modulation interferometry for compact on-axis optical heads. Sensors 21, 1708. 10.3390/s21051708 (2021).33801264
25. Heinzel G Álvarez MD Pizzella A Brause N Delgado JJE Tracking length and differential-wavefront-sensing signals from quadrant photodiodes in heterodyne interferometers with digital phase-locked-loop readout Phys. Rev. Appl. 2020 14 054013 10.1103/PhysRevApplied.14.054013
Heinzel, G., Álvarez, M. D., Pizzella, A., Brause, N. & Delgado, J. J. E. Tracking length and differential-wavefront-sensing signals from quadrant photodiodes in heterodyne interferometers with digital phase-locked-loop readout. Phys. Rev. Appl. 14, 054013. 10.1103/PhysRevApplied.14.054013 (2020).
