Digital signal processsing functions for ultra- low frequency calibrations ACTA IMEKO | www.imeko.org December 2020 | Volume 9 | Number 5 | 374 ACTA IMEKO ISSN: 2221-870X December 2020, Volume 9, Number 5, 374 - 378 DIGITAL SIGNAL PROCESSSING FUNCTIONS FOR ULTRA- LOW FREQUENCY CALIBRATIONS Henrik Ingerslev1, Søren Andresen2, Jacob Holm Winther3 1 Brüel & Kjær, Danish primary laboratory for acoustics, Denmark, henrik.ingerslev@bksv.com 2 Hottinger, Brüel and Kjær, Denmark, soren.andresen@bksv.com 3 Brüel & Kjær, Danish primary laboratory for acoustics, Denmark, jacobholm.winther@bksv.com Abstract: The demand from industry to produce accurate acceleration measurements down to ever lower frequencies and with ever lower noise is increasing [1][2]. Different vibration transducers are used today for many different purposes within this area, like detection and warning for earthquakes [3], detection of nuclear testing [4], and monitoring of the environment [5]. Accelerometers for such purposes must be calibrated in order to yield trustworthy results and provide traceability to the SI-system accordingly [6]. For these calibrations to be feasible, suitable ultra low-noise accelerometers and/or signal processing functions are needed [7]. Here we present two digital signal processing (DSP) functions designed to measure ultra low-noise acceleration in calibration systems. The DSP functions use dual channel signal analysis on signals from two accelerometers measuring the same stimuli and use the coherence between the two signals to reduce noise. Simulations show that the two DSP functions are estimating calibration signals better than the standard analysis. The results presented here are intended to be used in key comparison studies of accelerometer calibration systems [8][9], and may help extend current general low frequency range from e.g. 100mHz down to ultra-low frequencies of around 10mHz, possibly using somewhat same instrumentation. Keywords: low-noise; coherent power; coherent phase; calibration; dual-channel; ultra-low frequencies 1. INTRODUCTION In the field of dual channel signal analysis there are some very powerful functions for analysing signals, such as the well-known frequency response function and coherence. But there are also other functions like the coherent power function (COP) and non-coherent power function which are very powerful for decomposing noisy signals into the coherent part and the non-coherent part [10][11][12]. Consider an accelerometer calibration setup with two accelerometers mounted close to each other and measuring the same stimuli. They will both measure the acceleration of the shaker, but since they are different sensors with different conditioning, they will have different noise. The two signals will have a coherent part which is the acceleration signal and a noncoherent part which is the noise. Hence, the COP can be a powerful tool for extracting the signal from the noise and thereby increase the measuring accuracy of the power. A similar function for increasing the measuring accuracy of the phase by separating the coherent phase from the non-coherent phase is also derived in the next section, called the coherent phase (or argument) function (COA). For the COA to work in a proper manner, it is crucial that the signal applied to the shaker is a continuous signal, like a sine or a multi-sine, and that the frequencies of the sines are very precise and phase synchronized with the frequencies of the Fourier transformation, to prevent the phase from drifting or even make jumps. More details on this will be given in section 1.2. The two DSP functions analysed in this article may prove relevant to be used for e.g. key comparison of calibration systems down to extremely low frequencies of around 10mHz where noise becomes a real challenge [7][8][9]. The degree to which the COP, and the COA can separate a signal into coherent and noncoherent parts increases with the length of the measurement, and generally depends on parameters like how many time-samples the measurement is divided into, how long each time-sample is, the sampling rate, and the Fourier transform used. 2. DIGITAL SIGNAL PROCESSING FUNCTIONS - THEORY In this section the theory for two DSP functions is outlined. The two functions are based on dual http://www.imeko.org/ mailto:henrik.ingerslev@bksv.com mailto:soren.andresen@bksv.com mailto:jacobholm.winther@bksv.com ACTA IMEKO | www.imeko.org December 2020 | Volume 9 | Number 5 | 375 channel signal analysis and can give a better estimate of signals in very noisy environments, than standard signal analysis. The first function is the COP for estimating the power or amplitude of the signal. And the second is the COA for estimating the phase. Both functions rely on the coherence between the two signals. 1.1. Coherent power function Consider two sensors both measuring the same stimuli and positioned close enough for their mutual transfer function to be considered unity. As illustrated in Figure 1 the signal without noise called 𝑢(𝑡), and the noise from each sensor called 𝑛(𝑡) and 𝑚(𝑡) yields the output signals from the two sensors, called 𝑎(𝑡) and 𝑏(𝑡). Now consider 𝑗 = 1…𝑁 discrete time-samples measured with the two sensors: 𝑎𝑗(𝑡𝑖) = 𝑢𝑗(𝑡𝑖) + 𝑛𝑗(𝑡𝑖) (1) 𝑏𝑗(𝑡𝑖) = 𝑢𝑗(𝑡𝑖) + 𝑚𝑗(𝑡𝑖) (2) Here 𝑡𝑖 is the discrete time in each time-sample, 𝑢𝑗 is the discrete signal, and 𝑛𝑗 and 𝑚𝑗 are discrete noise in each time-sample. By discrete Fourier transformation of equation (1) and (2) we get 𝐴𝑗(𝑓𝑘) = 𝑈𝑗(𝑓𝑘) + 𝑁𝑗(𝑓𝑘) (3) 𝐵𝑗(𝑓𝑘) = 𝑈𝑗(𝑓𝑘) + 𝑀𝑗(𝑓𝑘) (4) where 𝑓𝑘 is the discrete frequency, 𝑈𝑗, 𝑁𝑗, and 𝑀𝑗 are the discrete Fourier transforms of 𝑢𝑗, 𝑛𝑗, and 𝑚𝑗 respectively. Figure 1: Illustration of the signals used in the dual channel signal analysis and derivation of the two DSP functions. 𝑢(𝑡) is the signal we want to measure, and 𝑛(𝑡) and 𝑚(𝑡) are the noise contributions from each sensor, which then yields the two output signals 𝑎(𝑡) and 𝑏(𝑡). The cross spectrum is given by 𝑆𝐴𝐵(𝑓𝑘) = 1 𝑁 ∑ 𝐴𝑗(𝑓𝑘)𝐵𝑗 ∗(𝑓𝑘) 𝑁−1 𝑗=0 (5) for 𝑁 → ∞, and * denotes complex conjugate. By inserting equation (3) and (4) in equation (5), and using that 𝑈𝑗 , 𝑁𝑗 , and 𝑀𝑗 are uncorrelated the cross spectrum can be given by 𝑆𝐴𝐵(𝑓𝑘) = 1 𝑁 ∑ 𝑈𝑗(𝑓𝑘)𝑈𝑗 ∗(𝑓𝑘) 𝑁−1 𝑗=0 ≝ 𝑆𝑈𝑈(𝑓𝑘) (6) where 𝑆𝑈𝑈(𝑓𝑘) is the power spectrum of the signal without noise, i.e. the coherent power. Therefore, the COP can in this setup be given by: 𝐶𝑂𝑃(𝑓𝑘) = 𝑆𝐴𝐵(𝑓𝑘) (7) 1.2. Coherent phase function Consider the following function, for 𝑁 → ∞: 𝐷𝐴𝐵(𝑓𝑘) = 1 𝑁 ∑ 𝐴𝑗(𝑓𝑘)𝐵𝑗(𝑓𝑘) 𝑁−1 𝑗=0 (8) It looks like the cross spectrum from equation (5), but without the complex conjugation. This “non- conjugated cross spectrum” is very useful for deriving a function for measuring the phase of the coherent signal. We can similarly to the derivation of the COP insert equation (3) and (4) in (8), and use that 𝑈𝑗, 𝑁𝑗, and 𝑀𝑗 are uncorrelated. Hence, the non-conjugated cross spectrum can now be given by, where we have omitted the 𝑓𝑘 dependence to make room. 𝐷𝐴𝐵 = 1 𝑁 ∑ 𝑈𝑗𝑈𝑗 𝑁−1 𝑗=0 (9) = 1 𝑁 ∑|𝑈𝑗| 2 exp⁡(2𝑖∠𝑈𝑗) 𝑁−1 𝑗=0 (10) ≃ 1 𝑁 ∑|𝑈𝑗| 2 𝑁−1 𝑗=0 exp(2𝑖 1 𝑁 ∑ ∠𝑈𝑗 𝑁−1 𝑗=0 ) (11) = 𝑆𝑈𝑈exp(2𝑖∠𝑈̅̅ ̅̅ ) (12) From equation (10) to (11) we have approximated the summation of vectors with length |𝑈𝑗| 2 and angle 2∠𝑈𝑗 by vectors with correct length but all with the mean angle 2∠𝑈̅̅ ̅̅ . Figure 2 illustrates this approximation and shows that for relatively small changes in angle from sample to sample this is a good approximation. In calibration applications the signal measures the acceleration of the shaker. And by applying a sine or multi-sine with frequencies exactly the same and phase synchronized with the Fourier frequencies to the shaker, the phase from sample to sample can be kept steady without drifting and the approximation will therefore be very good in such calibration applications. http://www.imeko.org/ ACTA IMEKO | www.imeko.org December 2020 | Volume 9 | Number 5 | 376 In equation (12) ∠𝑈̅̅ ̅̅ is the mean phase of the coherent signal. Therefore, the coherent phase function can be given by 𝐶𝑂𝐴 = ∠𝑈̅̅ ̅̅ . And by rewriting equation (12), the COA can be given by: 𝐶𝑂𝐴(𝑓𝑘) = 1 2 Imag(ln(𝐷𝐴𝐵(𝑓𝑘))) (13) We have in the derivation of equation (13) used that the signal power 𝑆𝑈𝑈 is purely real and can therefore be omitted. Figure 2: Schematic graphs illustrating the approximation from equation (10) to (11).The arrows represent 𝑈𝑗 and the x- and y- axes are the real and imaginary part of 𝑈𝑗 (a) shows the correct summation where each arrow has length |𝑈𝑗| 2and angle 2∠𝑈𝑗 as in equation (10), and (b) shows the approximative summation where each arrow has correct length but a mean angle 2∠𝑈̅̅ ̅̅ as in equation (11). As can be seen from (a) to (b), if 𝑈𝑗 does not change to much between samples the approximation is good. 3. SIMULATIONS In this section we test the COP and the COA functions on simulated data. We calculate the discrepancy between the functions estimate of the signal amplitude and phase, and the true values. And we compare this with standard signal analysis which is to average the amplitude and phase over all the samples. When measuring the cross spectrum or the non- conjugated cross spectrum for finding the COP and the COA we only measure a finite number of samples 𝑁, which therefore only yields an estimate of the COP and the COA. Hence, the more samples the more precise the estimate will be, and in the following the simulations is based on a 102400 s (~ 28 h) time sample divided into 𝑁 = 1024 samples of 100 s each, with 2048 discrete measurement points in each sample, and a Fourier transform from 10mHz to 10.24 Hz with a 10 mHz step. How well the COP and COA works is estimated by representing 𝑎𝑗, and 𝑏𝑗 from equation (1) and (2) with simulated data. Here the signal 𝑢𝑗 is a multi-sine with 𝑙 = 1…𝑀 frequencies (𝑓𝑙) and phases (𝜙𝑙), all with amplitude 𝑈0: 𝑢𝑗(𝑡𝑖) = 𝑈0∑sin⁡(2𝜋𝑓𝑙𝑡𝑖 + 𝜙𝑙) 𝑀 𝑙=1 (14) The noise 𝑛𝑗, and 𝑚𝑗 are random generated white noise with amplitude 𝑁0. And the signal to noise ratio is given by: 𝑆𝑁𝑅 = 𝑈0 𝑁0 (15) The frequencies (𝑓𝑙) of the multi-sine in equation (14) must be precisely the same or synchronized to a subset of the Fourier transformation frequencies (𝑓𝑘) from equation (3) and (4), otherwise the COA will drift. This requirement is easily met in the simulations presented here, since the subset of frequencies (𝑓𝑙) can be set to be identical to some of the frequencies (𝑓𝑘). But in real measurements this requirement might be challenging to meet. Figure 3: Based on simulated data the coherent power function and coherent phase function is tested for its strength for estimating the amplitude and phase of sine waves in noisy data. The graph shows the deviation of the amplitude and phase from the true vales. (a) shows the mean amplitude, ⁡0.5(𝐴 + 𝐵), and the coherent amplitude, √𝐶𝑂𝑃 . (b) shows the mean phase 0.5(∠𝐴 + ∠𝐵), and the coherent phase 𝐶𝑂𝐴. Figure 3(a) shows the discrepancy between the signal amplitudes 𝑈0 from equation (14) and the coherent amplitude, defined as √𝐶𝑂𝑃(𝑓𝑘) , for a signal to noise ratio of 𝑆𝑁𝑅 = 0.32 in red circles. And for comparison we also plot the mean amplitude in blue circles, that is, 0.5(𝐴(𝑓𝑘) + 𝐵(𝑓𝑘)) where 𝐴(𝑓𝑘) = ∑|𝐴𝑗(𝑓𝑘)| and 𝐵(𝑓𝑘) = ∑|𝐴𝑗(𝑓𝑘)| is the average amplitudes over all samples. And we plot only at the frequencies of the signal, i.e. at 𝑓𝑘 = 𝑓𝑙. It is clearly seen that the COP function estimates the amplitude better that the mean amplitude in the full frequency range. http://www.imeko.org/ ACTA IMEKO | www.imeko.org December 2020 | Volume 9 | Number 5 | 377 Similarly, Figure 3(b) shows the discrepancy between the phase 𝜙𝑙 and the COA for a signal to noise ratio of 𝑆𝑁𝑅 = 0.024 in red circles. And for comparison we also plot the mean phase defined as 0.5(∠𝐴 + ∠𝐵) where ∠𝐴 = ∑arg⁡(𝐴𝑗(𝑓𝑘)) and ∠𝐵 = ∑arg⁡(𝐵𝑗(𝑓𝑘)) are the average phases over all samples. And it is seen that the COA estimates the phase better that the mean phase. Figure 4:Simulated deviation in amplitude and phase versus signal to noise ratio, in the upper graph for the mean amplitude and the coherent amplitude defined as √𝐶𝑂𝑃, and in lower graph for the mean phase and the coherent phase COA. Each deviation point plotted here is an average over all the frequencies from Figure 3. The deviations plotted in the four curves in Figure 3 seems to be independent of frequency. Therefore, by averaging the deviation for each curve in Figure 3 over all the frequencies 𝑓𝑙, we get a mean deviation for the COP and the mean amplitude at 𝑆𝑁𝑅 = 0.32, and a mean deviation for the COA and the mean phases at 𝑆𝑁𝑅 = 0.024. We have done this for a range of signal to noise ratios from 𝑆𝑁𝑅 = 0.01 to 𝑆𝑁𝑅 = 10 and plotted it in Figure 4. It shows that the COP is better than the mean amplitude from about 𝑆𝑁𝑅 = 1, and the COA is better than the mean phase from about 𝑆𝑁𝑅 = 0.04. This shows that the “critical” SNR where the COP and COA functions has lower deviation than the simple mean values is different for COP and COA. And though the data in Figure 4 depends highly on the length of the time sample and Fourier transform used, this trend of a lower critical SNR for the COA function than the COP functions were always seen and needs further investigations. 4. SUMMARY We have derived two DSP functions for very accurate measurements of amplitude and phase from two accelerometers measuring the same stimuli. We have tested the two DSP functions on simulated data and our findings based on the simulations shows promise to the functions as good tools for accurately measuring amplitudes and phases of a multi-sine wave in a noisy environment. These findings may prove useful for key comparisons of accelerometer calibration systems down to ultra-low frequencies, since for such measurements noise becomes a huge problem as the frequencies approaches 10mHz. Hence, by replacing the accelerometer used in key comparisons by two accelerometers and by using the two DSP functions described here, the frequency range in key comparisons may be possible to extend down to ultra- low frequencies of around 10mHz. 5. ACKNOWEDGEMENTS We would like to acknowledge Torben Rask Licht and Lars Munch Kofoed for fruitful discussions. We would also like to acknowledge the EU EMPIR project “Metrology for low-frequency sound and vibration”, 10ENV03 Infra-AUV, as the initial reason for the approach on optimising method of analysis in noisy environments, especially found relevant when doing low frequency calibration. 6. REFERENCES [1] T. Bruns, S. Gazioch: “Correction of shaker flatness deviations in very low frequency primary accelerometer calibration, IOP Metrologia, vol. 53, no. 3, pp. 986 (2016). [2] J. H. Winther, T. R. Licht, “Primary calibration of reference standard back-to-back vibration transducers using modern calibration standard vibration exciters”, Joint Conference IMEKO TC3, TC5, TC22, Helsinki, (2017). [3] G. Marra, C. Clivati et al., “Ultra-stable laser interferometry for earthquake detection with terrestrial and submarine optical cables”, Science, vol. 361, Issue 6401, pp. 486-490, (2018). [4] P. Gaebler, L Ceranna et al., “A multi- technology analysis of the 2017 North Korean nuclear test”, Solid Earth, 10, 59–78, (2019). [5] R. A. Hazelwood, P. C. Macey, S. P. Robinson, L. S. Wang, “Optimal Transmission of Interface Vibration Wavelets - A Simulation of Seabed Seismic Responses”, J. Mar. Sci. Eng. 6, 61-79, (2018). [6] L. Klaus, M. Kobusch, Seismometer Calibration Using a Multi-Component Acceleration Exciter, IOP Conf. Series: Journal of Physics: Conf. Series 1065 (2018). http://www.imeko.org/ ACTA IMEKO | www.imeko.org December 2020 | Volume 9 | Number 5 | 378 [7] EU EMPIR project “Metrology for low- frequency sound and vibration”, (2020). [8] BIPM Key comparison, CCAUV.V-K5 (2017), https://www.bipm.org/kcdb/comparison?id=453 [9] EURAMET Key comparison, AUV.V-K5 (2019), https://www.bipm.org/kcdb/comparison?id=1631 [10] H. Herlufsen, “Dual channel FFT analysis part I”, B&K Technical Review, no.1, (1984). [11] J. S. Bendat, A. G. Piersol, “Random Data”, Wiley-Interscience, (1986). [12] J. S. Bendat, A. G. Piersol, “Engineering applications of correlation and spectral analysis”, Wiley, New York, (1993). http://www.imeko.org/ https://www.bipm.org/kcdb/comparison?id=453 https://www.bipm.org/kcdb/comparison?id=1631