SS-OCT Resampling Method and System Based on Undersampling Clock Recovery

Abstract

This invention discloses an SS-OCT resampling method and system based on undersampled clock recovery, belonging to the field of signal processing technology. It synchronously acquires a k-clock signal and an OCT interferometer signal under undersampled conditions; performs a Hilbert transform on the k-clock signal to obtain an analytic signal, from which an aliased phase sequence is extracted; performs a first-order difference operation on the aliased phase sequence; constructs several candidate frequencies for each sampling moment; and uses a dynamic programming algorithm to select the optimal frequency from the candidate frequencies to obtain the recovered instantaneous frequency sequence. The continuous phase sequence of the k-clock signal is reconstructed through cumulative integration; and the continuous phase sequence is used to perform equal-phase-interval resampling on the OCT interferometer signal to obtain an equal-wavenumber-interval OCT resampled signal. This method can accurately recover the instantaneous frequency of the k-clock signal under undersampled conditions, completing the resampling process and thus reducing system hardware costs.

Description

Technical Field

[0001] This invention relates to the field of signal processing technology, and more specifically to an SS-OCT resampling method and system based on undersampled clock recovery. Background Technology

[0002] During operation, a swept-frequency optical coherence tomography (OCT) system requires the simultaneous recording of two signals. The first signal is the OCT interferometric signal, which carries depth information of the sample under test. The second signal is the k-clock signal, generated by the Mach-Zehnder interferometer inside the swept-frequency laser, whose instantaneous frequency is linearly related to the wavenumber. Resampling the OCT signal at equal wavenumber intervals using the zero-crossing points or phase information of the k-clock signal is a necessary step to eliminate swept-frequency nonlinearity and ensure axial resolution.

[0003] Existing resampling schemes fall into two categories. One is the direct hardware sampling clock drive method, where the k-clock signal output by the swept laser is directly connected to the external clock input pin of the analog-to-digital converter (ADC), and each k-clock pulse triggers an OCT signal acquisition. However, this method requires strict matching of clock frequency, amplitude, jitter, and duty cycle (e.g., 300MHz to 800MHz). The absence of the sweep interval clock can easily lead to ADC instability. Furthermore, impedance matching, signal integrity, and delay calibration requirements are stringent, and the hardware interface needs to be redesigned for different laser parameters, resulting in poor adaptability.

[0004] The second method is software resampling that satisfies the Nyquist sampling theorem. This approach acquires the k-clock signal at a sampling rate no less than twice the highest k-clock frequency, extracts the instantaneous phase using Hilbert transform, and then matches it with a preset ideal equal-wavenumber interval phase sequence. Resampling is then completed using SINC interpolation. When the highest k-clock frequency reaches 800MHz, a sampling rate of at least 1.6GSPS is required. This requirement significantly increases the performance demands on the analog-to-digital converter, leading to a substantial increase in the cost of the acquisition card.

[0005] Under the aforementioned limitations, when the system uses a conventional 1GSPS analog-to-digital converter to acquire k-clock signals with a maximum frequency of 800MHz, if the sampling rate (1GSPS) is less than twice the highest k-clock frequency (800MHz) (i.e., 1GSPS < 1.6GSPS), the Nyquist sampling theorem cannot be satisfied, and severe frequency aliasing will inevitably occur in the k-clock signal. Aliasing causes the instantaneous phase directly extracted from the acquired data to lose its uniqueness; the aliased phase at each sampling point may correspond to two different original frequencies, making it impossible to directly determine the true instantaneous frequency. Without an accurate instantaneous frequency, subsequent equal-wavenumber interval resampling based on the k-clock phase cannot be performed, and the resolution and signal-to-noise ratio of the OCT image will be severely degraded.

[0006] Therefore, how to uniquely recover the true instantaneous frequency from the aliased k-clock signal under under-sampling conditions and achieve high-precision resampling has become a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0007] The purpose of this invention is to provide an SS-OCT resampling method and system based on undersampled clock recovery, which can accurately recover the instantaneous frequency of the k-clock signal during undersampling and complete the equal wavenumber interval resampling of the OCT signal, thereby reducing the system hardware cost.

