Thanks to its intrinsic advantages of high simplicity, low cost, and the same resolution as that of the typical Michelson interferometer, self-mixing interferometry (SMI), also named optical feedback interferometry (OFI), has been well-established in various metrology applications, and not only in traditional industrial fields, such as vibration displacement, velocimetry, and absolute distance finding^[1–5], but also in laser welding monitoring, particle sizing, or label-free biomedical sensing^[6–8].

Currently, SMI has attracted much attention in microparticle detection in chemical and biomedical fields, since this simple technique is very suitable for an integrated platform with microfluidic chip and other devices like fuel nozzles^[9–14]. Particularly, SMI can be a promising candidate for real-time single-particle detection for small dosage agent experiments in flow cytometry or particle mixing. Moreira et al., for the first time, observed the presence of single submicroscale particles in a 320 µm channel using SMI technology. They applied a simple bandpass filter to visualize the particle bursts in the noisy raw signal and measured the single-particle velocity by fast Fourier transform (FFT)^[15]. Contreras et al. managed to classify the single polystyrene sphere size for a higher signal-to-noise ratio (SNR). They developed a novel SMI scheme with edge-filter-enhanced self-mixing interferometry (ESMI)^[16], and they approached around 2 orders of magnitude, even at a 10 m operating distance. Zhao et al. presented a fringe counting method for single-microparticle sizing based on the Hilbert transform (HT)^[17]; the resolution was still a half wavelength. Unlike the continuous SMI signal from a bulk translating target, particle-induced signals represent a discrete waveform burst. This innovative topic is still challenging: the SMI signal is lower in amplitude and temporally discrete on a millisecond scale. The noises from the measurement environment can significantly reduce SNR and spoil the performance of the sensor system. The resolution in this specific scenario is limited to a basic half wavelength; what is more, the retrieved fringe pattern still presents the signal ambiguity in the waveform edge, where some fringes are missing.

Besides measuring temporal signal waveform fringe number^[18] and frequency spectral width^[7], novel precise alternatives for particle detection are still required. As a kind of phase observation tool, the phase-unwrapping method (PUM) has been exhaustively studied in SMI sensor systems^[19–27], and the spatial resolution even can reach $\lambda /67$^[28]. This processing method normally involves two steps: the first step is phase retrieval, and the second step is to estimate the feedback factor and linewidth enhancement factor for phase profile correction and displacement reconstruction. Apart from the typical vibration or displacement measurements, SMI signal phase assessment was also employed in 2D scanning and surface profile mapping. Lacot et al. proposed a new low-noise and phase-sensitive detection scheme for scattering-type scanning near-field optical microscopy (sSNOM). They successfully implemented high-resolution phase imaging of silicon-on-insulator (SOI) optical waveguides^[29]. Considering the high resolution, PUM can be a good candidate for particle detection; however, there appears to be no related literature existing on this topic.

The paper is organized as follows. First, a theoretical framework of SMI effect in the presence of a translating particle is demonstrated, and during the particle passage, the dependence of SMI signal phase variation on particle size is investigated. Second, for effective noise elimination and phase retrieval, continuous wavelet transform (CWT) is directly applied to the original SMI signal. And then a simple HT-based PUM is employed for phase calculation and unwrapping, with the resultant phase variation corresponding to single particle extraction. Finally, using a 532 nm solid-state laser SMI system, signal bursts induced by different-sized polymer beads are acquired and processed by our algorithm, and the measured phase variation values are compared with the simulation results.

When a laser shoots a beam onto a remote target and a portion of the backreflected or backscattered light re-enters the laser cavity, the interaction between the coupled feedback light and the emitting laser field evokes a laser output power fluctuation. The modulated laser output power $P(t)$ in the presence of self-mixing effect can be expressed by Eq. (1), $$P(t)=P_0[1+m(t)\cos \varphi ].$$(1)

$P_0$ is the original laser output power, and $m$ is the modulation index, which depends on the scattering intensity. $\varphi$ denotes the phase in the presence of the SMI effect, whose value can be calculated by the well-known excess phase equation^[30], $${\varphi }_0=\varphi +C(\cos \varphi +\arctan \alpha ).$$(2)

