Beat-Notes Acquisition of Laser Heterodyne Interference Signal for Space Gravitational Wave Detection

Abstract

In space gravitational wave detection missions, the laser heterodyne interference signal (LHI signal) has a high-dynamic characteristic due to the Doppler shift. Therefore, the three beat-notes frequencies of the LHI signal are changeable and unknown. This may further lead to the unlocking of the digital phase-locked loop (DPLL). Traditionally, fast Fourier transform (FFT) has been used as a method for frequency estimation. However, the estimation accuracy cannot meet the requirement of space missions because of the limited spectrum resolution. In order to improve the multi-frequency estimation accuracy, a method based on center of gravity (COG) is proposed. The method improves the estimation accuracy by using the amplitude of the peak points and the neighboring points of the discrete spectrum. For different windows that may be used for signal sampling, a general expression for multi-frequency correction of the windowed signal is derived. Meanwhile, a method based on error integration to reduce the acquisition error is proposed, which solves the problem of acquisition accuracy degradation caused by communication codes. The experimental results show that the multi-frequency acquisition method is able to accurately acquire the three beat-notes of the LHI signal and meet the requirement of space missions.

1. Introduction

Gravitational wave detection in space is based on laser heterodyne interferometry [1]. By measuring the phase change of the LHI signal, the distance change information between different satellite test masses can be calculated [2], so as to achieve gravitational wave signal inversion. However, because the displacement change caused by the gravitational wave signal is very weak, many other functions need to be implemented to ensure the measurement accuracy, for example, inter-satellite ranging and communication, and pilot tone jitter correction [3]. Therefore, the components of the LHI signal are very complicated, including the main carrier beat-note (carrier beat-note), two clock sideband beat-notes (two side beat-notes or side beat-notes), communication codes, and various noises [4]. The distance change information can be obtained by measuring the phase changes of three beat-notes.

The phasemeter is the payload for high-accuracy phase measurement [5]. The phasemeter applies the principle of the digital phase-locked loop to measure the phase change of three beat-notes. However, due to the relative speed between spacecraft, when the laser travels from one spacecraft to another, the Doppler shift [6,7] will shift the frequencies of three beat-notes. Take the LISA as an example [1,2,7]: the frequency range can reach 9 MHz [3], resulting in unknown frequencies of three beat-notes. This may lead to the locking time of the digital phase-locked loop being too long, or even out of lock [8], which will lead to phase information errors. Therefore, before measuring the phase, it is necessary to accurately calculate the frequency of the three beat-notes.

The design of the LISA frequency acquisition algorithm for the carrier beat-note was based on FFT [3]. However, due to the picket fence effect and spectrum leakage of the discrete spectrum [9,10], the frequency error is large [11], and can not meet the requirements of high accuracy. Taiji Plan is a space gravitational wave detection plan proposed by China Academy of Sciences. In the Taiji plan, the current set frequency acquisition accuracy is not less than 30 Hz, and this index may be further improved in the future.

In the fields of harmonic detection [12] and power systems [13,14,15], the interpolated discrete Fourier method has been proposed to analyze exponential signals and estimate frequency [16]. Additionally, the accuracy of frequency estimation is further improved by different interpolation methods [17,18,19,20]. However, these methods are difficult to realize and take up a lot of resources. In addition, based on discrete wavelet packet transformation [21], all-phase [22], decomposition filtering-based dual-window correction algorithms [23], autocorrelation [24], and neural network methods [25], the frequency estimation can be realized. The influence of white noise on frequency estimation has also been analyzed [26,27]. Notably, the COG method is often use for position measurements [28,29]. The real position of the discrete spectrum peak can be obtained by the amplitude ratio of the discrete spectrum. Therefore, applying the COG method after FFT can estimate frequencies. Compared with other methods, the COG method has the characteristics of easy implementation and simple operation, but the estimation accuracy is easily affected by communication codes or noises.

We therefore propose a high-accuracy multi-frequency acquisition method based on the COG method. The method improves the estimation accuracy by using the amplitude of the peak points and the neighboring points of the discrete spectrum. For different windows that may be used for signal sampling, a general expression for multi-frequency correction of the windowed signal is derived. Meanwhile, a method based on error integration to reduce the acquisition error is proposed, which solves the problem of acquisition accuracy degradation caused by communication codes. The proposed method is implemented in VHDL based on the Field Programmable Gate Array (FPGA) platform and experimentally verified.

The structure of the paper is as follows. Section 2 describes the composition of the LHI signal. Section 3 derives a general expression for multi-frequency correction. The error integration method is illustrated in Section 4. Section 5 shows the simulation results. Section 6 describes the implementation of the method on an FPGA platform. Section 7 describes the experimental facilities and test results. Section 8 gives the conclusions.