[0008] To achieve the above objectives, the following plan is adopted: On the one hand, the present invention provides an SS-OCT resampling method based on undersampled clock recovery, specifically including the following steps: S1. Simultaneously acquire the k-clock signal and the OCT interference signal using an analog-to-digital converter with a first sampling rate, wherein the first sampling rate is the undersampling frequency of the k-clock signal; S2. Perform Hilbert transform on the k-clock signal to construct an analytic signal and extract the aliased phase sequence from the analytic signal. Then, perform a first-order difference operation on the aliased phase sequence to obtain the aliased instantaneous frequency sequence. S3. Based on the Nyquist folding relation, several candidate frequencies are constructed for each sampling time using the aliased instantaneous frequency sequence to form a candidate sequence. The optimal frequency path is selected from the candidate sequence using a dynamic programming algorithm to obtain the recovered instantaneous frequency sequence. S4. Using the recovered instantaneous frequency sequence, the continuous phase sequence of the k-clock signal is reconstructed by cumulative integration. The continuous phase sequence is then used to perform equal-phase-interval resampling on the OCT interference signal to obtain an OCT resampled signal with equal wavenumber intervals.

[0009] In some specific implementation schemes, the specific process for obtaining the instantaneous frequency sequence of aliasing is as follows: S21. Extract the wrapped phase from the analytical signal, and then expand the wrapped phase to obtain a continuously expanded aliased phase sequence. S22. Perform a first-order difference operation on the aliased phase sequence, and multiply the result of the first-order difference operation by the first sampling rate to obtain the aliasing frequency. S23. Take the absolute value of the aliasing frequency to obtain the absolute value frequency of aliasing, and finally output the instantaneous frequency sequence of aliasing.

[0010] In some specific implementations, the candidate frequencies include two, namely the first candidate frequency and the second candidate frequency, denoted as: c1[i] = f_abs[i]; c2[i] = f_s – f_abs[i]; Where c1[i] represents the first candidate frequency at the i-th sampling time, c2[i] represents the second candidate frequency at the i-th sampling time, f_abs[i] represents the absolute value frequency of aliasing at the i-th sampling time in the aliasing instantaneous frequency sequence, and f_s represents the first sampling rate; Simultaneously, the two candidate frequencies are filtered for validity using the known frequency range of the k-clock signal, and the valid candidate frequencies are output.

[0011] In some specific implementation schemes, the specific process for obtaining the recovered instantaneous frequency sequence is as follows: S31. Set the number of state spaces in the dynamic programming algorithm to 2. Define the two candidate frequencies at each sampling time as the two states in the state space. Set a cost function to describe the cost of state transition between two adjacent sampling times. The cost function includes three components: frequency jump cost, direction penalty term, and acceleration penalty term. S32. Calculate the cost after each state transition based on the two states and the cost function, and calculate the cumulative cost at all sampling times using the recursive formula. S33. After completing the recursion of all sampling times, start backtracking from the state with the minimum cumulative cost in the last column to obtain the optimal candidate index for each sampling time, and then obtain the recovered instantaneous frequency sequence.

[0012] In some specific implementation schemes, the construction process of the direction penalty term is as follows: By analyzing the changes in the envelope slope of the instantaneous frequency sequence of aliasing, the boundaries of different stages of the sweep frequency period are estimated. Different directional penalties are applied to the boundaries of different stages, and a transition zone is set at each stage boundary to reduce the weight of the directional penalty within the transition zone.

[0013] In some specific implementations, the boundaries of different stages include rising segments and falling segments. When the candidate frequencies of multiple consecutive sampling times show an overall upward trend, it is defined as the rising segment and a first directional penalty is applied. When the candidate frequencies of multiple consecutive sampling times show an overall downward trend, it is defined as the falling segment and a second directional penalty is applied.

[0014] In some specific implementations, the frequency hopping cost Δf is defined as: Δf = |c_cur[i] – c_prev[i-1]| Where c_cur[i] represents the state at the i-th sampling time, and c_prev[i-1] represents the state at the (i-1)-th sampling time; Set the maximum allowed transition threshold Δf_max. Δf_max is adaptively adjusted according to the boundary of the current stage. If the frequency transition cost Δf exceeds the maximum allowed transition threshold Δf_max corresponding to the current stage boundary, then state transition is prohibited.

[0015] In some specific implementation schemes, the acceleration penalty term is defined as a second-order difference penalty: |(c_cur[i]–c_prev[i-1])–(c_prev[i-1]–c_prev2[i-2])| Where c_cur[i] represents the state at the i-th sampling time, c_prev[i-1] represents the state at the (i-1)-th sampling time, and c_prev2[i-2] represents the state at the (i-2)-th sampling time. When calculating the cumulative cost, the acceleration penalty term is multiplied by the coefficient γ and then included.

[0016] In some specific implementation schemes, step S4 is performed as follows: S41. Set an equal phase interval Δφ, and generate the target phase sequence using a continuous phase sequence; S42. By using linear interpolation, a mapping is established between the continuous phase sequence and the time axis, and the time point t_resampled corresponding to each target phase is calculated. S43. Interpolate the OCT interference signal to extract the amplitude value corresponding to the sampling time at time point t_resampled, and finally obtain the OCT resampled signal with equal wavenumber intervals.