${\varphi }_0$ is the initial phase without feedback effect, $C$ is the feedback coefficient, determining the feedback strength, and $\alpha$ is the linewidth enhancement factor^[31]. But when the target is a microscale moving particle, the backward scattering light acts as the feedback light, and the $C$ value is typically smaller than 1^[11]. The SMI system then operates in a weak feedback region, and the phase hysteresis can be negligible. Due to the well-known Doppler effect, the instant phase $\varphi$ value can be approximately given by^[32,33]$$\varphi (t)=2\pi f_d·t.$$(3)$f_d$ is the Doppler frequency that depends on the particle velocity and the incidence angle between the laser axis and the particle translation direction ${\theta }_{\mathrm{inc}}$, $$f_d=\frac{2V·\cos {\theta }_{\mathrm{inc}}}{\lambda }.$$(4)

Considering that the incident laser beam exhibits a Gaussian spatial intensity profile, and the scattering intensity properties of the particle are also strongly dependent on its position during the flight, both phenomena influence the modulation index $m$ value. Thus, we assume the $m$ value varies as a Gaussian function during the particle flowing inside the measurement volume, and the signal expression can be rewritten by Eq. (5)^[17], $$P(t)=P_0(1+m·\cos f_d·t)\exp [-\frac{{(t-t_0)}^2}{2{\tau }^2}].$$(5)

$\tau$ is the particle passage period inside the measurable region of the sensor system, whose dimension $M$ can be simplified as the combination of laser spot size $K$ and the particle diameter $D$. When $t=t_0$, the laser beam shoots onto the particle or cell center perpendicularly, and the particle center passes the light axis^[18], $$\tau =\frac{M}{V\sin \theta }=\frac{K+D}{V\sin \theta }.$$(6)

During the signal burst period $\tau$, the total phase variation $\mathrm{\Phi }$ is given by Eq. (7), as a linear function of $\tau$, without dependence on velocity, $$\mathrm{\Phi }=2\pi f_d·\tau .$$(7)

From Eqs. (6) and (7), the phase variation $\mathrm{\Phi }$ can be derived as a proportional function of the detectable region distance $M$ by Eq. (8). $K$ and $D$ are fixed in the given experiment; neither of them depends on the velocity, $$\mathrm{\Phi }=\frac{4\pi }{\lambda }·\frac{\cos \theta }{\sin \theta }·M=\frac{4\pi }{\lambda }·\frac{\cos \theta }{\sin \theta }(K+D).$$(8)

The particle diameter $D$ can be very easily calculated by the given equation, $$D=\frac{\lambda }{4\pi }·\frac{\sin \theta }{\cos \theta }\mathrm{\Phi }-K.$$(9)

We simulate the signal burst and the phase variation with the parameters in Table 1. The simulation results are shown in Fig. 1. Within the signal burst period (dashed line), the phase increases linearly as Eq. (7). The signal burst region can be exactly distinguished by the phase profile.

According to the above-mentioned section, the particle-induced SMI signal is in a temporally intermittent form during the particle passage. Moreover, the signal burst normally is affected by white Gaussian noise, impulsive noise, and DC pedestal^[34,35]. The signal burst of each particle is ambiguous and merged in the noises, so an efficient method of burst identification and denoising is demanded. The CWT is a well-established method for investigation of time and frequency details of signals whose frequency content varies over time^[23,36,37]. Eventually, this technique performs a correlation between a scaled/shifted wavelet basis and the SMI signal. Because the resultant amplitude corresponding to the uncorrelated and random noise is much lower than that of the desired signal, CWT can implement noise removal in a single step without additive complicated processing.

In this paper, CWT is applied to the SMI signal $\mathrm{Sig}(t)$. Based on CWT principles^[21], the CWT coefficient $W$ of the signal is expressed as follows: $$W(u,s)=\frac{1}{\sqrt{s}}{\int }_{-\infty }^{+\infty }\text{Sig}(t)\psi \left(\frac{t-u}{s}\right)\mathrm{d}t,$$(10)where $s$ is the scaling factor, describing the time length of single-period wavelet basis, and $u$ is the shifting factor. By changing the $t$ and $u$ values, CWT enables effective SMI signal decomposition into different frequency components.

