Anchor-Based, Real-Time Motion Compensation for High-Resolution mmWave Radar

In the modern domain of edge sensing and physically compact smart devices, mmWave radar has emerged as a prominent modality, simultaneously offering high-resolution perception capacity and accommodatingly small form factor. The inevitable presence of device motion, however, corrupts the received radar data, reducing target sensing capability and requiring active correction to address the resultant spectral “blurring”. Existing motion compensation techniques utilize computationally intensive post-processing algorithms and/or auxiliary hardware, aspects ill-suited for resource-limited edge devices requiring minimal system latency and complexity. Early works also often consider motion dynamics such as pure single-mode vibration, neglecting additional modes as well as non-harmonic motion content. We resolve both of these limitations by presenting a real-time-compatible, generalized complex motion compensation algorithm capable of correcting multicomponent platform trajectories involving both non-harmonic transients and multimode harmonic vibration. The proposed anchor-based approach achieves average SNR gains of 24.9 dB and 19.7 dB across transient duration and target velocity, respectively, and average multimode harmonic suppressions of 38.9 dB and 29.4 dB across vibration parameters and target velocity, respectively. These results, combined with minimal latency ($\leq \! 240$ ms), low algorithmic complexity, and the elimination of any additional auxiliary sensors, render the proposed method suitable for deployment in typical edge sensing applications.


I. INTRODUCTION
In the realm of edge sensing, mmWave radar has become a modality of choice, simultaneously providing high resolution, impressive maximum range and velocity, and robustness to inclement environmental conditions, yielding distinct advantages over alternative sensing modalities such as light detection and ranging (LiDAR) and cameras.Furthermore, operating wavelengths in the 60-GHz and 77-GHz bands permit minimally sized antennas conducive to on-chip integration, a necessity for area-constrained edge devices [1], [2].Given this consolidation of factors, it is therefore not surprising that mmWave radar silicon ICs are ubiquitous in the design of autonomous vehicles and drones, industrial robots, and wearable devices alike [3], [4], [5].Modern implementations, however, commonly rely upon high-performance back-end multicore processors to analyze the collected sensor data, a paradigm that entails high power/resource consumption with significant processing latency.The reality is that emerging edge devices are resource-limited with real-time response requirements; thus, any power/latency sacrifice is impractical.Navigation devices or vehicles, for instance, require response times of no more than a few hundred milliseconds to account for dynamic object presence in the environment (e.g., pedestrians), necessitating edge-localized decision-making.Edge sensor design, therefore, must simultaneously minimize both hardware/computational complexity and processing latency.
With the oft-used frequency modulated continuous wave (FMCW) radar, even the simplest case of single-mode harmonic motion generates a smeared Doppler spectrum exhibiting pairs of spectral sidebands about the primary signal lobe that could produce SNRs as low as -15 dB at mm-scale vibration amplitudes [23].Given a typical coherent detector false alarm rate of P f a = 10 −6 and detection probability of P d = 0.9, greater than 13-dB SNR is required at the receiver [24], thereby necessitating up to 25-30 dB harmonic motion suppression, with additional correction needed for any nonharmonic transient parasitics.

A. RELATED WORK
While motion compensation is certainly not a novel concept in the radar domain, existing high-resolution methods operate primarily in post-processing, utilizing the entire multidimensional set of data in compute-heavy pipelines.Holistic post-processing, though inarguably advantageous with respect to the resolution and SNR increase afforded by the vast quantity of data, is infeasible for resource-constrained edge devices requiring low-latency, real-time operation.Among existing works, synthetic aperture radar (SAR) motion compensation typically involves the use of computationally intensive algorithms such as the Doppler Keystone transform (DKT), range cell migration correction, and phase gradient autofocus [25], [26], [27], [28], [29], [30], [31].Other higher-performance methods such as parametric DKT [32], range envelope correction [33], short-time-Fourier-transform (STFT) histogram-based interference removal [34], range-dependent map drift correction [35], and ML algorithms [36], [37], [38], among others [25], [39], [40], [41], entail even greater computational load and/or have been exclusively deployed in post-processing.Such methods also tend to implement correction in the range domain only, for which typical range resolutions render mm-scale vibration correction impractical [31], [42], [43].
Moreover, even potentially low-latency models, in addition to many post-processing techniques, rely upon auxiliary hardware such as IMUs to estimate parasitic platform motion parameters, an undesired system hardware overhead and additional source of error for edge devices if used in isolation [39], [44], [45], [46], [47].A number of these promising algorithms also assume motion models such as pure singlemode platform vibration [25], [30], [46], [48] or long-scale translatory motion [49], neglecting the additional modes and high-frequency components present in most realistic applications.One valuable concept that has been recently exploited in the domain of automotive SAR is the idea of stationary ground control reference points, or anchor points, which [50] uses in a weighted least squares estimation to determine the Doppler frequency/velocity error terms associated with the vehicle motion.Combined with the subsequent time-domain phase correction, this algorithm operates on a real-time interval of 200 ms, though it is limited to compensation of constanttranslation-induced artifacts and would not be ideal for more complex motion models or high-frequency transients.
Among non-radar modalities, mobile smartphone optical image stabilization (OIS) techniques often fuse the data from MEMS sensors (e.g., gyroscopes, accelerometers, IMUs, etc.) to detect and correct for physical hand movement.This auxiliary sensor data is routed in direct feedback to local actuators, such as voice coil motors, to automatically adjust the optical lens and stabilize the collected image/video prior to digitization [51].While OIS mitigates image degradation and post-processing complexity, it nevertheless requires secondsscale exposure/integration times and increased hardware.Electronic image stabilization (EIS) and hybrid techniques, on the other hand, function by stitching together frames in the software post-processing, employing methods such as feature point matching [52], optical flow analysis [53], and block matching [54] to assess and compensate for the mobile motion.EIS thus reduces the system hardware involved, at the expense of increased post-processing complexity and greater latency.More advanced techniques have explored the use of magnetic shape memory actuators to simplify system hardware [55], as well as post-processing algorithms including combined global/local motion estimation [56], co-optimized frame stitching and stabilization [57], deep learning [58], and graph optimization [59].With all algorithms, however, complexity and latency tradeoffs remain.
Meanwhile, LiDAR and visual-inertial odometry systems often rely upon IMU fusion [60], [61], [62] and/or computationally intensive correction algorithms [63], [64], yielding the same downsides as radar solutions.The authors in [64], for instance, rely upon a pre-defined quadratic motion model for platform vibration within each chirp interval, adopting a two-step, semi-iterative method to optimally estimate and eliminate both the primary and quadratic vibration components.Both steps utilize an effective frequency-domain correction to correct for ranging errors, though the holistic process would not be suitable for resource-limited applications.
Previously, the authors proposed a real-time, Dopplerdomain motion compensation algorithm for resourceconstrained, edge-based radar sensors, fusing auxiliary data from an IMU [65].However, this work possessed two key limitations: 1) the assumption of purely single-mode vibratory platform motion, and 2) reliance upon an additional auxiliary sensor to determine the parasitic vibration parameters.While the former precludes application to sensor platforms experiencing much more complex, multicomponent motion trajectories, the latter places the system at the mercy of IMU measurement errors and increases the hardware complexity for each local edge device, prohibiting the development of a single-chip mmWave solution for diverse integration scenarios.This work resolves both limitations, as well as those presented by existing motion compensation algorithms.