[0017] The concept of this application is as follows: This paper addresses the technical problem in existing swept-frequency optical coherence tomography systems where frequency aliasing of the k-clock signal occurs due to the analog-to-digital converter (ADC) sampling rate being less than twice the highest frequency of the k-clock signal, thus hindering accurate recovery of the instantaneous frequency and resampling. This application directly selects an ADC with a sampling rate less than twice the highest frequency of the k-clock signal for sampling, actively creating frequency aliasing to reduce hardware costs. It constructs a combination of dual candidate frequencies and dynamic programming, building two candidate frequencies based on the absolute value f_abs of the aliasing frequencies, and using these candidate frequencies as the state space of the dynamic programming (fixed to two states). This combination solves the frequency ambiguity problem under undersampling and differs from the traditional multi-state Viterbi algorithm (used for time-frequency plane ridge extraction). Furthermore, a second-order difference penalty is added to the cost function of the dynamic programming to suppress abrupt changes in the frequency change rate, resulting in a smoother recovered frequency curve that conforms to physical inertia.

[0018] Furthermore, to accommodate the segmented nonlinear characteristics of actual frequency sweeping, most existing technologies assume that the k-clock frequency exhibits an ideal linear triangular wave variation. However, actual frequency-sweeping lasers exhibit asymmetric and nonlinear variations with rapid rise, slow rise, slow fall, and rapid fall. This embodiment employs an adaptive direction penalty mechanism, transition region softening treatment, and acceleration penalty mechanism, applying direction penalties of varying intensities at different stages and setting a transition region near the inflection point to reduce the penalty weight. This enables accurate tracking of the actual frequency change trajectory. Simulation results show that the frequency recovery error RMS is less than 100kHz, and the phase error RMS is less than 0.01 radians. Even under a 20dB signal-to-noise ratio condition, the algorithm can still correctly resolve ambiguity, demonstrating superior robustness compared to methods based on the ideal linear assumption.

[0019] Secondly, this application provides an SS-OCT resampling system based on undersampled clock recovery, comprising: SS-OCT signal source, used to generate k-clock signal and OCT interference signal; The acquisition module is used to synchronously acquire a k-clock signal and an OCT interference signal using an analog-to-digital converter with a first sampling rate at an undersampling frequency, wherein the first sampling rate is the undersampling frequency of the k-clock signal; The aliasing frequency extraction module is used to perform Hilbert transform on the k-clock signal, construct an analytical signal, extract the aliasing phase sequence from the analytical signal, and perform first-order difference operation on the aliasing phase sequence to obtain the aliasing instantaneous frequency sequence. The phase recovery module is used to construct several candidate frequencies for each sampling time based on the Nyquist folding relationship and the aliased instantaneous frequency sequence, forming a candidate sequence, and then using a dynamic programming algorithm to select the optimal frequency path from the candidate sequence to obtain the recovered instantaneous frequency sequence. The remapping module is used to reconstruct the continuous phase sequence of the k-clock signal by accumulating the recovered instantaneous frequency sequence, and then use the continuous phase sequence to perform equal-phase-interval resampling on the OCT interference signal to obtain an OCT resampled signal with equal wavenumber intervals.

[0020] The beneficial effects of this invention are as follows: This invention can operate normally even when the sampling rate of the analog-to-digital converter is less than twice the highest frequency of the k-clock (i.e., the k-clock signal is undersampled). It employs a combination of dual candidate frequency construction and dynamic programming. Candidate frequencies are constructed based on the aliased instantaneous frequency sequence, and the candidate frequencies are used as the state space of dynamic programming to recover the true frequency sequence of the aliased signal. Based on the recovered true frequency sequence, continuous phase is reconstructed, and then OCT signal sampling points are directly extracted at equal phase intervals. This solves the frequency ambiguity problem under undersampling. Compared with the traditional full-sampling scheme (which requires more than 1.6 GSPS), it significantly reduces hardware costs and lowers the data transmission bandwidth requirements of the acquisition card. Based on the analysis of the asymmetric and nonlinear changes of actual frequency-sweeping lasers, this invention divides the frequency sweeping process into multiple stages through segmented adaptive directional penalty, transition region softening treatment and acceleration penalty mechanism. Different directional penalties are applied in different stages, and a transition region is set near the turning point to reduce the penalty weight, which can accurately track the actual frequency change trajectory. Attached Figure Description

[0021] Figure 1 The flowchart of the SS-OCT resampling method based on undersampled clock recovery provided in this embodiment of the invention is as follows: Figure 2 This is a block diagram of the SS-OCT resampling system based on undersampling clock recovery provided in an embodiment of the present invention. Detailed Implementation

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the present invention or its application or use. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0023] Unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps described in these embodiments do not limit the scope of the invention.