2. Composition of the LHI Signal

In the LISA, the initial setting of the frequency of carrier beat-note is 11 MHz and the frequencies of side beat-notes are 11 MHz ± 1 MHz. The relative velocities of the two spacecraft can reach 10 ms⁻¹. When the laser wavelength is 1064 nm, the maximum Doppler shift is defined as:

$$\Delta \mathrm{f}=\mathrm{f}\times \left(\frac{\sqrt{\mathrm{c}+\Delta \mathrm{v}}}{\sqrt{\mathrm{c}-\Delta \mathrm{v}}}-1\right)\approx \frac{\Delta \mathrm{v}\cdot \mathrm{f}}{\mathrm{c}}=\frac{\Delta \mathrm{v}}{\mathrm{c}\lambda }=9\mathrm{MHz}$$

(1)

where $\Delta \mathrm{f}$ is the maximum Doppler shift, $\mathrm{f}$ is the laser frequency, $\mathrm{c}$ is the speed of light, $\mathsf{\Delta }\mathrm{v}$ is the relative velocity and $\mathsf{\lambda }$ is the laser wavelength.

Three frequencies of beat-notes are unknown due to Doppler shift, as shown in Figure 1. The frequency of carrier beat-note ranges from 2 MHz to 20 MHz, and the range of frequency variation of the LHI signal is 1 MHz to 21 MHz.

Reference [4] derived the single quadrant voltage formula of the LHI signal output through a Trans-Impedance Amplifier (TIA):

$$\begin{array}{c}\mathrm{V}={\mathrm{J}}_0^2\left({\mathrm{m}}_{\mathrm{s}\mathrm{b}}\right)\sin \left[2\mathsf{\pi }\left({\mathrm{f}}_{\mathrm{m}}+{\mathrm{f}}_{\mathrm{D}}\right)\mathrm{t}+{\mathsf{\phi }}_{\mathrm{m}}+{\mathrm{m}}_{\mathrm{p}\mathrm{r}\mathrm{n}}{\mathrm{c}}_{\mathrm{n}}+\mathrm{n}\left(\mathrm{t}\right)\right]\\ +{\mathrm{J}}_1^2\left({\mathrm{m}}_{\mathrm{s}\mathrm{b}}\right)\left\{\sin \left[2\mathsf{\pi }\left({\mathrm{f}}_{\mathrm{u}}+{\mathrm{f}}_{\mathrm{D}}\right)\mathrm{t}+{\mathsf{\phi }}_{\mathrm{u}}+\mathrm{n}\left(\mathrm{t}\right)\right]+\sin \left[2\mathsf{\pi }\left({\mathrm{f}}_{\mathrm{l}}+{\mathrm{f}}_{\mathrm{D}}\right)\mathrm{t}+{\mathsf{\phi }}_{\mathrm{l}}+\mathrm{n}\left(\mathrm{t}\right)\right]\right\}\end{array}$$

(2)

where ${\mathrm{J}}_{\mathrm{k}}\left(\mathrm{m}\right)$ is the first kind of k-order Bessel function, and ${\mathrm{m}}_{\mathrm{s}\mathrm{b}} \mathsf{\approx } 0.45$ is the phase modulation index. ${\mathrm{f}}_{\mathrm{m}},{\mathrm{f}}_{\mathrm{u}},{\mathrm{f}}_{\mathrm{l}}$ are the frequencies of three beat-notes, respectively, ${\mathrm{f}}_{\mathrm{D}}$ is the Doppler shift, and ${\mathsf{\phi }}_{\mathrm{m}},{\mathsf{\phi }}_{\mathrm{u}},{\mathsf{\phi }}_{\mathrm{l}}$ are the phases of three beat-notes, respectively; ${\mathrm{m}}_{\mathrm{p}\mathrm{r}\mathrm{n}}={0.1}^{} \mathrm{rad}$ is the communication modulation index, and ${\mathrm{c}}_{\mathrm{n}}$ is communication codes composed of binary [–1, 1] sequences; $\mathrm{n}\left(\mathrm{t}\right)$ stands for noise. The amplitude ratio of three beat-notes is about 18:1:1. The logarithmic spectrum of the LHI signal is shown in Figure 2.

3. Multi-Frequency Correction

3.1. Limitations of the Traditional Method

Traditional frequency estimation algorithms are generally based on the FFT, which converts a signal in the time domain into a discrete spectral sequence. The maximum value of the spectral sequence corresponds to the frequency of the signal. The expression is defined as:

$$\mathrm{X}\left(\mathrm{k}\right)={\displaystyle \sum _{\mathrm{n}=0}^{\mathrm{N}-1}{\mathrm{x}}_{\mathrm{n}}{\mathrm{e}}^{-\mathrm{j}\frac{2\mathsf{\pi }\mathrm{n}\mathrm{k}}{\mathrm{N}}}}, 0 \mathsf{\le } \mathrm{k} \mathsf{\le } \mathrm{N}-1$$

(3)

$$\mathrm{f}=\frac{\mathrm{find} \max \left[\mathrm{k},\mathrm{a}\mathrm{b}\mathrm{s}\left(\mathrm{X}\left(\mathrm{k}\right)\right)\right]}{\mathrm{N}}{\mathrm{f}}_{\mathrm{s}}$$

(4)

where $\mathrm{N}$ is the number of sampling points, ${\mathrm{x}}_{\mathrm{n}}$ is the sampling signal, $\mathrm{X}\left(\mathrm{k}\right)$ is the spectrum sequence, $\mathrm{abs}\left(\mathrm{X}\left(\mathrm{k}\right)\right)$ is the spectrum amplitude, ${\mathrm{f}}_{\mathrm{s}}$ is the sampling frequency, and $\mathrm{f}$ is the frequency of ${\mathrm{x}}_{\mathrm{n}}$. The function of $\mathrm{find}\max \left[\mathrm{k},\mathrm{A}\left(\mathrm{k}\right)\right]$ is to find the maximum value of $\mathrm{A}\left(\mathrm{k}\right)$ and return the serial number $\mathrm{k}$ corresponding to the maximum value. The calculation process of FFT generally adopts a base 2 butterfly algorithm, and $\mathrm{N}$ is set to an integer multiple of 2.

However, because the discrete spectrum is the sampling of the real spectrum, only the amplitude at the discrete points can be observed, and the real maximum value is often not obtained. Therefore, the frequency estimation method based on the FFT has error, and the maximum error is half of the spectral resolution, which is defined by:

$$\max \mathrm{ess}=\frac{\mathsf{\delta }\mathrm{f}}{2}=\frac{{\mathrm{f}}_{\mathrm{s}}}{2\mathrm{N}}$$

(5)

where $\max \mathrm{ess}$ is the maximum error and $\mathsf{\delta }\mathrm{f}$ is the spectral resolution.

Because the frequency range of the LHI signal is 21 MHz, according to Nyquist sampling theorem, ${\mathrm{f}}_{\mathrm{s}}$ should be greater than 42 MHz. In order to ensure the quality of the sampling signal, ${\mathrm{f}}_{\mathrm{s}}$ is set to 80 MHz. Additionally, the number of sampling points needs to be greater than 2²¹, because the error of frequency estimation needs to be less than 30 Hz. This leads to a serious waste of resources and a long operation time, and so it needs to be optimized.

3.2. Principles of COG

On the basis of the discrete spectrum, the frequency can be corrected by using the amplitude of the peak points of the spectrum sequence relative to the neighboring points. Therefore, the accuracy of frequency estimation can be improved by introducing a simple calculation after performing FFT with a smaller number of sampling points. The principle is shown in Figure 3.

By using the monotonicity of the main lobe function near the peak point, the corresponding relationship between the ratio function $\mathrm{f}\left(\Delta \mathrm{x}-1\right)/\mathrm{f}\left(\Delta \mathrm{x}\right)$ and $\Delta \mathrm{x}$ is constructed. By calculating the ratio of ${\mathrm{y}}_{\mathrm{k}}$ and ${\mathrm{y}}_{\mathrm{k}+1}$, $\Delta \mathrm{x}$ can be obtained in reverse. The expression is determined by:

$$\mathrm{g}\left(\mathrm{x}\right)=\frac{\mathrm{f}(\mathrm{x})}{\mathrm{f}(\mathrm{x}+1)},\mathrm{x}\in \left[-1,0\right]$$

(6)

$$\Delta \mathrm{x}={\mathrm{g}}^{-1}(\frac{{\mathrm{y}}_{\mathrm{k}}}{{\mathrm{y}}_{\mathrm{k}+1}})$$

(7)

$${\mathrm{f}}_{\mathrm{c}}=\left(\mathrm{k}-\mathsf{\Delta }\mathrm{x}\right)\frac{{\mathrm{f}}_{\mathrm{s}}}{\mathrm{N}}$$

(8)