$\psi$ is the wavelet basis. There are many types of wavelet basis, such as Daubechies (db), Morlet, Coiflets, and Haarhere. Considering the particle-induced frequency modulation signal is sinusoidal-like in Eq. (5), we chose Morlet, whose capability has been proven in SMI measurements^[38], as the wavelet basis in the data processing. The Morlet CWT basis with a frequency ${\omega }_0$ in the time domain is defined as follows: $$\psi (t)={\pi }^{-1/4}{e}^{i{\omega }_{0}t}{e}^{-t^2/2}.$$(11)

The HT can represent a signal in its analytical form over the orthogonal plane and retrieve the phase information very conveniently without modifying the amplitude. So HT is very suitable for phase unwrapping in SMI measurements^[17,19].

The original signal function $f(t)$ can be transformed into a complex analytical form as a combination of a real part $\xi (t)$ and an imaginary part $\mathrm{\Theta }(t)$, as in Eq. (12), $$f(t)=\xi (t)+j·\mathrm{\Theta }(t),$$(12)and $\mathrm{\Theta }(t)$ is calculated by performing a 90° phase shift of the real part by an HT^[17], $$\mathrm{\Theta }(t)=\mathrm{HT}[\xi (t)].$$(13)

Finally, the instant phase $\varphi (t)$ can be calculated using inverse tangent function, $$\varphi (t)=\arctan \frac{\xi (t)}{\mathrm{\Theta }(t)}.$$(14)

After the HT, using a MATLAB-customized phase-unwrapping routine, the realistic phase $\mathrm{\Phi }$ can be calculated after the unwrapping, $$\mathrm{\Phi }=\varphi (t)+2n\pi ,$$(15)where $n$ is an integer, and the value of $n$ is 0, 1, 2, ….

The data processing flow chart is illustrated in Fig. 2. The algorithm involves two steps: first, CWT for raw signal noise elimination and burst retrieval; second, HT for phase estimation and unwrapping function.

To validate the simulation and the capability of our algorithm, a series of particle sizing measurements are performed with the setup shown in Fig. 3. We select a compact diode-pumped solid-state laser (DPSS) at 532 nm wavelength (Thorlabs, DJ532-10) as the SMI laser source in the system. Due to the higher ratio of fluorescence-to-photon lifetime, a DPSS SMI configuration has higher sensitivity than the typical laser diode (LD) setup^[39]. The laser is controlled by a laser controller (Thorlabs, ITC4001), and the injection current and the laser package temperature are kept to be 240 mA and $24°\mathrm{C}\pm 0.1°\mathrm{C}$, respectively. The collimated output light is spliced into two beams by a 1:9 beam splitter (BS). One beam path (90%) as a detection path is focalized tightly by an aspherical lens (Thorlabs, CMD240-C) and shot onto the microparticle flux. The laser spot diameter $K$ is around 1 µm, and the incidence angle ${\theta }_{\mathrm{inc}}$ is $25\pm 1\mathrm{deg}$. The other beam path (10%) is directed to a fixed gain photodiode detector (Thorlabs, PDF10A2) for SMI signal retrieving. Afterward, the SMI voltage signal is digitized using a fast speed data acquisition card (NI, 6361USB). The sampling frequency and acquisition data number are 500 kHz and $2^{17}$, respectively. Polystyrene spheres (PSs) in seven different diameters (100, 500, 2, 3, 4, 5, and 6 µm) are employed as the scattering objects, and the particle stream is made by a homemade powder feeder. The dimension coefficients of variation (CVs) of the particles are better than 3%.