[0024] Furthermore, for clarity and brevity, descriptions of well-known functions and methods may have been omitted. Those skilled in the art will recognize that various changes and modifications can be made to the examples described herein without departing from the spirit and scope of this disclosure.

[0025] Techniques, methods, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, such techniques, methods, and equipment should be considered part of the specification.

[0026] In all examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0027] Example 1 like Figure 1 As shown, this embodiment provides an SS-OCT resampling method based on undersampled clock recovery, which specifically includes the following steps: S1. Simultaneously acquire the k-clock signal and the OCT interference signal using an analog-to-digital converter with a first sampling rate, wherein the first sampling rate is the undersampling frequency of the k-clock signal; Specifically, the first sampling rate is denoted as f_s, the highest instantaneous frequency of the k-clock signal is denoted as f_max, and the k-clock signal is in an undersampled state when f_max > f_s / 2. For example, the first sampling rate f_s = 1GHz, the highest instantaneous frequency f_max = 800MHz, and the lowest instantaneous frequency f_min = 10MHz.

[0028] S2. Perform Hilbert transform on the k-clock signal to construct an analytic signal and extract the aliased phase sequence from the analytic signal. Then, perform a first-order difference operation on the aliased phase sequence to obtain the aliased instantaneous frequency sequence. The specific process for obtaining the instantaneous frequency sequence of aliasing is as follows: S21. Extract the wrapping phase from the analytical signal, and then expand the wrapping phase to obtain a continuously expanded aliased phase sequence, denoted as φ_unwrap. Due to undersampling, this aliased phase sequence loses high-frequency folding information, and its derivative corresponds to the instantaneous frequency after aliasing. S22. Perform a first-order difference operation on the aliased phase sequence, and multiply the result of the first-order difference operation by the first sampling rate to obtain the aliasing frequency f_alias: f_alias=diff(φ_unwrap)×f_s / (2π); Here, diff represents the first-order difference operation.

[0029] S23. Take the absolute value of the aliasing frequency to obtain the absolute aliasing frequency f_abs=|f_alias|, and restrict f_abs to the interval [0,f_s / 2], and finally output the instantaneous aliasing frequency sequence.

[0030] S3. Based on the Nyquist folding relationship, construct several candidate frequencies for each sampling time using the aliased instantaneous frequency sequence. At the same time, use the known frequency range of the k-clock signal to filter the validity of several candidate frequencies, output the valid candidate frequencies, and form a candidate sequence.

[0031] There are two candidate frequencies, namely the first candidate frequency and the second candidate frequency, denoted as: c1[i] = f_abs[i]; c2[i] = f_s – f_abs[i]; Where c1[i] represents the first candidate frequency at the i-th sampling time, c2[i] represents the second candidate frequency at the i-th sampling time, f_abs[i] represents the absolute value frequency of aliasing at the i-th sampling time in the aliasing instantaneous frequency sequence, and f_s represents the first sampling rate; These two candidate frequencies correspond to different Nyquist regions where the original frequency might fall. Then, the candidate frequencies are filtered for validity using the known frequency range [f_min, f_max] of the k-clock signal, eliminating candidates outside this range. f_min and f_max can be obtained from the factory parameters or pre-calibration of the swept laser.

[0032] In actual data acquisition, due to noise or instantaneous frequency estimation errors, it is possible that both candidate frequencies at a certain sampling moment exceed the known frequency range, or only one candidate frequency is valid. To address these situations, this invention employs the following processing strategy: (1) When the number of valid candidate frequencies at a certain sampling time is 1, the algorithm takes the unique candidate as the candidate state at that time (that is, the state space at that time degenerates into a single state). Dynamic programming is executed normally. Since there is no need to choose between multiple states, the unique candidate is accepted directly. The cumulative cost is calculated as follows: calculate the transition cost of the two states at the previous time to the current unique state respectively, and take the smaller one as the cumulative cost at the current time.

[0033] (2) When the number of valid candidate frequencies at a certain sampling time is 0, an extrapolation is performed based on the frequency value f_est[i-1] recovered from the previous time and the trend of change during the frequency sweep phase (estimated by the instantaneous rate of change Δf_rate through historical frequency difference) to obtain a temporary candidate frequency: c_temp = f_est[i-1] + Δf_rate. If the extrapolated value exceeds the known frequency range [f_min, f_max], it is truncated to the nearest boundary value. At the same time, when calculating the cumulative cost at the current time, an additional penalty value (e.g., 50MHz) is added to the transition corresponding to the extrapolation to reflect the uncertainty of the estimate and prevent the extrapolated frequency from deviating too much from the physical reality.