where, $\mathrm{f}(\mathrm{x})$ is the main lobe function, $\mathrm{g}\left(\mathrm{x}\right)$ is the constructed monotonic function, $\mathrm{k}$ is the smaller index of discrete spectrum peak and the amplitude of the adjacent point. If ${\mathrm{y}}_{\mathrm{k}+1}$ is the discrete spectrum peak, then ${\mathrm{f}}_{\mathrm{c}}$ is the amplitude of the adjacent point. ${\mathrm{f}}_{\mathrm{c}}$ is the corrected value of ${\mathrm{f}}_{\mathrm{c}}$, ${\mathrm{f}}_{\mathrm{c}}$ is the corrected frequency. Additionally, when the value of $\mathsf{\Delta }\mathrm{x}$ is −0.5, ${\mathrm{f}}_{\mathrm{c}}$ is equal to ${\mathrm{y}}_{\mathrm{k}+1}$. In this way, the specific expression for $\mathrm{f}(\mathrm{x})$ is obtained to enable frequency correction.

3.3. Derivation of the Main Lobe Function

The LHI signal contains three beat-notes, and its spectrum is superimposed on the spectrum of three beat-notes, as shown in Figure 4.

By simplifying the LHI signal, the general expression of the superposition signal is obtained:

$$\mathrm{x}\left(\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{x}}_{\mathrm{m}}\left(\mathrm{t}\right)}={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cos \left(2\mathsf{\pi }{\mathrm{f}}_{\mathrm{m}}\mathrm{t}+{\mathsf{\phi }}_{\mathrm{m}}\right)}$$

(9)

where M is the number of signals, and ${\mathrm{f}}_{\mathrm{m}}$, ${\mathrm{A}}_{\mathrm{m}}$ and ${\mathsf{\phi }}_{\mathrm{m}}$ are the frequency, amplitude, and phase of the mth signal, respectively. $\mathrm{x}\left(\mathrm{t}\right)$ is expressed as a complex sine:

$$\mathrm{x}\left(\mathrm{t}\right)=\mathrm{Re}\left[{\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cdot {\mathrm{e}}^{\mathrm{j}\left(2\mathsf{\pi }{\mathrm{f}}_{\mathrm{m}}\mathrm{t}+{\mathsf{\phi }}_{\mathrm{m}}\right)}}\right]$$

(10)

The real signal $\mathrm{x}\left(\mathrm{t}\right)$ can be obtained by ignoring the imaginary part of the complex function $\widehat{\mathrm{x}}$. The complex signal $\widehat{\mathrm{x}}$ is denoted as:

$$\widehat{\mathrm{x}}\left(\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cdot {\mathrm{e}}^{\mathrm{j}\left(2\mathsf{\pi }{\mathrm{f}}_{\mathrm{m}}\mathrm{t}+{\mathsf{\phi }}_{\mathrm{m}}\right)}}$$

(11)

The sampling truncation of $\widehat{\mathrm{x}}\left(\mathrm{t}\right)$, ${\mathrm{f}}_{\mathrm{s}}$ is the sampling frequency, $\mathrm{N}$ is the number of sampling points, the sampled signal is represented as:

$$\widehat{\mathrm{x}}\left(\mathrm{k}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cdot {\mathrm{e}}^{\mathrm{j}\left(2\mathsf{\pi }\frac{{\mathrm{f}}_{\mathrm{m}}}{{\mathrm{f}}_{\mathrm{s}}}\mathrm{k}+{\mathsf{\phi }}_{\mathrm{m}}\right)}},\mathrm{k}\in \left[0,\mathrm{N}-1\right]$$

(12)

Denote ${\mathsf{\omega }}_{\mathrm{m}}$ as the digital angular frequency. ${\mathsf{\omega }}_{\mathrm{m}}$ is determined by:

$${\mathsf{\omega }}_{\mathrm{m}}=2\mathsf{\pi }\frac{{\mathrm{f}}_{\mathrm{m}}}{{\mathrm{f}}_{\mathrm{s}}}$$

(13)

Substitute Equation (13) into Equation (12) and $\mathrm{x}\left(\mathrm{k}\right)$ is defined as:

$$\widehat{\mathrm{x}}\left(\mathrm{k}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cdot {\mathrm{e}}^{\mathrm{j}\left({\mathsf{\omega }}_{\mathrm{m}}\mathrm{k}+{\mathsf{\phi }}_{\mathrm{m}}\right)}},\mathrm{k}\in \left[0,\mathrm{N}-1\right]$$

(14)

Typically, a window is used to sample the signal in order to minimize the spectrum leakage caused by data truncation. Signal sampling without an extra window is equivalent to adding a rectangular window. The general expression for the time domain of a window function is:

$$\mathrm{w}\left(\mathrm{k}\right)={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{I}}{\left(-1\right)}^{\mathrm{i}}\cdot {\mathsf{\alpha }}_{\mathrm{i}}\cdot \cos \left(\frac{2\mathsf{\pi }\mathrm{i}}{\mathrm{N}}\cdot \mathrm{k}\right)},\mathrm{k}\in \left[0,\mathrm{N}-1\right]$$