Thousands of monodispersed signal bursts in different particle sizes are captured and processed by our above-mentioned algorithm. The raw signal burst without processing in 2 µm is denoted by the black line in Fig. 4(a). As the figure shows, the signal burst is merged in the strong noise floor, leading to great difficulty in interference fringe observation. To reduce the noise, we applied both bandpass filter (BPF) and Morlet-based CWT upon the raw signal, respectively. Firstly, we analyze the signal by CWT, and the scalogram of the real part of the CWT coefficient $W$ is then depicted in Fig. 4(b) to visualize both the spatial and temporal components. It can be seen that the maximum $W$ value in light yellow color is located in the temporal range from 21.3 to 23.5 ms, representing the particle flight period inside the laser illumination area, and the frequency is in a limited range from 2 to 4 kHz.

We set the threshold to be 0.2 of the maximum $W$ value and eliminate the signal level under the threshold value. Besides the CWT, we also use a fifth-order Butterworth BPF as reference^[17] for signal denoising. We set the filter cut-off frequency range to be from 2 to 4 kHz, the same as the scalogram shows. The blue line and red line denote the denoised bursts from the CWT and the BPF in Fig. 4(a), respectively. The bursts are added $\pm 0.4\mathrm{V}$ amplitude offset for more convenient observation.

From the figures, it can be seen that even though the BPF can eliminate the noise effectively in the red line, the burst margin regions are still greatly ambiguous. There are some supernumerary noise fluctuations remaining in the dashed line circles, which can be mistaken for extra signal interference fringes. Compared with the BPF, a well-defined impulse waveform corresponding to a single particle can be distinguished easily by CWT (blue line) without margin ambiguity.

Afterward, the denoised signals are processed by the HT-PUM algorithm to calculate the phase profile. As shown in Fig. 5, obviously, both phase profiles present good linear trend as the simulation. However, compared with the one from the BPF, the linear profile CWT bends at the burst period endpoint, which is exactly the same as in the scalogram in Fig. 4(b). The total phase variation from the starting point $t=21.3\mathrm{ms}$ to the ending point $t=23.5\mathrm{ms}$ is around 31 rad, which is in good agreement with the calculation from Eq. (8). From the results, we can conclude that even without the scalogram prediction, CWT can precisely extract the signal scope in both time domain and the phase profile regardless of the noise.

By using the threshold CWT algorithm, the particle signals in other different particle sizes are also illustrated in Fig. 6, indicating the robustness of denoising ability with CWT.

Signal burst retrievals from each particle size are repeated around 30 times; the mean values of the measured phase $\mathrm{\Phi }$ as a function of the particle diameter (blue marks) and the simulation by Eq. (8) (black line) are depicted in Fig. 7. It should be noted that the measurement results present a good linear profile with respect to the linear fitting (red line), which is expected in simulation. The measurement and simulation results are generally in good agreement over the whole given particle size range, with a slight deviation. The measured laser spot size $K$ value is a little smaller than the simulated one. This can be explained by the uncertain quantities of the measurement value before the particle sizing experiment. The slope of the measurement trend (the red line) is greater than the simulated one; the incident angle deviation is the most credible reason for this phenomenon. Unfortunately, both $K$ and ${\theta }_{\mathrm{inc}}$ discrepancies are difficult to avoid in practice.

In this work, we present a simple and capable PUM algorithm for single-particle size characterization. CWT and HT are involved for noise elimination, individual particle-induced burst retrieval, and phase calculation. The resultant phase variation corresponding to single-particle modulation is extracted efficiently compared with the typical bandpass filtering. The experimental results showed good linear trend of SMI signal phase variation with the particle size over a practical valuable measurement range from 100 nm to 6 µm, which corresponds to typical living cells and bacteria. The good consistency between the simulation results and measured signals has proven the reliability of our theoretical model and the capability of phase observation method in microparticle detection. Thus, the advantage of the new approach with our CWT-HT-PUM algorithm can be more pronounced for particle sizing in the biomedical and chemical fields. Unlike in displacement or vibration measurements, a piezoelectric translator (PZT) enables ultraprecise target motion on a 0.1 nm scale. In our measurement, the artificial particle diameter cannot be regulated less than 1 µm, and the particle dimension inhomogeneity is inevitable during fabrication. So the minimum particle sizing limitation of our method cannot be validated perfectly for the moment. More experiments with finer particle size deviation could be performed in future.