[0034] (3) If there are no valid candidates for multiple consecutive sampling instants (e.g., more than 5), it is determined that the signal quality of this sweep period is unavailable. The algorithm outputs an error flag and terminates the processing of the current period to avoid generating unreliable resampling results.

[0035] (4) When the number of valid candidate frequencies at a sampling instant is 2, the following dynamic programming algorithm is used to select the optimal frequency path from the candidate sequence to obtain the restored instantaneous frequency sequence. The specific process is as follows: S31. Set the number of state spaces of the dynamic programming algorithm to 2. Define the two candidate frequencies at each sampling instant as the two states of the state space. Set a cost function to describe the cost of the state transition between two adjacent sampling instants, that is, from the state prev at sampling instant i - 1 to the state cur at sampling instant i. The values of prev and cur correspond to the two candidate frequencies at each sampling instant. The cost function includes three components: frequency jump cost, direction penalty term, and acceleration penalty term; 1. The frequency jump cost Δf is defined as: Δf = |c_cur[i] – c_prev[i - 1]| where c_cur[i] represents the state at the i-th sampling instant, and c_prev[i - 1] represents the state at the (i - 1)-th sampling instant; 2. The construction process of the direction penalty term is as follows: By analyzing the change in the envelope slope of the aliased instantaneous frequency sequence f_abs, the boundaries of different stages of the sweep period are estimated. The stage boundaries include: Rising segment. When the candidate frequencies at multiple consecutive sampling instants show an overall rising trend, it is defined as the rising segment. According to the change speed of the rising trend, the rising segment is further divided into a fast rising segment and a slow rising segment; Falling segment. When the candidate frequencies at multiple consecutive sampling instants show an overall falling trend, it is defined as the falling segment. According to the change speed of the falling trend, the falling segment is divided into a fast falling segment and a slow falling segment; Different direction penalties are imposed on different stage boundaries. In the rising segment, if c_cur[i] < c_prev[i - 1] (frequency decreases), the first direction penalty P_up (e.g., 50 MHz) is imposed; in the falling segment, if c_cur[i] > c_prev[i - 1] (frequency increases), the second direction penalty P_down (e.g., 50 MHz) is imposed. [[ID=2s]]

[0036] And a transition zone is set around each estimated stage boundary to reduce the direction penalty weight in the transition zone. For example, the direction penalty weight in the transition zone is reduced to 0.3 times the original value to adapt to the non-linear transition in actual frequency sweeping.

[0037] 3. The acceleration penalty term penalty_acc is defined as a second-order difference penalty: penalty_acc=|(c_cur[i]–c_prev[i-1])–(c_prev[i-1]–c_prev2[i-2])| Where c_cur[i] represents the state at the i-th sampling time, c_prev[i-1] represents the state at the (i-1)-th sampling time, and c_prev2[i-2] represents the state at the (i-2)-th sampling time. When calculating the cumulative cost, the acceleration penalty term is multiplied by a coefficient γ (e.g., 0.1) and then included.

[0038] Set the maximum allowable transition threshold Δf_max, which is adaptively adjusted according to the boundary of the current stage: a larger threshold (e.g., 20MHz) is used in the fast rise and fast fall stages, and a smaller threshold (e.g., 5MHz) is used in the slow rise and slow fall stages. If the frequency transition cost Δf exceeds the maximum allowable transition threshold Δf_max corresponding to the current stage boundary, then state transition is prohibited.

[0039] S32. Calculate the cost after each state transition based on the two states and the cost function. Calculate the cumulative cost across all sampling times using a recursive formula. The recursive formula means that when selecting one state from the two candidate frequencies corresponding to the previous time step, the cumulative cost is minimized when transitioning to a state corresponding to a candidate frequency at the current time step. The recursive formula is: ; Where cost(cur,i) represents the minimum cumulative cost from the initial time to the i-th sampling time, and the state selected at sampling time i is cur; cost(prev,i-1) represents the minimum cumulative cost from the initial time to the (i-1)-th sampling time, and the state selected at sampling time i-1 is prev; the value of prev is any one of the two candidate frequencies corresponding to sampling time i-1; the value of cur is any one of the two candidate frequencies at time i; its single-step transition cost includes three terms: Δf(prev, cur, i) represents the frequency jump cost when sampling time i transitions from state prev to state cur; penalty_trend(prev, cur, i) represents the direction penalty term when sampling time i transitions from state prev to state cur; and penalty_acc(prev, cur, i) represents the acceleration penalty term when sampling time i transitions from state prev to state cur.