(15)

$${\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{I}}{\mathsf{\alpha }}_{\mathrm{i}}=1}$$

(16)

where I is the number of terms in the window function, N is the number of window function points (signal sampling points), and ${\mathsf{\alpha }}_{\mathrm{i}}$ is the coefficient of the ith term. For example, ${\mathsf{\alpha }}_0=1$, indicates a rectangular window. ${\mathsf{\alpha }}_0=0.5,{\mathsf{\alpha }}_1=0.5$, indicates a Hann window. ${\mathsf{\alpha }}_0=0.42,{\mathsf{\alpha }}_1=0.5,{\mathsf{\alpha }}_2=0.08$, indicates a Blackman window. ${\mathsf{\alpha }}_0=0.35875,{\mathsf{\alpha }}_1=0.48829,{\mathsf{\alpha }}_2=0.14128,{\mathsf{\alpha }}_3=0.01168$, indicate a Blackman–Harris window.

A windowed signal is expressed as:

$${\mathrm{x}}_{\mathrm{N}}\left(\mathrm{k}\right)=\widehat{\mathrm{x}}\left(\mathrm{k}\right)\cdot \mathrm{w}\left(\mathrm{k}\right)$$

(17)

${\mathrm{x}}_{\mathrm{N}}\left(\mathrm{k}\right)$ is the product of $\mathrm{x}\left(\mathrm{k}\right)$ and $\mathrm{w}\left(\mathrm{k}\right)$ in the time domain, the discrete Fourier transformation of ${\mathrm{x}}_{\mathrm{N}}\left(\mathrm{k}\right)$ is convolved in the frequency domain, and the Fourier transformation of ${\mathrm{x}}_{\mathrm{N}}\left(\mathrm{k}\right)$ is defined as:

$${\mathrm{X}}_{\mathrm{N}}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)=\mathrm{F}[\widehat{\mathrm{x}}\left(k\right)\cdot \mathrm{w}\left(\mathrm{k}\right)]=\frac{1}{2\mathsf{\pi }}\mathrm{X}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)\ast \mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)=\frac{1}{2\mathsf{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\mathrm{X}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\theta }}\right)\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\left(\mathsf{\omega }-\mathsf{\theta }\right)}\right)\mathrm{d}\mathsf{\theta }}$$

(18)

where $\mathrm{F}$ is the Fourier transform function, and $\ast$ is the convolution operator. The spectral amplitude of the windowed signal is expressed as:

$$\left|{\mathrm{X}}_{\mathrm{N}}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)\right|=\left|{\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cdot \mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\left(\mathrm{\omega }-{\mathsf{\omega }}_{\mathrm{m}}\right)}\right)}\cdot {\mathrm{e}}^{\mathrm{j}{\mathsf{\phi }}_{\mathrm{m}}}\right|={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{m}}\cdot \left|\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\left(\mathsf{\omega }-{\mathsf{\omega }}_{\mathrm{m}}\right)}\right)\right|}$$

(19)

The above formula ignores the influence of the initial phase of the signal, so this method is susceptible to phase noise. The spectral function of the window is determined by:

$$\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)={\mathrm{e}}^{-\mathrm{j}\frac{\mathrm{N}}{2}\mathsf{\omega }}\left\{{\mathrm{a}}_0\mathrm{D}\left(\mathsf{\omega }\right)+\frac{1}{2}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{I}}{\left(-1\right)}^{\mathrm{i}}{\mathsf{\alpha }}_{\mathrm{i}}\left[\mathrm{D}\left(\mathsf{\omega }-\frac{2\mathsf{\pi }\mathrm{i}}{\mathrm{N}}\right)+\mathrm{D}\left(\mathsf{\omega }+\frac{2\mathsf{\pi }\mathrm{i}}{\mathrm{N}}\right)\right]}\right\}$$

(20)

where $\mathrm{D}\left(\mathsf{\omega }\right)$ is the Dirichlet kernel, which is defined as:

$$\mathrm{D}\left(\mathsf{\omega }\right)=\frac{\sin \frac{\mathrm{N}\mathsf{\omega }}{2}}{\sin \frac{\mathsf{\omega }}{2}}{\mathrm{e}}^{\mathrm{j}\frac{\mathsf{\omega }}{2}}$$

(21)

The spectrum modulus of $\mathrm{w}\left(\mathrm{k}\right)$ is denoted as:

$$\left|\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)\right|={\mathsf{\alpha }}_0\frac{\sin \frac{\mathrm{N}\mathsf{\omega }}{2}}{\sin \frac{\mathsf{\omega }}{2}}+\frac{1}{2}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{I}}{\mathsf{\alpha }}_{\mathrm{i}}\left[\frac{\sin \left(\frac{\mathrm{N}\mathsf{\omega }}{2}-\mathsf{\pi }\mathrm{i}\right)}{\sin \left(\frac{\mathsf{\omega }}{2}-\frac{\mathsf{\pi }\mathrm{i}}{\mathrm{N}}\right)}+\frac{\sin \left(\frac{\mathrm{N}\mathsf{\omega }}{2}+\mathsf{\pi }\mathrm{i}\right)}{\sin \left(\frac{\mathsf{\omega }}{2}+\frac{\mathsf{\pi }\mathrm{i}}{\mathrm{N}}\right)}\right]}$$

(22)

Substitute Equation (22) into Equation (19). For the lth frequency component, $\mathsf{\Delta }{\mathsf{\omega }}_{\mathrm{i}\mathrm{l}}$ is expressed as the frequency difference between the ith frequency component and the lth frequency component, and the expression for obtaining the lth main lobe function is represented as:

$${\mathrm{f}}_{\mathrm{l}}\left(\mathsf{\omega }\right)={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{i}}\cdot \left|\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\left(\mathsf{\omega }-\mathsf{\Delta }{\mathsf{\omega }}_{\mathrm{i}\mathrm{l}}\right)}\right)\right|}={\mathrm{A}}_{\mathrm{l}}\cdot \left|\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)\right|+{\displaystyle \sum _{\mathrm{i}=1,\mathrm{i}\ne \mathrm{l}}^{\mathrm{M}}{\mathrm{A}}_{\mathrm{i}}\cdot \left|\mathrm{W}\left({\mathrm{e}}^{\mathrm{j}\left(\mathsf{\omega }-\mathsf{\Delta }{\mathsf{\omega }}_{\mathrm{i}\mathrm{l}}\right)}\right)\right|}$$

(23)

$$\mathsf{\Delta }{\mathsf{\omega }}_{\mathrm{i}\mathrm{l}}={\mathsf{\omega }}_{\mathrm{i}}-{\mathsf{\omega }}_{\mathrm{l}}$$

(24)

Let $\mathsf{\omega }=2\mathsf{\pi }\mathrm{x}/\mathrm{N}$, substitute $\mathsf{\omega }$ into Equation (23), the general formula ${\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right)$ is determined by:

$${\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right)={\displaystyle \sum _{\mathrm{r}=1}^{\mathrm{R}}\left({\mathrm{A}}_{\mathrm{r}}\frac{\sin \left(\mathrm{x}\mathsf{\pi }+\frac{\mathrm{N}\mathsf{\pi }\mathsf{\Delta }{\mathrm{f}}_{\mathrm{m}\mathrm{r}}}{{\mathrm{f}}_{\mathrm{s}}}\right)}{\sin \left(\frac{\mathrm{x}\mathsf{\pi }}{\mathrm{N}}+\frac{\mathsf{\pi }\mathsf{\Delta }{\mathrm{f}}_{\mathrm{m}\mathrm{r}}}{{\mathrm{f}}_{\mathrm{s}}}\right)}\left[{\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{I}}{\left(-1\right)}^{\mathrm{i}}\frac{{\mathsf{\alpha }}_{\mathrm{i}}}{\cos \left(\frac{\mathrm{i}\mathsf{\pi }}{\mathrm{N}}\right)}}\right]\right)}$$

(25)

where $\mathrm{R}$ is the number of signal frequencies, ${\mathrm{A}}_{\mathrm{r}}$ is the signal amplitude of the rth frequency component, $\mathrm{N}$ is the number of sampling points, $\mathsf{\Delta }{\mathrm{f}}_{\mathrm{m}\mathrm{r}}$ is the difference between the mth and the rth frequency component, ${\mathrm{f}}_{\mathrm{s}}$ is the sampling frequency, $\mathrm{I}$ is the number of terms of the cosine combined window function, and ${\mathsf{\alpha }}_{\mathrm{i}}$ is the ith coefficient. Theoretically, Equation (25) is suitable for cosine window functions, and it has the same correction accuracy for different window functions. However, due to the impact of data truncation in the calculation process, the number of terms of the window function does not easily become too large.

4. Error Integration

LISA utilizes the direct sequence spread spectrum communication method. It employs the pseudo-random noise code (PRN code) as the communication code to phase modulate the carrier laser. The modulation information is then read out via demodulation, and distance measurement is completed. However, the phase of the signal is unknown during frequency acquisition. The PRN code in the signal acts as a type of phase noise. This creates a deviation between the actual main lobe function and the theoretical function, as shown in Figure 5. Consequently, there will be errors when estimating the actual frequency using the theoretical main lobe function.