B. KEY CONTRIBUTIONS
The proposed complex motion compensation algorithm offers the following distinguishing contributions: 1) Demonstrates real-time, generalized complex platform motion compensation for FMCW mmWave radars, addressing composite multimode harmonic motion and nonharmonic transient motion trajectories.This expanded motion model permits applicability to most typical sensor motion scenarios while boasting a real-time-compatible correction latency of less than 240 ms.With respect to the nonharmonic motion component, the proposed algorithm achieves an average 24.9-dBSNR improvement across transient duration and 19.7 dB across target velocity.Meanwhile, with respect to the multimode vibration component, it demonstrates an average 38.9-dB suppression across vibration frequency/amplitude and 29.4 dB across target velocity.Constant false alarm rate (CFAR) evaluation on the pre-and post-correction range-Doppler maps indicate marked improvement in detection accuracy, corroborating the achieved quantitative suppression.

2) Proposes a unique radar anchor point data fusion
technique in the multimode harmonic motion correction step to either complement or eliminate the need for an auxiliary IMU, thus increasing system redundancy and reliability or improving integration flexibility via a reduced hardware overhead.Relative to the performance using an IMU, the alternative anchor-based method yields negligible (< 1-dB) suppression deviation, validating its equivalent utility.3) Achieves motion compensation with minimal computational complexity, permitting feasible integration into a resource-constrained edge IC.Low complexity is achieved via 1) a limited processing window size of N = 2560 chirps, and 2) a novel two-step adaptive incremental correction method, with results demonstrating virtually no (< 0.6-dB) suppression reduction relative to conventional higher-complexity techniques.Collectively, the wide generalizability, low algorithmic and hardware complexity, and real-time nature yield a distinct advantage over existing edge perception systems, with equivalent or superior motion suppression performance.While motion compensation is often associated with SAR, the proposed approach instead targets edge sensing applications for which hardware/resource complexity and response time is critical, including short-range automotive radar target detection and pedestrian micro-Doppler analysis/classification [3], human activity recognition [21], [66], early-warning collision detection and proximity sensing [66], industrial self-navigating autonomous robots, and wearable personal electronics (e.g., for gesture sensing, augmented reality, health monitoring, etc.) [67], among others.Nonetheless, the methods presented are also extendable to SAR applications.Note that the utilized processing interval of 240 ms is selected for compatibility with the minimum real-time response capacity typically required by the aforementioned edge applications [68], [69], [70].
This paper is structured as follows: Section II details the utilized complex motion model and analyzes the theoretical resultant Doppler spectra.Section III describes the proposed motion compensation algorithm, while Section IV provides an outline of the hardware measurement setup used to evaluate the correction efficacy.Compensation results and analysis are presented in Section V, and Section VI concludes with a brief summary and discussion of future work.

II. COMPLEX MOTION MODEL
While complex motion on a radar platform inherently encompasses a wide variety of three-dimensional movement trajectories, those witnessed in typical sensing applications may be broadly grouped into two categories: 1) multimode harmonic motion, and 2) nonharmonic transient motion.The former is prevalent with automotive or aerial radar platforms, for which engine vibration is primarily a composition of multiple vibration frequencies/modes.The latter, meanwhile, accounts for all other nondeterministic motion trajectories or intermittent pulse-like aberrations, prevalent in autonomous navigation applications (e.g., rapid road-or environment-induced bumps and motion artifacts) and in industrial machinery encountering occasional faults.Fig. 1 graphically summarizes this multicomponent complex motion model in a typical autonomous navigation application with a target pedestrian; the depicted environmental anchor point, mentioned in Section I, is utilized in the correction process, as described in detail in Section III-A.
To investigate the expected nature of the mmWave data with such platform motion, we construct a theoretical complex motion model comprising the two aforementioned components: 1) a superposition of harmonic modes with varying amplitude/frequency in differing orientations in 3D space, and 2) non-harmonic, intermittent transient spur motion.

A. MULTIMODE HARMONIC MOTION
Parasitics consisting of multiple distinct modes in typical radar sensing scenarios often stem from complex vibration induced by engines or motors.While theoretically, the number  of such modes may be limitless, most applications, such as automotive sensing and navigation, robotics, and industrial machinery, are typically witnessed to involve no more than three or four significant primary modes.Spectral components at all other frequencies are generally of much lower magnitude and of limited impact on the signal [11], [12], [13], [14], [15], [17], [18], [19], [20].In our complex motion model, and in the correction algorithm detailed in Section III, we assume that these ∼3 primary modes occur between 0 Hz and 70 Hz, with corresponding harmonics and intermodulation spurs permitted to fall outside this range.Apart from referencing [11], [12], [13], [14], [15], [17], [18], [19], [20], these assumptions on the number and frequency of the primary modes were also validated by measuring the vibration of an actual automotive engine using a high-resolution accelerometer.Fig. 2 depicts the resulting triple-axis vibration spectra for two different engine RPMs, verifying that the significant peaks, or primary vibration frequency modes, do indeed occur at frequencies less than 70 Hz and in limited quantity (≤ 3), with any additional peaks representing harmonics or intermodulation spurs.