[0040] α, β, and γ are weighting coefficients; in this embodiment, α=1, β=1, and γ=0.1. During initialization, the cumulative cost of the valid states at the first sampling time is set to 0. After completing the recursion for all sampling times, backtracking begins from the state with the minimum cumulative cost in the last column to obtain the optimal candidate index for each sampling time, thereby obtaining the recovered instantaneous frequency sequence f_est.

[0041] S33. After completing the recursion of all sampling times, start backtracking from the state with the minimum cumulative cost in the last column to obtain the optimal candidate index for each sampling time, and then obtain the recovered instantaneous frequency sequence.

[0042] S4. Using the recovered instantaneous frequency sequence, the continuous phase sequence of the k-clock signal is reconstructed through cumulative integration. This continuous phase sequence is then used to perform equal-phase-interval resampling on the OCT interferometer signal to obtain an equal-wavenumber-interval OCT resampled signal. The specific process is as follows: Phase reconstruction: Using the recovered instantaneous frequency sequence f_est, the continuous phase sequence φ_recon of the k-clock signal is reconstructed by cumulative integration. Specifically, the phase of the first sampling moment in the continuous phase sequence is initialized to the first value of the aliased phase sequence, that is, let φ_recon[1]=φ_unwrap[1]. Then, for i=1 to M (M=N-1, N is the original number of sampling points), execute φ_recon[i+1]=φ_recon[i]+2π·f_est[i] / f_s.

[0043] S41. Set an equal phase interval Δφ, and generate the target phase sequence using a continuous phase sequence; S42. By using linear interpolation, a mapping is established between the continuous phase sequence and the time axis, and the time point t_resampled corresponding to each target phase is calculated. S43. Interpolate the OCT interference signal to extract the amplitude value corresponding to the sampling time at time point t_resampled, and finally obtain the OCT resampled signal with equal wavenumber intervals.

[0044] Specifically, the typical value of Δφ is 2π, that is, one point is sampled every k-clock period to generate the target phase sequence φ_ticks=φ_recon[1]:Δφ:φ_recon[N]. A mapping is established between φ_recon and the time axis t by linear interpolation, and the time point t_resampled corresponding to each target phase is calculated. Finally, the amplitude value at time t_resampled is extracted from the OCT signal by spline interpolation or linear interpolation to obtain the OCT resampled signal with equal wavenumber intervals.

[0045] The method in this embodiment can be implemented on different hardware, for example: 1. Offline processing mode A data acquisition card equipped with a dual-channel 1GSPS analog-to-digital converter (such as the AlazarTech ATS9373) is connected to a swept-frequency laser (such as the Axsun SS-OCT engine, with a k-clock frequency range of 10MHz to 800MHz). During acquisition, the laser outputs swept-frequency light with a repetition frequency of 100kHz and a cycle duration of 10μs. The data acquisition card simultaneously records the k-clock signal and the OCT interference signal, with 10,000 sampling points per cycle.

[0046] The acquired data was transmitted to a personal computer, and the algorithm of this invention was implemented using MATLAB software. Specific parameter settings were as follows: sampling rate f_s = 1e9 (1GHz), minimum instantaneous frequency of k-clock f_min = 10e6 (10MHz), maximum instantaneous frequency f_max = 800e6 (800MHz), candidate frequency soft constraint lower limit f_min_soft = 5e6 (5MHz), upper limit f_max_soft = 1e9 (1GHz), equal phase interval Δφ = 2π (one complete cycle), maximum allowable jump threshold of 20MHz in the fast segment and 5MHz in the slow segment, directional penalty rise segment P_up = 50MHz, fall segment P_down = 50MHz, transition region length 20ns, acceleration penalty coefficient γ = 0.1. After running the algorithm, an OCT resampled signal with equal wavenumber intervals was output, obtaining approximately 4.05 × 10^6 resampled points per cycle (calculated based on the total phase change). After removing DC, adding a Hamming window, and performing a Fourier transform on the signal, an A-scan depth image was obtained. Compared with the traditional solution using a 2GSPS acquisition card, the depth images obtained by this solution have no significant difference in resolution, while the hardware cost is reduced by about 40%.

[0047] 2. Real-time processing mode (FPGA implementation) The algorithm of this invention was optimized for fixed-point processing and deployed on a Xilinx Zynq UltraScale+ FPGA platform. The Hilbert transform was implemented using an FIR filter, the phase expansion employed a hardware-friendly recursive algorithm, and the dynamic programming two-state recursion used a parallel pipeline architecture. The FPGA operates at a clock speed of 250MHz, and processing one 10μs period of k-clock data takes approximately 12μs (including pipeline delay), meeting the real-time processing requirements of a 100kHz A-scan rate. The recovered frequency sequence and the resampled OCT signal were transmitted to the host computer via a PCIe interface for image reconstruction. This implementation case verifies the real-time processing feasibility of the algorithm of this invention.