In fact, the communication codes used in inter-satellite ranging communication are one or more fixed groups, and the modulation order is fixed. By studying the error of corrected frequency under different sampling starting conditions, it is found that the error of corrected frequency is sinusoidal with the delayed sampling of the signal. Figure 6 shows the relationship between the delayed sampling of the signal and the correction error when the signal is modulated with 2 different groups of communication codes.

By integrating the corrected frequency within the error fluctuation period, the influence of the signal-coupled communication codes or phase noise can be suppressed. The calculation formula is:

$$\stackrel{-}{{\mathrm{f}}_{\mathrm{c}}}=\frac{1}{\mathrm{U}}{\displaystyle {\int }_0^{\mathrm{U}}{\mathrm{f}}_{\mathrm{c}}\left(\mathrm{u}\right)}\mathrm{d}\mathrm{u}=\frac{1}{\mathrm{Q}}{\displaystyle \sum _{\mathrm{q}=1}^{\mathrm{Q}}{\mathrm{f}}_{\mathrm{c}}\left(\mathrm{q}\right)}$$

(26)

where $\mathrm{U}$ is the period of error fluctuation, $\mathrm{Q}$ is the sampling number within one period, and $\stackrel{-}{{\mathrm{f}}_{\mathrm{c}}}$ is the integral corrected frequency.

5. Simulation

The signal shown in Equation (27) is simulated. First, the frequency estimation error of the traditional method is calculated. When the sampling frequency is 80 MHz and the number of sampling points is 65,536, the maximum sampling error is 611 Hz, which is equivalent to half of the spectral resolution. Then, the multi-frequency acquisition method is simulated by using three different window functions, and the acquisition errors of the three methods are similar, which verifies the universality of the proposed method. Finally, the maximum frequency error with and without error integration is compared. The accuracy can be further improved by using error integration, and the probability of the error being less than 1 Hz is about 90 percent.

$$\mathrm{s}\left(\mathrm{k}\right)=0.9\sin \left(2\mathsf{\pi }\frac{{\mathrm{f}}_{\mathrm{m}}}{{\mathrm{f}}_{\mathrm{s}}}\mathrm{k}+0.1{\mathrm{m}}_{\mathrm{p}\mathrm{r}\mathrm{n}}\right)+0.05\sin \left(2\mathsf{\pi }\frac{{\mathrm{f}}_{\mathrm{u}}}{{\mathrm{f}}_{\mathrm{s}}}\mathrm{k}\right)+0.05\sin \left(2\mathsf{\pi }\frac{{\mathrm{f}}_{\mathrm{l}}}{{\mathrm{f}}_{\mathrm{s}}}\mathrm{k}\right)$$

(27)

$${\mathrm{f}}_{\mathrm{m}}=2\mathrm{MHz}+\mathsf{\Delta }\mathrm{f}\times \left(\mathrm{M}-1\right)$$

(28)

where ${\mathrm{f}}_{\mathrm{m}},{\mathrm{f}}_{\mathrm{u}},{\mathrm{f}}_{\mathrm{l}}$ are the frequency of the three beat-notes, ${\mathrm{f}}_{\mathrm{u}}={\mathrm{f}}_{\mathrm{m}}+1\mathrm{MHz}$, ${\mathrm{f}}_{\mathrm{l}}={\mathrm{f}}_{\mathrm{m}}-1\mathrm{MHz}$, ${\mathrm{f}}_{\mathrm{s}}$ = 80 MHz, N = 65,536, ${\mathrm{m}}_{\mathrm{p}\mathrm{r}\mathrm{n}}$ is communication codes. The method was tested ${\mathrm{f}}_{\mathrm{m}}$ from 2 MHz to 19,997,989.4 Hz. Additionally, the number of test points M = 3383. $\mathsf{\Delta }\mathrm{f}$ = 5321.7 Hz.

Table 1. The maximum error of the frequency acquisition simulation.

Method		Error of fm	Error of fu	Error of fl
FFT		610.5531 Hz	610.5632 Hz	610.7658 Hz
frequency correction	Hann	1.5336 Hz	25.7858 Hz	30.6588 Hz
	Blackman	1.8053 Hz	28.5197 Hz	43.9575 Hz
	Blackman–Harris	2.3936 Hz	35.9808 Hz	37.8291 Hz
error integration	Hann	0.3325 Hz	4.0745 Hz	4.5455 Hz
	Blackman	0.2945 Hz	4.0925 Hz	4.1662 Hz
	Blackman–Harris	0.3708 Hz	2.9534 Hz	3.3880 Hz

shows the maximum error of the frequency acquisition simulation.

