Wavelet-based saturated absorption line detection for laser frequency locking

Owing to the presence of noise and the Doppler background, accurate saturated absorption (SA) peak automatic identification technology poses a significant challenge for laser frequency tuning and locking. To address this issue, a novel peak identification algorithm for the SA spectrum is proposed. First, a Gaussian filter based on a Gaussian continuous wavelet transform is proposed to mitigate the spectral high-frequency noise. Subsequently, a hybrid method combining a first-order Gaussian continuous wavelet transformation and adaptive threshold judgment was designed for multi-peak boundary segmentation. Finally, we obtained the target peak and its sweeping voltage based on an adaptive nonlinear fitting algorithm, which was almost unaffected by the peak asymmetry caused by the Doppler background.


Introduction
Atomic gyroscopes and magnetometers represent a new generation of highly precise quantum sensors, which have attracted increasing attention [1,2].In quantum-sensing systems, lasers are crucial pump and detector light sources [3].However, frequency fluctuations in lasers can deteriorate the atomic pumping rate and polarization in an atomic sensor system [4,5].The laser frequency is locked to a fixed-frequency reference, such as the resonance frequencies of a high-Q Fabry-Perot (FP) cavity [6] or the transition lines of atoms.Despite ultra-stable cavities being usually made of ultralow thermal Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.expansion coefficient materials and equipped with expensive high-vacuum systems, they are unsuitable for integrating a quantum sensor prototype [7].With the development of integrated photonics technology, miniaturized Brillouin lasers based on integrated optical waveguides have achieved hertzlevel linewidth outputs, providing new prospects for atomic sensors.However, integrated Brillouin lasers cannot simultaneously balance the frequency and high-power output to achieve efficient pumping rate of inertial atoms [8].And the thermal stability of a microwave waveguide resonator is another challenge in maintaining the long-term stability of the laser frequency [9].
Furthermore, the atomic saturated absorption (SA) transition line generated by the intrinsic properties of atoms and their features is commonly used in atomic sensing devices [10,11].The localization of SA peaks is achieved through manual observation with the assistance of an oscilloscope, which cannot guarantee the long-term stability of the system because it cannot automatically detect and relock [12,13].Therefore, exploring an automatic peak detection method is critical for ensuring long-term stability.
Peak detection technology has been extensively studied in various fields.For instance, peak recognition algorithms have been developed for sensing pressure or temperature using optical fiber sensors [14,15].R-peak and T-peak-seeking algorithms of electrocardiograms have been studied for cardiac anomaly surveillance using moving average filters, discrete wavelet transform-modulus maxims, and other methods [16,17].However, the peak recognition method is usually proposed based on corresponding spectral features.The SA peaks are distinct from those described above.Multiple peaks are present in a single sweep, and an irregular Doppler background is superimposed on the SA spectrum.Furthermore, the amplitudes of some weak peaks are much smaller than those of the Doppler background peak; thus, identifying a specific atomic transition peak among the numerous true and interfering false peaks is challenging.
To mitigate the deleterious effects of the Doppler background, Guo et al [18] proposed a technique for determining the SA peak based on an error signal.The error signal is obtained using the lock-in-amplification technology; therefore, a phase shift is inevitably introduced during signal processing.Li designed vector machines to identify peak features [19].Furthermore, this method eliminated the interference of the Doppler background superimposed on the spectrum.However, this depended on the offline training results of the samples.In addition to software methods, optical path improvements can be utilised to eliminate the Doppler background.The use of additional hardware may increase the system's complexity and cost.
This study introduces a new SA peak detection algorithm utilised for laser locking.To circumvent the additional phase shift induced by traditional low-pass filters, a Gaussian filter based on the Gaussian continuous wavelet transform (GCWT) was devised to smooth the noise and jitter of the spectrum.Moreover, a first-order Gaussian continuous wavelet transform (FGCWT) with an adaptive scale factor selection method is proposed to suppress the interference of the Doppler background to determine the peak boundaries.Furthermore, this technique transforms the complex recognition of weak multioverlapping peaks into simple localisation of strong single peaks.Finally, the SA peak position was determined by fitting the SA Lorentzian contour features, which enabled the identification of the actual resonance frequency produced by the atomic hyperfine-level transition rather than the maximum peak in the test data.