[0048] It is understood that this embodiment can operate normally at a sampling rate less than twice the highest k-clock frequency (1GSPS sampling rate processes 800MHz signals). Compared to traditional full-sampling schemes (requiring 1.6GSPS or higher), it can reduce the cost of the analog-to-digital converter by approximately 50%, while also reducing the data transmission bandwidth requirements of the acquisition card. Hardware costs are significantly reduced. This application achieves frequency recovery and resampling entirely through software algorithms, eliminating the need for additional anti-aliasing filters, multi-sampling clocks, or dedicated wavenumbers, and can be directly applied to existing SS-OCT systems using 1GSPS acquisition cards.

[0049] Example 2 like Figure 2 As shown, this embodiment 2 provides an SS-OCT resampling system based on undersampling clock recovery, implementing the method of embodiment 1, including: The SS-OCT signal source includes a linearly swept laser source, a Mach-Zehnder interferometer, and a photoelectric conversion module connected in sequence. The photoelectric conversion module sends the k-clock signal and the OCT interference signal to the acquisition module, respectively. The acquisition module is used to synchronously acquire a k-clock signal and an OCT interference signal using an analog-to-digital converter with a first sampling rate at an undersampling frequency, wherein the first sampling rate is the undersampling frequency of the k-clock signal; The aliasing frequency extraction module is used to perform Hilbert transform on the k-clock signal, construct an analytical signal, extract the aliasing phase sequence from the analytical signal, and perform first-order difference operation on the aliasing phase sequence to obtain the aliasing instantaneous frequency sequence. The phase recovery module is used to construct several candidate frequencies for each sampling time based on the Nyquist folding relationship and the aliased instantaneous frequency sequence, forming a candidate sequence, and then using a dynamic programming algorithm to select the optimal frequency path from the candidate sequence to obtain the recovered instantaneous frequency sequence. The remapping module is used to reconstruct the continuous phase sequence of the k-clock signal by accumulating the recovered instantaneous frequency sequence, and then use the continuous phase sequence to perform equal-phase-interval resampling on the OCT interference signal to obtain an OCT resampled signal with equal wavenumber intervals.

[0050] This application enables frequency recovery of k-clock signals with a maximum frequency of 800MHz and completes equal wavenumber interval resampling of OCT signals using only a 1GSPS sampling rate, thereby reducing system hardware costs.

[0051] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Based on the technical essence of the present invention, any simple modifications, equivalent substitutions, and improvements made to the above embodiments within the spirit and principles of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A resampling method for SS-OCT based on undersampling clock recovery, characterized in that, Specifically, the following steps are included: S1. Simultaneously acquire the k-clock signal and the OCT interference signal using an analog-to-digital converter with a first sampling rate, wherein the first sampling rate is the undersampling frequency of the k-clock signal; S2. Perform Hilbert transform on the k-clock signal to construct an analytic signal and extract the aliased phase sequence from the analytic signal. Then, perform a first-order difference operation on the aliased phase sequence to obtain the aliased instantaneous frequency sequence. S3. Based on the Nyquist folding relation, several candidate frequencies are constructed for each sampling time using the aliased instantaneous frequency sequence to form a candidate sequence. The optimal frequency path is selected from the candidate sequence using a dynamic programming algorithm to obtain the recovered instantaneous frequency sequence. S4. Using the recovered instantaneous frequency sequence, the continuous phase sequence of the k-clock signal is reconstructed by cumulative integration. The continuous phase sequence is then used to perform equal-phase-interval resampling on the OCT interference signal to obtain an OCT resampled signal with equal wavenumber intervals.

2. The SS-OCT resampling method based on undersampling clock recovery according to claim 1, characterized in that, The specific process for obtaining the instantaneous frequency sequence of aliasing is as follows: S21. Extract the wrapped phase from the analytical signal, and then expand the wrapped phase to obtain a continuously expanded aliased phase sequence. S22. Perform a first-order difference operation on the aliased phase sequence, and multiply the result of the first-order difference operation by the first sampling rate to obtain the aliasing frequency. S23. Take the absolute value of the aliasing frequency to obtain the absolute value frequency of aliasing, and finally output the instantaneous frequency sequence of aliasing.

3. The SS-OCT resampling method based on undersampling clock recovery according to claim 2, characterized in that, There are two candidate frequencies, namely the first candidate frequency and the second candidate frequency, denoted as: c1[i] = f_abs[i]; c2[i] = f_s – f_abs[i]; Where c1[i] represents the first candidate frequency at the i-th sampling time, c2[i] represents the second candidate frequency at the i-th sampling time, f_abs[i] represents the absolute value frequency of aliasing at the i-th sampling time in the aliasing instantaneous frequency sequence, and f_s represents the first sampling rate; Simultaneously, the two candidate frequencies are filtered for validity using the known frequency range of the k-clock signal, and the valid candidate frequencies are output.