The simulation results without error integration are shown in Figure 7. The acquisition errors of three kinds of window function are similar, the maximum acquisition error of the carrier beat-note is about 2.4 Hz, and the maximum acquisition error of the two side beat-notes is about 36 Hz. Compared with 610 Hz of the FFT method, the precision of two side beat-notes is improved by an order of magnitude and the carrier beat-note is improved by two orders of magnitude. This method is only compared with the FFT method because it is an improvement on FFT. Of course, this method is not the most accurate method, but it is suitable for space gravitational wave detection. This is because this method requires less computation, so it is suitable for real-time processing.

The simulation results of error integration are shown in Figure 8. The acquisition errors under the three window functions are similar. The maximum acquisition error of the carrier beat-note is below 0.4 Hz, and the maximum acquisition error of the two side beat-notes are below 4.6 Hz. Compared with the simulation results without error integration, the acquisition accuracy is increased by about an order of magnitude. Additionally, the probability of two side beat-notes between −1 Hz and 0 is about 0.9. Figure 9 shows the probability distribution of frequency acquisition errors, where P represents the probability.

6. Implementation of the Algorithm

The multi-frequency acquisition algorithm consists of clock module, communication module, signal windowing sampling (SWS) module, discrete spectral amplitude calculation (DSAC) module, peak search (PS) module, frequency correction (FC) module, and error integration (EI) module, as shown in Figure 10.

The clock module divides the 20 MHz external clock into the 80 MHz system clock and 320 MHz AD sampling clock; the communication module uses an RS 232 serial port to read the operation instructions sent by the PC and receive the calculation results for the three frequencies; the SWS module adopts the ADI’s high-speed ADC chip, with a sampling rate of 80 MHz and a 16-bit digital signals output, and a serial data rate of 320 Mbps. The window is Hanning, and the sequence of 16-bit window functions is stored in ROM. Windowing is realized by multiplication, and the 32-bit windowed signal is the output. The DSAC module uses the existing algorithm packages: the FFT algorithm and CORDIC algorithm calculate the modulus of the FFT results; the PS module calculates the peak value and peak index of the spectrum by peak search algorithm, and the flow chart of the algorithm is shown in Figure 11. The FC module stores the frequency correction curve in ROM, and calculates the inverse solution of the amplitude ratio using a division algorithm. Finally, the EI module further reduces the acquisition error based on Equation (26).

7. Experimental Facilities and Test Results

The experimental facilities of the frequency acquisition experiment included an LHI signal simulation system for space gravitational wave detection, a frequency acquisition algorithm experimental circuit board, an external reference clock, a power supply, and a PC, as shown in Figure 12.

The signal simulation system was used to generate the simulated LHI signal. The signal is shown in Equation (27). ${\mathrm{f}}_{\mathrm{m}}$ ranged from 2 MHz to 20 MHz. The number of test frequency M was 34. The frequency change interval $\mathsf{\Delta }\mathrm{f}$ was 532,170 Hz. Each frequency was acquired 10 times. Figure 13 shows the results of the frequency acquisition experiment. The maximum error of the carrier beat-note was less than 1 Hz. Additionally, the maximum error of two side beat-notes was less than 10 Hz. The acquisition time was 125 ms. Comparing Figure 8 and Figure 13, the experimental acquisition error was slightly larger than the theoretical simulation value. This was caused by the experimental circuit and the truncation error of AD sampling, and further analysis of this influence is needed in the future. At present, the maximum acquisition error meets the task requirement of less than 30 Hz. However, in the future, with the improvement of index requirements, higher requirements may be put forward for acquisition accuracy.

8. Conclusions

In this paper, a multi-frequency acquisition metod of LHI signals is presented, and a system platform was built. A signal simulation system was used to output an LHI simulation signal, which was analyzed using the frequency acquisition method in FPGA, and then the frequency was obtained using a PC. In addition, the simulation results of frequency acquisition under different window functions show that the frequency acquisition errors are basically the same, which verifies the universal applicability of this method. The experimental results show that the method can reduce the fluctuation of communication codes, and acquire signals from 1 MHz to 21 MHz. Compared with FFT, the maximum acquisition error of this method is reduced to 10 Hz, which meets the requirements of current frequency acquisition tasks. However, due to the noise of the experimental circuit or the truncation error of AD sampling, the acquisition error of the experimental results was slightly larger than that of the theoretical simulation. Therefore, in future work, it is necessary to analyze the influence of noise and clarify the influence mechanism, so as to apply the method better. In addition, we need to further optimize the structure and computational complexity of the method, so as to improve the computational efficiency in engineering while maintaining accuracy. Finally, it is worth noting that the method is not only suitable for space gravitational wave detection, but also for other occasions of high-accuracy acquisition of multi-frequency signals, such as harmonic detection of power signals.