The SAS peaks-seeking method
To facilitate the analysis, the spectral data contained in a sweep period can be represented by the following expression: where y sas (t) represents the saturation spectral data within the sampling period.x peak−i (t) is the ith SA peak.f Doppler (t) is the The intensities and linewidths of the SA peaks were smaller than those of the Doppler background.Furthermore, peak recognition was achieved by the following three key steps, which are based on the time-frequency properties of the SA signal components.First, spectrum denoising was employed to mitigate the impact of noise.Second, peak boundary segmentation was conducted to differentiate the SA peaks from other peaks.Finally, SA peak modelling and localisation were performed to identify peak positions.

Spectrum de-noising with GCWT
To eliminate the interference of noise on peak recognition, the time-frequency characteristics of the spectrum were analysed using the short-time Fourier transform method.Figure 1 shows the time-frequency curves of the 780 nm and 795 nm SA spectra, indicating that the frequency range of the SA peaks is confined to within 100 Hz, and that the noise in the spectrum is composed of Gaussian noise.
The traditional low-pass filter always introduces an additional frequency shift owing to phase delay characteristics when filtering out noise, which must be avoided in highprecision peak recognition.Considering the advantages of zero phase shift of the Gaussian function, the GCWT can be formulated in equation ( 2): where y sas (t) denotes the input signal.g 1 (t) = e − t 2 2 .a and τ denote the dilation parameter and location parameter, called scale factor and the time translation, respectively.Equation ( 2) is the sum of y sas (t) after adding Gaussian windows at different positions.In this study, because there is no mutation point in the SA spectrum, the uniform smoothing weight is selected with the time transfer parameter τ = 1. Figure 2 illustrates the relationships among the scale factor (a), the SAS signal, and the transformed result GCWT.As the scale factor increases uniformly in the range of 0-20, the roughness of the spectrum gradually reduces from smoothing; however, when the scale factor increases, the weak peaks in the spectrum are blurred.
Therefore, the signal-to-noise ratio (SNR) and root mean squared error (RMSE) were calculated between the wavelets at different scales to select the optimal scale factor.SNR and MSE are the most used criteria for assessing noise reduction performance and are presented in figure 3.They are defined as where sA spectrum -denoise is the sum of the amplitudes of the spectrum after GCWT and sA spectrum is the sum of the amplitudes of the sampled original spectrum.A spectrumi is the sampled ith original spectrum data and A spectrum -denoisei is the calculated i-th spectrum data with GCWT.A larger SNR and a smaller MSE indicate a better noise reduction effect.As the scale factor increased, the AC amplitude of the low-frequency signal increased; however, the high-frequency detail was smoothed to a smaller value because of the reduction in the Gaussian window.The SNR and MSE results are shown in figure 3. The SNR initially increased and decreased with the scale-factor increment, while the MSE showed the opposite trend.The optimal scale factor yielded the highest SNR and smallest MSE, which is consistent with figure 2. In addition, to illustrate the phase-delay-free characteristics of the GCWT, figure 4 shows the noise reduction effect of the original spectrum using this method and a traditional fourth-order FIR lowpass filter.Moreover, the GCWT denoising method reduces phase delay and improves peak-finding accuracy.