B. TRANSIENT SPUR MOTION
Any platform motion that is not harmonic in nature, neglecting continuous lateral translation, is likely to be a sharp transient or pulse movement.The assumed narrow durations are often the result of gearbox or bearing faults, road-induced bumps, etc., and are typically characterized by high-frequency components [12].We make no assumptions on frequency content or amplitude of the transient spur in either the analysis that follows, or in the correction algorithm itself.

C. THEORETICAL COMPLEX MOTION RADAR SPECTRA
To analyze the expected effect of complex platform motion on the received radar data, we construct a mathematical model comprising the described components.For simplicity, we consider dual-mode harmonic motion along the coordinate axis (e.g., z-axis) corresponding to the radar line of sight (LOS), but the analysis easily generalizes to higher-order modes.The amplitude, angular frequency, and initial phase of vibration mode i are denoted by A v,i , ω v,i , and φ 0,i , respectively.Meanwhile, the nonharmonic transient pulse signal component is modeled by a windowed sinusoidal peak with time delay t d,p , duration τ p , and amplitude A p .While this selected model emulates the rapid motor revolution used in the actual hardware setup (Section IV), any pulse function (e.g., Gaussian) impacts the resultant data similarly.The platform motion trajectory can therefore be expressed as a function of time via the superposition of these components: The theoretical impact of such motion upon the radar data may be analyzed directly, as shown in the Appendix, yielding a Doppler spectrum characterized by multimode-harmonicmotion-induced sidebands, containing both harmonic and intermodulation spurs, superposed upon a transient-induced raised noise floor.Fig. 3 depicts the theoretical motion waveforms and corresponding radar spectra using the above model for dual-mode harmonic motion with a single transient, confirming the expected frequency characteristics.While this example utilizes a stationary target for simplicity, nonzero target velocity would simply induce a frequency shift.

D. TAKEAWAYS FROM MOTION MODEL ANALYSIS
Examination of the theoretical received radar spectra resulting from the superposed multimode harmonic motion and transient informs the appropriate approach to inline parasitics correction.The sidebands resulting from the multimode harmonic motion, for instance, may be inverted via the two-step process of 1) estimating the constituent modes, and 2) formulating a corresponding frequency-domain deconvolution kernel to apply to the data.This kernel is similar in spirit to the single-mode vibration compensation method from [65], but of necessarily greater complexity to suppress all additional modes, harmonics, and intermodulation spurs.The broadband dynamic range reduction induced by the transient, however, is not trivially corrected in the frequency domain and thus requires some level of time-domain treatment.

III. COMPLEX MOTION COMPENSATION ALGORITHM
The proposed complex motion compensation algorithm borrows somewhat from the single-mode vibration correction concepts introduced in the authors' previous work [65], but with several significant performance-enhancing changes and additions conducive to greater motion generalizability and hardware integrability.In particular, we demonstrate 1) utilization of a stationary anchor point, rather than an auxiliary IMU (used here for ground-truth performance comparison only), to achieve commensurate compensation performance without any additional hardware, and 2) a real-time motion compensation algorithm capable of significant transient and multimode sideband suppression with simultaneously low computational complexity suitable for integration into a resource-limited edge sensor device.Real-time operation is achieved by operating on nonoverlapping temporal data windows of 240 ms, corresponding to N = 2560 chirps, the minimum interval required for sufficient Doppler resolution.
The proposed approach is decomposed into two independent steps to address the transient and multimode harmonic motion components separately.As analyzed in Section II-C, since the transient motion components will not manifest themselves as separable tones in a purely frequency-domain data representation, they require a time-frequency analysis to suppress.Following the transient removal, the multimode harmonic motion is addressed in the frequency domain via a deconvolution kernel capable of inverting the sidebands induced by the various vibration modes and their intermodulation components.Fig. 4 depicts a high-level diagram of the proposed algorithm, described in the following sections.

A. USING A STATIONARY ANCHOR POINT
The first step of the proposed algorithm extracts the anchor point RX signal, as shown in Fig. 4.While the results of [65] demonstrate that an IMU can be effectively used for vibration parameter estimation, the presence of an auxiliary device increases edge system hardware complexity, precluding the development of a convenient single-chip integrated radar sensor/processor, while also introducing its own measurement error.As an alternative, we propose a novel anchor-based motion compensation technique.An anchor point may be defined as a reference object in the environment exhibiting a Doppler profile with distinct and known properties in the absence of platform motion.Note that although we demonstrate an anchor-based approach that achieves equivalent performance to that with an IMU and renders the latter optional, future exploration using the anchor in conjunction with an IMU or other auxiliary sensor could boost performance and system reliability/robustness even more via the increased redundancy afforded by two independent measurements.
The most convenient anchor point is a stationary target, which is simple enough to identify in a typical sensing scenario (e.g., the ground, a building, a signpost, a static point target, etc.), and for which the ideal Doppler profile is simply a single peak at zero Doppler.Since any superposed parasitic platform motion that corrupts the primary target's Doppler profile will corrupt the anchor point profile in a commensurate manner, the platform motion characteristics may be inferred by simply comparing the received anchor spectrum with the ideal single-peak spectrum, without the hardware overhead of an auxiliary inertial sensor.When the frequency and amplitude information embedded in this inferred parasitic motion is effectively "subtracted" out, only the desired target profile remains, including any inherent target micro-Doppler features.
Note that the motion artifact similarity in the target and anchor spectra assumes sufficiently high radar-cross-sectiondependent signal returns from both objects to permit easy identification in the presence of clutter.The utilized colinearity with the radar LOS is unnecessary, however, and in any case, angular offsets would not impact the detected vibration frequencies in the anchor spectrum and may be easily calibrated out.Although such angular offsets would affect the relative signal amplitudes of different modes, this has only a small impact on the maximum achievable suppression and may be resolved in future work by referencing multiple anchor points, as [50] does, and performing a best-fit across all.Such a future approach would also be expedient for highly dynamic (e.g., marine) environments that lack stationary anchor points.In particular, a best-fit across multiple non-stationary anchors would extract the parasitics common to all, representing the platform-induced motion.In this work, we evaluate the proposed motion compensation algorithm for both the anchor-and IMU-based estimation methods, independently, to demonstrate that equivalent performance is achievable solely using the stationary anchor.