4. The SS-OCT resampling method based on undersampling clock recovery according to claim 1, characterized in that, The specific process of obtaining the recovered instantaneous frequency sequence is as follows: S31. Set the number of state spaces in the dynamic programming algorithm to 2. Define the two candidate frequencies at each sampling time as the two states in the state space. Set a cost function to describe the cost of state transition between two adjacent sampling times. The cost function includes three components: frequency jump cost, direction penalty term, and acceleration penalty term. S32. Calculate the cost after each state transition based on the two states and the cost function, and calculate the cumulative cost at all sampling times using the recursive formula. S33. After completing the recursion of all sampling times, start backtracking from the state with the minimum cumulative cost in the last column to obtain the optimal candidate index for each sampling time, and then obtain the recovered instantaneous frequency sequence.

5. The SS-OCT resampling method based on undersampling clock recovery according to claim 4, characterized in that, The construction process of the direction penalty term is as follows: By analyzing the changes in the envelope slope of the candidate sequence, the boundaries of different stages of the sweep frequency period are estimated. Different directional penalties are applied to the boundaries of different stages, and a transition zone is set at each stage boundary to reduce the weight of the directional penalty within the transition zone.

6. The SS-OCT resampling method based on undersampling clock recovery according to claim 5, characterized in that, The boundaries of different stages include rising segments and falling segments. When the candidate frequencies of multiple consecutive sampling times show an overall upward trend, it is defined as the rising segment and a first directional penalty is applied. When the candidate frequencies of multiple consecutive sampling times show an overall downward trend, it is defined as the falling segment and a second directional penalty is applied.

7. The SS-OCT resampling method based on undersampling clock recovery according to claim 5, characterized in that, The frequency jump cost Δf is defined as: Δf = |c_cur[i] – c_prev[i-1]| Where c_cur[i] represents the state at the i-th sampling time, and c_prev[i-1] represents the state at the (i-1)-th sampling time; Set the maximum allowed transition threshold Δf_max. Δf_max is adaptively adjusted according to the boundary of the current stage. If the frequency transition cost Δf exceeds the maximum allowed transition threshold Δf_max corresponding to the current stage boundary, then state transition is prohibited.

8. The SS-OCT resampling method based on undersampling clock recovery according to claim 4, characterized in that, The acceleration penalty term is defined as a second-order difference penalty: |(c_cur[i]–c_prev[i-1])–(c_prev[i-1]–c_prev2[i-2])| Where c_cur[i] represents the state at the i-th sampling time, c_prev[i-1] represents the state at the (i-1)-th sampling time, and c_prev2[i-2] represents the state at the (i-2)-th sampling time. When calculating the cumulative cost, the acceleration penalty term is multiplied by the coefficient γ and then included.

9. The SS-OCT resampling method based on undersampling clock recovery according to claim 4, characterized in that, The specific process of step S4 is as follows: S41. Set an equal phase interval Δφ, and generate the target phase sequence using a continuous phase sequence; S42. By using linear interpolation, a mapping is established between the continuous phase sequence and the time axis, and the time point t_resampled corresponding to each target phase is calculated. S43. Interpolate the OCT interference signal to extract the amplitude value corresponding to the sampling time at time point t_resampled, and finally obtain the OCT resampled signal with equal wavenumber intervals.

10. An SS-OCT resampling system based on undersampling clock recovery, characterized in that, include: SS-OCT signal source, used to generate k-clock signal and OCT interference signal; The acquisition module is used to synchronously acquire a k-clock signal and an OCT interference signal using an analog-to-digital converter with a first sampling rate at an undersampling frequency, wherein the first sampling rate is the undersampling frequency of the k-clock signal; The aliasing frequency extraction module is used to perform Hilbert transform on the k-clock signal, construct an analytical signal, extract the aliasing phase sequence from the analytical signal, and perform first-order difference operation on the aliasing phase sequence to obtain the aliasing instantaneous frequency sequence. The phase recovery module is used to construct several candidate frequencies for each sampling time based on the Nyquist folding relationship and the aliased instantaneous frequency sequence, forming a candidate sequence, and then using a dynamic programming algorithm to select the optimal frequency path from the candidate sequence to obtain the recovered instantaneous frequency sequence. The remapping module is used to reconstruct the continuous phase sequence of the k-clock signal by accumulating the recovered instantaneous frequency sequence, and then use the continuous phase sequence to perform equal-phase-interval resampling on the OCT interference signal to obtain an OCT resampled signal with equal wavenumber intervals.

Images