Peaks boundary segmentation
In addition to the interference of white noise, the concaveconvex Doppler background formed by the atomic Doppler effect interferes with the accurate identification of absorption peaks.The mechanism of the SA peaks is the susceptibility absorption of hyperfine transitions; however, the Doppler background is the Doppler broadening of the atoms' thermal motion.Therefore, the linewidth of the peaks was similar to that of the Lorentzian peak, while the Doppler peak had a broad Gaussian profile [20].Therefore, when conducting differential operations on the denoised spectrum, the curvature of the SAS peak was greater than that of the Doppler background.This outstanding feature distinguishes between the Doppler background and the SA peaks.Direct differential operations can result in weaker signals and additional spectral jitter.Based on their differential properties, Where W g (x (a, τ )) is the GCWT mentioned in the previous section, and W g 1 (x (a, τ )) is the FGCWT.The FGCWT can smooth the Doppler background noise without introducing high-frequency noise.
The SA peak position corresponds to the zero position after the FGCWT.By identifying extreme values, the SA peaks in the spectrum can be divided into several single peaks.Figure 5 shows a matrix diagram of the SA spectrum after FGCWT.Scanning scale factor a, the smaller scale factor in the FGCWT, facilitates the identification of weak peaks, but it is susceptible to false peaks arising from jitter or noise.When the scale factor was increased, the identification of the overlapping peaks became blurred, and the Doppler curvature was amplified.The peak boundary blurring phenomenon, evident in the overlapping SA peaks in the 780 nm spectrum, increases with an increase in the wavelet scale factor.
As shown in figure 6, when a is small, the 780 nm spectrum has six extreme values of the FGCWT, corresponding to the boundary of the three SA peaks.However, because of the blurring effect of the overlapping peaks, only three extreme values could be distinguished when a was equal to 20, and the boundary of the three overlapping SA peaks could no longer be accurately identified.
Consequently, an adaptive threshold is employed to search for the optimal scale factor of the FGCWT wavelet matrix.The optimisation steps are as follows: A. The scale factors, a, are scanned from small to large, and the amplitude threshold and width threshold at the current scale factors are computed as follows: where mean (FGCWT a=1 ) is the average value of the FGCWT results at a = 1 and std (FGCWT a=1 ) is the standard deviation of the FGCWT results at a = 1.peak width is a fixed value used to determine the optical path structure.B. According to the mountain prominence concept, whether it is noise or a peak, there are adjacent valleys and peaks.The distance from the valley bottom to the peak was defined as the prominence.After FGCWT, the prominence of the noise was smaller than that of the extreme value of the absorption peak after transformation.The valleys and peaks of each peak were searched using a heap-sort algorithm.Figure 7 shows the original spectrum and the FGCWT spectrum after the previous steps.The red and blue dotted lines mark the peak and the valley bottom of the peak, respectively.C. The prominence of each peak was compared with the amplitude threshold.The distance between two adjacent peaks is defined as the width and is calculated for comparison with the width threshold.The two extreme widths of the FCGWT peaks were determined based on the widths of the absorption peaks.The jitter interference occurs when the peak prominence is below the amplitude threshold.Conversely, it is an effective extreme value for the FGCWT.Subsequently, adjacent peak distances and width thresholds were compared.According to the previous analysis, when the scale factor is increased to gradually blur the signal, the adjacent peak width of the extreme value increases.If the calculated width is smaller than the width threshold, the denoised spectrum is transformed using FGCWT after increasing the scale factor.The prominence and width were calculated and compared repeatedly until a threshold was reached.D. The positions of adjacent peaks whose spacing satisfies the threshold condition are marked as a group of SA peak boundaries.The zero position between the boundaries was the initial position of the SA peak.E. The number of zeros in the extreme value boundary was counted, and it was determined whether it matched the number of true peaks in the SA spectrum.If the number of zeros exceeds the number of actual peaks, the scale factor is increased to execute the above steps.