B. TRANSIENT MOTION CORRECTION 1) OVERVIEW
The notion of transient spur correction has been well studied; myriad applications in the domains of communications, image processing, and sensing have investigated methods of such rapid "impulse noise" removal [71].As demonstrated with the analysis of Section II-C, transient motion components will not manifest themselves as single tones in a Fourier spectrum of the observation window data, rendering pure frequency-domain correction (e.g., via low-pass/bandpass filters) infeasible.On the other hand, pure time-domain approaches, such as median filters [72], [73] and poolingbased filters [74], are often simple to implement, but their lack of frequency selectivity is a drawback for radar Doppler signals in the event that the transient/impulse frequency content coincides with or falls below that of the target Doppler or vibration lobes.Time-domain ML approaches, meanwhile, are unnecessarily complex [75], [76].The upshot, therefore, is that a combined time-frequency analysis is required.
With this in mind, we operate in the Wavelet domain, relying upon coefficient thresholding to isolate the signal associated with the transient across both time and frequency [77], [78], [79].Unlike an STFT, a Wavelet transform uses an adaptive frequency resolution, with greater resolution at lower frequencies.Since the target Doppler and vibration signals are low frequency, Wavelet analysis provides significantly greater resolution than the STFT for an equal number of frequency bins.Wavelet-based approaches have gained significant traction for the purpose of transient detection and/or filtering [80], with numerous works in the domains of seismic monitoring [81], [82], arrhythmia detection [83], [84], and machinery fault and quality event identification [85], [86] showcasing the performance gains afforded by Wavelet domain processing.Promising results have also been demonstrated for Doppler radar transient singularity detection specifically, often via the discrete Wavelet transform (DWT) and continuous Wavelet transform (CWT) [87], [88].The proposed approach employs a thresholded 1D dual-tree complex Wavelet transform (DTCWT) to suppress the transient spur signal components.The multi-channel 1D DTCWT operates on the complex time-domain signal in each radar range bin and provides several advantages over alternative Wavelet-based methods.Contrary to a standard DWT, the DTCWT is suitable for complex signals, with only a 2x computation increase.Furthermore, it is approximately shiftinvariant; coefficient energy distributions across scales do not depend upon the transient's time of occurrence, yielding greater performance stability across temporal location.Finally, unlike with the DWT, the DTCWT-based signal decomposition/reconstruction is not subject to aliasing artifacts in the reconstructed signals in the event of coefficient thresholding between analysis and synthesis [89], [90], [91].Meanwhile, relative to the complex-signal-compatible CWT, the DTCWT's base-2 scales are more conveniently implemented in hardware than the CWT exponential scales, affording significantly lower computational complexity.To maximize filter simplicity without sacrificing performance, we utilize a nearsymmetric biorthogonal filter pair with scaling filter of length 5 and Wavelet filter of length 7 for the first decomposition level and an orthogonal Q-shift Hilbert wavelet filter pair of length 10 for subsequent levels.
2) APPROACH Fig. 5 summarizes the transient correction algorithm.The length-N complex input signal in the mth processing window (s m ) is first decomposed into its analytic (s + m ) and antianalytic (s − m ) components via a pair of Hilbert transforms to account for the asymmetric nature of its Fourier spectrum about zero Doppler.The signal components are subsequently transformed via a DTCWT, using K = log 2 (N ) decomposition levels, yielding the pair of sets V + and V − (together, Visualizing the coefficient magnitude scalogram allows for trivial identification of the transient, localized in both time and frequency and separate from the anchor-, target-, and vibration-induced signals. Inverted hard thresholding is then applied to the union set V ± , nulling the coefficients exceeding a predefined magnitude threshold T , to suppress the transient.This thresholding is simultaneously scale-isolated, operating only in the Wavelet frequency scales not containing either a) the anchorcentric Doppler neighborhood, or b) the target-induced signal.Neglecting the former scales ensures preservation of the multimode vibration sidebands flanking the anchor point zero-Doppler frequency bin, thus permitting accurate anchor-based mode estimation, as discussed in Section III-A.Given the maximum expected primary mode frequency f max v = 70 Hz, the set of coefficients in these scales may be defined by: Meanwhile, neglecting the latter scale ensures preservation of the target signal, corresponding to the coefficient set: where the sgn function helps denote (c + i ) j or (c − i ) j for positive or negative target Doppler, respectively, and where the target velocity Doppler frequency f d,target is estimated using the maximum peak from the corrected spectrum in the previous processing window m−1.This f d,target approximation is valid since even a slightly deviated velocity in the span of a single window likely maps to the same scale, given the coarse, base-2 logarithmic frequency discretization utilized by the DTCWT.Unlike with the anchor Doppler bin, suppressing the f d,target -centric vibration sidebands in this step would be welcome indeed.However, typically, these sidebands are inevitably contained within the target Doppler scale given the coarse decomposition granularity, which yields a wide signal frequency band extent of up to 160 Hz, depending upon the target velocity.This inclusion of the parasitic sidebands thus necessitates the subsequent multimode harmonic motion correction step.The final coefficient set isolated for the thresholding operation is therefore: To conclude the transient correction step, the transientfiltered time-domain signal s m is reconstructed via an inverse DTCWT using the thresholded Wavelet coefficients and a final analytic/anti-analytic component recomposition.
Note that the proposed method, while effective, is nevertheless limited in its suppression capacity by the DTCWT's decomposition granularity gradient, as the enforced scale isolation may preclude optimal impulse filtering across a wide band if the transient and target Doppler Wavelet scales coincide at high frequencies.The other tradeoff with the proposed time-frequency-domain approach relates to its complexity; while superior performance-wise to pure frequency-domain band-pass or time-domain median filters, it nevertheless incurs a greater computational load.