Peaks localization with a fitting algorithm
In general, the mathematical model of the SA peak is a Lorentz function that does not consider bias and includes parameters such as peak position, width, and amplitude.However, the asymmetry caused by the bias must be considered when establishing a mathematical model of the SA peak superimposed on a Doppler background.This asymmetry can result in a peak shift.
Consequently, with the information on the extreme and zero points obtained in section 2.2, the new model of the SA peak is built as follows: where t i denotes the peak position corresponding to the actual resonance frequency of the transition.γ i is the peak's full width at half maximum (FWHM).where k i is the ith peak intensity, a i is the Doppler slope, and b i is the Doppler bias.
In conjunction with the initial values obtained in section 2.2, the relationship between the peak parameters in equation ( 4) was established.The maximum values on the left and right can be expressed as follows: Subsequently, according to the peak amplitudes of the left and right boundary positions equal to the FWHM of the peak, equation ( 6) is obtained: The magnitude of the extremum as equations ( 7) and ( 8) Finally, the SA peak model parameters were computed by solving these four equations, the residual between the predicted model results and the test data was calculated based on equation ( 9), and the predicted results were obtained using the Levenberg-Marquardt algorithm [21].The iteration stops when the residual is less than 10 −10 .

Off-line simulations
We performed offline simulations of the collected 780 nm and 795 nm SA spectra.Figure 8 shows the convergence efficiency of the LM fitting for the 780 nm spectrum and 795 nm spectrum.All the test data for the SA peaks were accurately fitted with the predicted data.The identified absorption peaks are the reference frequency points commonly used in frequency stabilisation experiments.Despite having several weak peaks on the right side of the f peak in the 780 nm spectrum, owing to the weak amplitude and poor symmetry of these peaks, the dynamic range of the demodulated error signal is too narrow to be locked stably.Therefore, they were regarded as false peaks.The weak peak between the d and e peaks in the spectrum of the 795 nm spectrum is also a false peak for the same reason.In the process of identifying the boundary of the algorithm, after the FGCWT of the false peak, its prominence is much smaller than that of the normal peak owing to asymmetry, and the peak spacing of the adjacent peaks is greater than the width threshold; therefore, it can also be automatically determined as a false peak.
To further evaluate the effectiveness of the algorithm, it was used to detect the laser mode-hopping point.The modehopping and peak points are marked and enlarged in Boxes as shown in figure 9.As the mode-hopping point satisfied the amplitude threshold but not the width threshold.Therefore, the mode-hopping point was accurately identified as a false peak.
This study compared the fitting effects of an unbiased Gaussian function model, an unbiased Lorentz function model, and a biased Lorentz model through simulations, as shown in figure 10.The characteristic parameters of the three models were determined based on the initial values obtained  in section 2.2.The prediction methods employ the same LM algorithm that can adaptively adjust the parameters to converge rapidly and accurately.It is shown than the fitting accuracy of the biased Lorentz model was greater than that of the unbiased Lorentz model, and that of the unbiased Lorentz model was greater than that of the unbiased Gaussian model.

Experiment setup and results
The experimental setup for frequency stabilisation and measurement is shown in figure 11(a).After the isolator, the laser was divided into two beams using a splitter plate (50%:50%).One of the reflected light beams was incident on the Rb atomic  vacuum cell, which was used to optically pumped rubidium atoms for achieving hyperfine level transitions.And the other transmitted light beam was incident on the FP etalon (Thorlabs SA210) used to calibrate the laser frequency, which is a commonly used method for precise calibration of laser frequency [22].The light pumped through the vacuum cell was divided into two beams by a high-transmittance, low-reflection mirror (75%:25%).The transmitted light was reflected by a plane mirror, returned to the cell and received by the photodetector.Another light beam was incident on a wavelength meter (WS7-30) with 1 MHz resolution and 30 MHz absolute accuracy.
First, a commercial analogue frequency-locking controller was used for peak finding and frequency locking.The locking results are presented in figure 11(c).Subsequently, a selfmade digital frequency-locking controller with an automatic peak-finding algorithm was verified.Triangular sweep and SA spectral signals were collected in the FPGA using a 16-bit analogue-to-digital converter.The spectral signal was saved and transmitted to the RAM of the DSP processor by testing the initial trigger of the period-triangular sweeping signal.After the saved spectral signal ran the peak-seeking algorithm, the triangular wave signal corresponding to the peak-seeking position was output to the laser current source as an accurate tuning value.
The locking results are presented in figures 11(b) and (c).The locking process of automatic peak finding method is much quicker than the analogue controller.As the analogue controller requires continuous adjustment of the laser current to approach the frequency transition.The frequency fluctuation after automatic peak seeking and frequency locking was 3.7 MHz h −1 , and that after manual peak seeking and frequency-locking was 4.7 MHz.Despite the frequency stability being affected by many factors, the frequency-locking process of automatic peak finding is more agile and has better repeatability than the manually tuned peak finding.