C. MULTIMODE HARMONIC MOTION CORRECTION
With any nonharmonic transient motion components filtered out, the only remaining significant motion artifacts are due to the multimode vibration, which creates harmonic and intermodulation sidebands about the primary signal lobes.These sidebands are corrected in the frequency domain, allowing for computationally cheap O(N ) pointwise multiplication with the formulated deconvolution kernel.The proposed method entails: 1) O(N ) signal frequency content estimation via an adaptive, incremental orthogonal matching pursuit (OMP) technique free of any O(N log(N )) discrete/fast Fourier transform (DFT/FFT) operations, and 2) deconvolution kernel generation and application.

1) STEP 1: SIGNAL FREQUENCY CONTENT ESTIMATION VIA ADAPTIVE ORTHOGONAL MATCHING PURSUIT
Estimation of the constituent vibration frequencies and amplitudes using the stationary anchor point is achieved via OMP [92], [93], which provides several advantages over alternative spectral estimation techniques.A simple baseline FFT, for instance, would not provide sufficient frequency resolution in the event that vibration modes occur in adjacent Doppler bins.OMP, on the other hand, may be implemented upon an oversampled complex Fourier basis, a technique that [94] refers to as discretized OMP, yielding greater resolution and motion suppression capacity than an FFT.This Doppler resolution gain is critical given the real-time reliance upon short, 240-ms windows of data.Coincidentally, a direct FFT analysis in each successive processing window also entails greater complexity than the proposed approach, as discussed in Section III-D.While alternative high-resolution spectral estimation subspace methods (e.g., ESPRIT, MUSIC, atomic norm soft-thresholding (AST), (Newtonized) Lasso optimization, Newtonized OMP, etc.) may achieve slight performance gains relative to the utilized oversampled OMP technique [94], [95], [96], [97], [98], [99], OMP is advantageous in its simplicity and fast implementation.
However, for real-time sensing on a resource-limited edge device, algorithm complexity reduction is critical.Thus, we propose three alterations to conventional oversampled OMP: r Bounded frequency grid: The OMP complex exponential Fourier atom dictionary space is bounded to reflect the maximum expected mode frequency, f max v = 70 Hz.r Two-step incremental OMP: The OMP procedure is decomposed into a two-step process, with a coarse estimation step, followed by an incremental fine estimation step.The coarse step operates on a DFT-defined, non-oversampled Doppler grid, yielding rough initial estimates of the signal frequency content.This is succeeded by a fine step operating on the oversampled complex Fourier basis, but with a dictionary containing only the atoms in the near neighborhood of the coarseestimated frequencies.By limiting the search space on the oversampled grid and utilizing a minimal sparsity level of L = 3, complexity is significantly reduced.
r Adaptive estimation: Rather than an independent OMP application to each successive length-N processing window, a closed-loop weighting mechanism (Fig. 6) is utilized to adapt to either fast-or slow-changing signal frequency content over time, thus avoiding the computation cost of a repeated coarse estimation step.The output of the harmonic mode content analysis in the mth data window is a vector of frequency estimates x m ∈ R L .
Coarse Estimation Step: As illustrated in Fig. 6, rather than perform an OMP-or FFT-based sparse signal approximation, the coarse estimation step in the current processing window m computes a linearly weighted average of the frequencies estimated in the previous M processing windows.Weights w = {w i } ∈ R M are determined by the rate of change of the signal frequency content, computed as a function f of the RMS error x between the most recent frequency estimate vector x m and the coarse approximation x from which it was derived.Sufficiently small change x < ε l in the signal approximation maps to a uniformly distributed weight vector, while a sufficiently high RMS deviation x ≥ ε h results in maximum weight assigned to the most recent window, thus accounting for the rapidly changing constituent signal modes.
The weight distribution for ε l ≤ x < ε h is simply determined via a linear gradient between the two edge cases, i.e., in terms of the indicator function 1{. ..}: , (5) for i ∈ [1, M].This computationally cheap linearly weighted average isolates the initial L coarse DFT-resolution frequency indices, represented by the individual spectral peaks of Fig. 6.Given the enforced frequency grid bounds, this vector of L frequency indices,

Fine Estimation
Step: Following the coarse estimation step, the fine estimation OMP step takes as an input the transient-filtered anchor point signal (s anchor ) m in the current processing window m and operates on a β-oversampled Fourier atom Doppler grid, with a frequency-bounded orthonormal dictionary D containing only the complex Fourier atoms, temporally indexed by n ∈ [0, (βN − 1)], corresponding to a P-radius neighborhood, ρ = [−P, −(P − 1), . . ., P], about each of the Doppler indices in k coarse , i.e., This truncated, fine-resolution dictionary space is indicated with the red shaded regions in the spectrum of Fig. 6.To minimize computation without sacrificing performance, we utilize β = 2 and P = 2.The fine estimation step then computes the optimal incremental nonlinear approximation to (s anchor ) m using D, with a maximum of one selected atom (ρ l ) i within each set {k i + ρ l }.In other words, the transient-filtered anchor signal input is "matched" to a Fourier-dictionary-based superposition of atoms given by: Hence, given an ideal multimode harmonic motion signal: where ω v,i represents the angular vibration frequency of the ith mode, the fine estimation step yields a magnitude-L set of dictionary atom indices, κ = {κ i : i ∈ [0, (L − 1)]}, corresponding to frequencies f fine vib in Figs. 4 and 7, where: These indices correspond to the set of Fourier atoms: 2) STEP 2: DECONVOLUTION KERNEL GENERATION AND APPLICATION Having estimated the sparse Fourier atoms comprising the complex multimode harmonic motion anchor signal, the proposed pipeline proceeds with the sideband compensation process.Fig. 7 illustrates the proposed algorithm.To formulate the deconvolution kernel needed to invert the parasitic sidebands, the algorithm first computes the harmonic and intermodulation components associated with these atoms, which manifest themselves quite prominently in the Doppler spectrum of the IF signal returns, as demonstrated in Section IV.For computational simplicity, only the second harmonic and first-order intermodulation components are computed.Hence, to the set of atoms A (m) OMP generated by the frequency estimation stage, we append the additional L + L(L − 1) = L 2 harmonic and intermodulation Fourier components: for i, l ∈ [0, (L − 1)] and i = l.Reindexing with μ, we thus obtain the final magnitude-(L + L 2 ) union set given by: H∪IM The zero-Doppler-centered deconvolution kernel H (m) 0 ( jω) in processing window m is constructed via a cascaded multinotch filter composed of |A (m)  est | second-order, digital IIR filter sections, with notches at the frequencies corresponding to the atoms of A (m)  est and quality factors determined by the anticipated Doppler frequency migration.Intuitively, the Q-factor determines each notch bandwidth, which must be dynamically adapted to invert the induced Doppler leakage at each sideband lobe.Referring to the parameter definitions in Table 1, the migration term f dm,i , corresponding to atom e j2πμ i n/(βN ) , is maximized at f = B/2, and can be computed using the peak vibration mode velocity |v max v,i |=|a v,i |/ω v,i , given the associated acceleration |a v,i | and frequency ω v,i = 2πμ i /(βNT c ).
The notch bandwidth ω,i for atom e j2πμ i n/(βN ) is simply twice this Doppler migration frequency, yielding a Q-factor: where the scaling factor α Q is selected to account for window resolution and to optimize performance.Now, while a v,i is easily provided by the accelerometer in the case of sensor fusion utilization, relying solely upon the anchor point, as proposed, necessitates analysis of the anchor sideband amplitude S v,i corresponding to mode i in the Doppler spectrum.However, since S v,i is linearly proportional to |a v,i |, i.e., S v,i = γ |a v,i |, the effect upon Q i is trivial: The formulated kernel is then shifted in Doppler according to ( f d,target ) m−1 for the target located in range bin r target ∈ [1, ADC ], generating the final deconvolution filter: Frequency-domain application of the constructed deconvolution kernel to the transient-filtered target signal (s target ) m , followed by reconversion to the time-domain, yields the final corrected signal (s * target ) m in window m: The complete 2D array of corrected data is obtained via: 1) Duplication of H (m) 0 ( jω) across range bins, shifted in Doppler for each range bin r ∈ [1, ADC ] according to the associated ( f d,target ) m−1 to obtain H (m)  r ( jω).This generates a 2D multinotch deconvolution kernel across both the range and Doppler dimensions, permitting generalization of the algorithm to multiple targets and anchor points in the scene.
2) Simple frame concatenation and regularization for all m along the time/Doppler dimension.

D. MULTIMODE HARMONIC MOTION CORRECTION TIME COMPLEXITY ANALYSIS
To analyze the complexity reduction achieved via the proposed adaptive two-step OMP algorithm, we begin by noting that a non-adaptive implementation would utilize a fine estimation OMP step identical to that of the adaptive version, with a time complexity of O(N|D| + Nk + Nk 2 + k 3 ) in the kth OMP iteration [100].In the proposed algorithm, however, the dictionary size in the fine estimation step is limited to only the near neighborhood of the L coarse-estimated Fourier atoms, and so |D| N. Furthermore, since k ≤ L, it follows that k N. Hence, our reduced-complexity fine estimation step entails only O(N ) complexity.
The initial coarse estimation step, however, differs between the adaptive and non-adaptive methods.In particular, the adaptive method coarse step simply entails a weighted average of the L frequency estimates in the previous M processing windows, yielding an O(LM ) ∼ O(1) time complexity.On the other hand, a non-adaptive implementation would require independent spectral analysis in the initial coarse step.Yet, toward this end, even a simple FFT, the least-intensive method of frequency estimation, requires O(N log(N )) time complexity.
Combining the coarse and fine estimation steps in the two cases thus yields an O(N ) complexity for the proposed adaptive algorithm, compared with worse O(N log(N )) complexity for a non-adaptive implementation.Since the ensuing frequency domain correction simply involves an O(N ) pointwise product, the initial coarse step of a non-adaptive implementation would indeed dominate the deconvolution kernel generation and application process.Thus, the adaptive method is significantly more resource efficient.