Conclusion
In this study, we proposed a novel automatic SA peak detection method.The CGWT and FGCWT are employed to eliminate noise, and the Doppler background is superimposed on the spectrum, which reduces the misjudgment caused by false peaks.Simultaneously, an intelligent fitting algorithm was studied to further improve the accuracy of the peak-finding algorithm.This study demonstrates that a new method combining CGWT and FGCWT with Doppler background superimposition significantly improves the accuracy of identifying and isolating SA peaks.This method is applicable not only to SA stabilization, but also to PDH stabilization and other automatic frequency stabilization systems that require automatic identification of frequency references.Future research must focus on reducing computational complexity, increasing sampling rate, and practical applications in signal processing used for intelligent loss-of -lock judgment.

G Li et al
accessible or reusable by other researchers.The data that support the of this study are available upon reasonable request from the authors.

Figure 1 .
Figure 1.The time-frequency diagrams of the 780 nm and 795 nm SAS spectra.Figures (a) and (c) are the time domain diagrams of the 780 nm SAS signal and 795 nm SAS signal.Figures (b) and (d) are their time-frequency diagrams obtained by short-time Fourier.

Figure 2 .
Figure 2. The three-dimensional graph representing the relationship of the spectrum, wavelet scale factor, and the GCWT.(a) The three-dimensional graph of 780 nm spectroscopy.(b) The three-dimensional graph of 795 nm spectroscopy.

Figure 3 .
Figure 3.The relationship between spectral Gaussian noise reduction performance and scale factor.The blue line represents the signal-noise-ratio varying with the scale factor.The red line represents the root mean squared error varying with the scale factor.

Figure 4 .
Figure 4.The frequency shift introduced during the denoising process.(a) The zero-frequency shift of the GCWT and the large frequency shift of the FIR low pass filter on the 780 nm spectrum.(b) The zero-frequency shift of the GCWT and the large frequency shift of the FIR low pass filter on the 795 nm spectrum.

Figure 5 .
Figure 5.The FGCWT matrix of the spectrum with different scale factors.(a) The FGCWT matrix diagram of the 780 nm spectrum.(b) The FGCWT matrix diagram of the 795 nm spectrum.

Figure 6 .
Figure 6.The FGCWT under different scale factors.(a) The FGCWT spectrum under different scale factors of 780 nm spectrum.(b) The FGCWT spectrum under different scale factors of 795 nm spectrum.

Figure 7 .
Figure 7.The prominence of the FGCWT spectrum.The red line is the part of the SA peaks.The blue line is the FGCWT of the SA peaks.The pentagram and the circle respectively mark the valley and peak.The black cross marks the false peak that does not satisfy the threshold condition.

Figure 8 .
Figure 8. Nonlinear fitting results of the SAS peaks.(a) Nonlinear fitting results of 795nm SAS peaks and the enlarged fitting result of e peak.(b) Nonlinear fitting results of 780nm SAS peaks and the enlarged fitting result of f peak.

Figure 9 .
Figure 9. Peak identification of SA spectra with mode-hopping points.The red line is the SAS curve with a mode-hopping point, and the blue line is its corresponding FGCWT curve.

Figure 10 .
Figure 10.Fitting results with different

Figure 11 .
Figure 11.(a) Experimental setup of the saturated absorption frequency stablization.(b) The frequency locking result with the manual peak tuning and locking algorithm.(c) The frequency locking result with the automatic peak finding algorithm.