IV. SYSTEM MEASUREMENT SETUP
To generate complex platform trajectories consisting of both multimode harmonic motion and nonharmonic transients, the hardware setup of Fig. 8 is designed.As illustrated, the Texas Instruments (TI) AWR1843 mmWave radar module is secured inside an open-face enclosure, with a TDK InvenSense ICM-20948 IMU situated directly behind the TX/RX antenna array to provide continuous accelerometer readings of the induced vibration for ground-truth estimation and algorithm performance comparison.Mounted around the periphery of the enclosure are four eccentric rotating mass (ERM) motors, three of which are driven to provide continuous multimode vibration with up to three modes and oriented along various spatial planes.Of these three motors, two supply low-frequency vibration with a high amplitudeto-frequency ratio, while the third provides high-frequency vibration with a lower amplitude-to-frequency ratio.As described in Section II-A, the number and frequency range of these motors was selected to reflect the typical vibration profiles induced by engines in automotive, robotics, and industrial sensing/navigation applications.The final fourth motor is driven by pulses with varying duration to generate intermittent transients with amplitudes sufficient to significantly impair the radar data.Combining all motors thus yields complex platform motion trajectories with up to three harmonic modes and occasional high-amplitude pulses.
Meanwhile, the scene of interest contains two targets: 1) a primary target of interest with programmable velocity along the radar LOS (z-axis), as controlled by a precise linear stage, and 2) a stationary corner reflector target acting as an environmental anchor point.A TI CC2652R1 microcontroller receives continuous data from the IMU, while triggering the ERM motor driver circuits appropriately to induce the complex platform motion.Finally, a central embedded application (laptop) synchronizes all data collection from the radar and IMU with the motor and linear stage control signals.Table 1 details the relevant scene characteristics and operating parameters utilized for all hardware units.Fig. 9 depicts an example measured radar Doppler spectrum for the platform undergoing triple-mode harmonic motion with a moving target (+0.3m/s), along with the corresponding accelerometer spectrum.As shown, the three primary modes present in the accelerometer spectrum correlate precisely with  the sidebands present around both the moving target and stationary anchor point signal lobes, demonstrating the utility of the IMU for accurate ground-truthing.Furthermore, the location of the sidebands in the Doppler spectrum match that predicted from the theoretical analysis in Section II-C, with frequency spacings around the target and anchor point lobes including both harmonics of the primary modes, as well as the intermodulation components.

V. MEASUREMENT RESULTS AND DISCUSSION
To evaluate the effectiveness of the proposed complex motion compensation algorithm, it is necessary to consider the transient and multimode harmonic motion components separately, since their contrasting manifestations in the resulting radar spectra and range-Doppler profiles require different assessment metrics.In particular, the transient motion compensation efficacy is measured by considering the achieved SNR improvement in the worst-case single radar frame of duration T f , which yields an unbiased, interval-independent suppression lower bound.On the other hand, since the multimode harmonic motion manifests itself as sidebands in the Doppler spectra, its suppression is assessed via the improvement in spurious free dynamic range (SFDR) to capture the achieved sideband power reduction.

A. TRANSIENT MOTION SUPPRESSION RESULTS
Fig. 10 depicts example CWT scalograms showcasing the performed transient spur suppression, for various illustrative combinations of either one or multiple modes, with either stationary or moving targets.While the transients are easily distinguishable in the original signals, the corrected images demonstrate a significant degree of filtering.To quantify this suppression, we first note that, since the transients are generated using an ERM, for which amplitude is directly proportional to frequency, a longer-duration transient has inherently lower magnitude and thus, generates limited noise, yielding less room for SNR improvement.This ERM-motorinduced duration-amplitude relationship across measurements is plotted in Fig. 11(a).Hence, it is most appropriate to visualize the suppression (a range of 17.9-30.8dB, with an average of 24.9 dB, across transient duration) normalized against transient amplitude, as in Fig. 11(b), demonstrating an SNR gain that increases with duration.This trend reflects the fact that a longer duration generates noise across a wider time span, permitting greater room for improvement.
Fixing the transient duration and sweeping across target velocity yields the curve of Fig. 11(c   reduction at greater velocity magnitudes (and hence, greater main lobe Doppler frequencies), may be attributed to a corresponding increased main lobe colocation with the transient Wavelet scales, rendering transient removal less effective given the use of signal scale isolation.In this case, portions of the transient would be corrected in the harmonic motion  compensation step via the frequency-localized multinotch deconvolution filter, yielding equivalent total suppression.

B. MULTIMODE HARMONIC MOTION SUPPRESSION RESULTS
The performance of the multimode harmonic motion compensation step is evaluated on the transient-filtered signals to independently assess its efficacy, and results are compared with techniques utilizing IMU-enabled and/or non-adaptive frequency content evaluation.Fig. 12 qualitatively illustrates the correction for various mode content with both stationary and moving targets, indicating significant suppression of all harmonic and intermodulation sidebands in the corrected spectra.Investigating the suppression extent across vibration frequency/amplitude for the single-mode vibration case with each of the three motors yields the plots of Fig. 13.The results for all vibration mode combinations across frequency are reported in Table 2, yielding a composite average 38.9-dBSFDR improvement with the proposed method, relative to 39.9-dB and 39.4-dB averages for the IMU-based and nonadaptive implementations, respectively.From both Fig. 13 and Table 2, it is clear that the proposed adaptive, anchor-based correction algorithm achieves nearly equivalent SFDR improvement to the IMU-based and/or non-adaptive realizations, verifying that reduced hardware and computational complexity may be achieved without sacrificing performance.Note that the higher performance for modes M [3] and M [2,3] are a consequence of the fact that the high-frequency Motor 3 spurs represent the dominant sidebands in these cases only, and the inherently lower spectral leakage at these higher frequency offsets relative to the primary signal lobe permit greater room for SFDR improvement.
Finally, Fig. 14 illustrates the achieved SFDR improvement across target velocity for all three mode combinations, with an average 29.4-dB sideband suppression for the anchor-based, adaptive algorithm.Again, the IMU-based (29.7 dB) and nonadaptive (29.2 dB) methods offer similar performance despite increased hardware/computational complexity, validating the utility of the proposed method.Relative to existing works, the presented approach is uniquely advantageous in its collective real-time nature, low computational and hardware complexity, and complex motion generalizability, with either comparable or superior suppression performance, as summarized in Table 3.Note that [25] and [27] present motion compensation approaches for SAR applications, which require neither realtime nor low-latency operation.Nevertheless, these works still provide an apt point of reference for the achievable motion suppression relative to computational complexity, verifying the presented algorithm's potential for edge sensing applications.

C. CFAR DETECTION RESULTS
We conclude by applying the proposed complex motion correction algorithm to a cell-averaging (CA) CFAR detector to evaluate the detection accuracy gains afforded by the platform motion suppression.Fig. 15 depicts example 2D range-Doppler profiles prior to and following application of the proposed algorithm.While the original range-Doppler maps exhibit evident multimode-harmonic-motion-induced sidebands and transient-induced noise spread, the corrected maps illustrate successful localization of the anchor and primary targets in the scene.This enhanced target localization is corroborated by the quantitative CFAR detection results presented in Fig. 16, which plots the average false discovery rate (FDR) as a function of the CFAR threshold factor for all different mode and velocity combinations under constant training band and guard band size conditions.Indeed, it is clear that the correction algorithm significantly improves the CFAR performance, with the associated FDRs dropping to zero at a drastically increased rate relative to the uncorrected detection curves.Note that transient parameters are not swept in the plots of Fig. 16, since the CFAR detection results do not

VI. CONCLUSION
This work has presented a real-time algorithm for complex platform motion compensation in edge-based mmWave radar systems.The proposed approach's combined low computational complexity, minimal hardware (i.e., no auxiliary IMU), real-time operation compatibility, and wide generalizability to typical sensor motion trajectories yield a distinct advantage over existing methods for resource-limited edge devices, with comparable or superior correction performance.Evaluation using a CA CFAR detector corroborates the algorithm's efficacy for real-time detection and tracking applications.Now, although the achievements of this work represent a significant advancement in the domain of high-resolution edge sensing, the proposed algorithm remains to be implemented on an actual edge device.Furthermore, while an extension to multiple targets and a translating platform is realizable, the current work does not demonstrate corresponding results, opting instead to spotlight the primary algorithmic contributions.Measurements and results also utilize a single RX channel, thereby ignoring the angular dimension afforded by an antenna array, which would be necessary to resolve and correct multiple targets at equal range, as well as the event of Doppler signal overlap with the sidebands of an adjacent target.Future work can expect: 1) extended measurements to address the aforementioned limitations, and 2) evaluation of the algorithm's physical resource efficiency (e.g., power, area, etc.) and the achievable detection/tracking effectiveness when deployed in real time on an edge IC.A hardware implementation, combined with the substantial analytical framework already presented, would pave the way for a next generation of ultra-efficient, high-performance mmWave edge devices.

APPENDIX
Here, we derive the theoretical complex motion Doppler spectra resulting from the motion trajectory described by (1).Given a target with retreating velocity v z at an initial static range R 0 , the dynamic range is given by: R(t ) = R 0 + v z t + A v,1 sin ω v,1 t + φ 0,1 The Doppler spectrum of the radar signal can be evaluated by assuming zero-degrees azimuth and elevation angles for the motion displacement in the radar platform reference frame, with nonzero angles simply inducing trigonometric scaling factor constants.For a fast-chirp FMCW radar with a start frequency f c and operating wavelength λ, and assuming a target reflectivity ρ, the slow-time component of the RX signal, used for Doppler processing, may be written as: s RX (t ) = ρe = ρe The corresponding downconverted IF signal is simply: Here, ψ (t ) represents the signal component due to the multimode harmonic motion, while ξ (t ) represents the component to the transient pulse.The remaining exponential factor is a function of target characteristics and simply induces a scaling and frequency shift in the resultant spectra.

FIGURE 1 .
FIGURE 1. Complex motion example: vehicle-mounted radar undergoing multicomponent motion, corrected in real-time using the anchor data.

FIGURE 4 .
FIGURE 4. Top-level diagram of the proposed complex motion compensation algorithm.

FIGURE 6 .
FIGURE 6. Diagram of the proposed two-step, incremental, adaptive OMP frequency content estimation procedure.

FIGURE 7 .
FIGURE 7. Proposed method of deconvolution kernel generation.

FIGURE 8 .
FIGURE 8. (a) Hardware setup designed to generate multidimensional complex motion with synchronized mmWave radar and accelerometer data collection, (b) photo of motor-moutned radar enclosure system, (c) close-up photo of mmWave and IMU boards, and (d) controller hardware photo.

FIGURE 9 .
FIGURE 9. (a) Example measured Doppler spectrum for triple-mode harmonic motion, and (b) the corresponding accelerometer spectrum.
Fig.10depicts example CWT scalograms showcasing the performed transient spur suppression, for various illustrative combinations of either one or multiple modes, with either stationary or moving targets.While the transients are easily distinguishable in the original signals, the corrected images demonstrate a significant degree of filtering.To quantify this suppression, we first note that, since the transients are generated using an ERM, for which amplitude is directly proportional to frequency, a longer-duration transient has inherently lower magnitude and thus, generates limited noise, yielding less room for SNR improvement.This ERM-motorinduced duration-amplitude relationship across measurements is plotted in Fig.11(a).Hence, it is most appropriate to visualize the suppression (a range of 17.9-30.8dB, with an average of 24.9 dB, across transient duration) normalized against transient amplitude, as in Fig.11(b), demonstrating an SNR gain that increases with duration.This trend reflects the fact that a longer duration generates noise across a wider time span, permitting greater room for improvement.Fixing the transient duration and sweeping across target velocity yields the curve of Fig.11(c), indicating an average SNR improvement of 19.7 dB.Note that the performance

FIGURE 10 .
FIGURE 10.Pre-and post-correction example CWT images illustrating the performed transient spur suppression.

FIGURE 12 .
FIGURE 12. Pre-and post-correction example spectra illustrating the performed multimode harmonic motion suppression.

FIGURE 13 .
FIGURE 13.Vibration mode suppression across frequency in the case of single-mode harmonic motion, for each of the three motors (a)-(c).

FIGURE 16 .
FIGURE 16.CFAR detection performance: average false discovery rate (FDR) vs. CFAR threshold for all motion and target velocity combinations.

TABLE 3 . Performance Summary and Comparison With Existing Motion Suppression Works
2 sin ω v,2 t + φ 0,2 ).Meanwhile, the transient component, ξ (t ), can be analyzed in the frequency domain by first rewriting as follows:As with the multimode harmonic motion case, the complex exponential can be rewritten using a Bessel function series expansion.Meanwhile, a rectangular window function translates to a sinc function in the Fourier domain, yielding: τp (t−td,p) , t − t d,p < τ p , 1, t − t d,p ≥ τ p .