1 Introduction↩︎

Diffusion models [1]–[3] have emerged as a powerful class of generative frameworks for time series modeling, owing to their strong ability to approximate complex data distributions. By progressively denoising latent representations initialized from Gaussian noise, these models learn a reverse process that captures rich temporal dependencies. As a result, diffusion-based approaches have achieved state-of-the-art performance across a range of time series tasks, including imputation [4], [5], forecasting [6], [7], and generation [8], [9].

Despite their success, existing diffusion models for time series still suffer from several notable limitations: (i) Insufficient reconstruction of high-frequency details [10]–[12]. In standard time-domain diffusion, each diffusion step injects isotropic Gaussian noise whose energy is uniformly distributed across frequencies. However, real-world time series usually exhibit highly imbalanced spectral energy distributions, with most energy concentrated in low-frequency components. This mismatch makes high-frequency reconstruction more challenging and limits the model’s ability to recover fine-grained temporal details. (ii) Imbalanced modeling of global trends and local dynamics [13]–[15]. Global trends and local dynamics often compete during reconstruction. Prioritizing long-term structures may lead to over-smoothed outputs, while emphasizing short-term fluctuations can disrupt global consistency. (iii) Frequency-insensitive modeling [16], [17]. Existing diffusion models typically use frequency-agnostic step embeddings and inject noise uniformly across frequency bands. As a consequence, they have limited ability to adapt their denoising behavior to the distinct recovery characteristics of low- and high-frequency components, which can lead to over-smoothed reconstructions when local variations are weak or sparsely observed.

Fig. 1 provides a representative imputation example on PhysioNet. The time-domain diffusion baseline, CSDI, captures the coarse level of the sequence but attenuates the sharp local variation in the highlighted region, resulting in an over-smoothed reconstruction. The spectrum computed over the same window further shows that the reconstructed signal underestimates middle- and high-frequency amplitudes compared with the ground truth. This example indicates that ignoring frequency dynamics can lead to observable reconstruction errors, motivating diffusion models that explicitly account for frequency-domain denoising.

In summary, these limitations can be attributed to two fundamental causes. (i) Time-domain-only noise modeling. In most existing diffusion models, both noise injection and denoising are performed exclusively in the time domain. This design overlooks the imbalanced spectral energy distribution of real-world time series, where high-frequency components often have weaker energy and are therefore harder to recover. Incorporating frequency-domain modeling enables explicit characterization of band-wise spectral statistics, allowing high-frequency transients and local dynamics to be more effectively captured during denoising. (ii) Lack of frequency-specific attention in step embeddings. Existing diffusion step embeddings do not incorporate frequency-dependent information and therefore fail to capture how different frequency bands should be emphasized or suppressed at different denoising stages. This prevents the model from adaptively allocating capacity between low- and high-frequency components, exacerbating the loss of fine-grained details and weakening multi-scale representation.

To address these challenges, we propose HyFAD, a Hybrid time-frequency Diffusion model with Frequency-Aware embedding for time series imputation. Unlike conventional time-domain diffusion approaches, HyFAD incorporates frequency-domain dynamics into the diffusion process, enabling explicit differentiation across spectral bands and improving high-frequency reconstruction. Specifically, HyFAD formulates a coupled time-frequency diffusion process, where frequency-domain and time-domain noise are injected sequentially in the forward process and removed in a mirrored order during reverse denoising. Furthermore, we analyze the relationship between diffusion steps and spectral components and introduce a frequency-aware step embedding that provides step-dependent spectral guidance, allowing the model to focus on appropriate frequency bands throughout denoising.

Our contributions are summarized as follows:

1. We propose HyFAD, a coupled time-frequency diffusion framework for time series imputation. The proposed forward and reverse processes explicitly model both time-domain and frequency-domain noise within a unified DDPM framework, enabling coarse-to-fine reconstruction.

2. We introduce a frequency-aware diffusion step embedding that provides step-dependent spectral guidance by modeling the correspondence between diffusion steps and frequency bands, allowing adaptive focus on appropriate spectral components during denoising.

3. We conduct extensive experiments on multiple real-world benchmark datasets, where HyFAD consistently achieves state-of-the-art performance, demonstrating the effectiveness of the proposed method.

2 Preliminaries↩︎

2.1 Discrete Fourier Transform↩︎

Given a real-valued discrete time series \(\mathbf{x} = [x_0, \ldots, x_{L-1}] \in \mathbb{R}^L\) with length \(L\), its unitary Discrete Fourier Transform (DFT), denoted by \(\tilde{\mathbf{x}} \in \mathbb{C}^L\), is defined as: \[\tilde{\mathbf{x}}_k = \frac{1}{\sqrt{L}}\sum_{\tau=0}^{L-1}\mathbf{x}_{\tau}\exp\left(- \frac{k2\pi i}{L}\tau\right),\] where \(k\) is the frequency index, \(\tau\) is the time index, and \(i\) denotes the imaginary unit. As \(\mathbf{x}\) is real-valued, \(\tilde{\mathbf{x}}\) exhibits conjugate symmetry i.e., \(\tilde{\mathbf{x}}_k = \tilde{\mathbf{x}}_{L-k}^*\), making almost half of the coefficients redundant. To obtain a compact representation that preserves isometry (i.e., satisfying Parseval’s identity) in the real domain, we define an orthonormal real DFT (rDFT) transform.

Specifically, we exploit the spectral symmetry to discard the redundant negative frequencies. To offset the energy reduction caused by this truncation, we apply a scaling factor of \(\sqrt{2}\) to the real and imaginary components of the intermediate frequencies. The elements of the transformed vector \(\mathbf{y}\) are arranged as follows: \[\mathbf{y}_m = \begin{cases} \Re(\tilde{\mathbf{x}}_0)&m=0\\ \Re(\tilde{\mathbf{x}}_m)&m=\frac{L}{2}~(\text{L is even})\\ \sqrt{2}\Re(\tilde{\mathbf{x}}_m)&0<m<\lceil\frac{L}{2}\rceil\\ \sqrt{2}\Im(\tilde{\mathbf{x}}_{k(m)})&m\ge N_r, \end{cases}\] where \(\Re(\tilde{\mathbf{x}})\) and \(\Im(\tilde{\mathbf{x}})\) denote the real and imaginary parts of \(\mathbf{x}\) and \(\mathbf{y} = [y_0,\ldots,y_{L-1}]^\top \in \mathbb{R}^L\) is the real-valued frequency representation constructed from \(\Re(\tilde{\mathbf{x}})\) and \(\Im(\tilde{\mathbf{x}})\). The index \(m=0,\ldots,L-1\) denotes the entry index of \(\mathbf{y}\), \(k(m)=m-N_r+1\) maps entries in the imaginary-part block of \(\mathbf{y}\) to the corresponding positive-frequency index and \(N_r = \lfloor\frac{L}{2}\rfloor\) is the number of frequency bins in the complex spectrum.

It is worth noticing that rDFT is a linear transformation, i.e., \(\mathbf{y} = \mathbf{W}\mathbf{x}\), where \(\mathbf{W} \in \mathbb{R}^{L \times L}\) is the transformation matrix and orthogonal (i.e., \(\mathbf{W}\mathbf{W}^T = \mathbf{I}\)). This implies that the rDFT operation is invertible, allowing for reconstructing the time-domain representation \(\mathbf{x}\) from the frequency-domain representation \(\mathbf{y}\) with \(\mathbf{x}=\mathbf{W}^T\mathbf{y}\).

2.2 Time Series Diffusion in the Frequency Domain↩︎

While most existing diffusion-based approaches [18], [19] for time series focus exclusively on modeling in the time domain, [13] explores the diffusion process in the frequency domain by introducing a noise scaling matrix \(\Lambda\) to the frequency-domain diffusion process under DFT. \(\Lambda\) plays to essential roles: (i) Rescaling spectral coefficients to make frequency-domain noise isotropic across all frequency components, (ii) Ensuring energy consistency between the time and frequency domains under the frequency domain transform. For an input sequence \(\mathbf{x}\in\mathbb{R}^{K\times L}\), the diagonal noise scaling matrix \(\Lambda\in\mathbb{R}^{L\times L}\) is defined as : \[\Lambda=\begin{cases} 1& n=1 \text{ or n=\frac{L}{2}}~(L~\text{is even}) \\ \frac{1}{\sqrt{2}}&\text{otherwise} \end{cases} \label{eq:g}\tag{1}\] Different from [13], we set \(\Lambda=\mathbf{I}\) since rDFT in Sec.2.1 is employed in our implementations instead of DFT for isotropic noise and energy consistency.

2.3 Problem Definition and Notations↩︎

Problem Statement (Time Series Imputation). Given a time series with missing value \(\tilde{\mathbf{X}}=(\mathbf{X}, \mathbf{M}, \mathbf{T})\), the goal of time series imputation is to learn an imputation function \(\mathcal{M}_{\theta}\), such that \[\bar{\mathbf{X}} = \mathcal{M}_{\theta}(\tilde{\mathbf{X}}),\] where \(\bar{\mathbf{X}}\in \mathbb{R}^{K \times L}\) is the imputed time series, \(\bar{\mathbf{X}}_{i,j}\) denotes the imputation output if \(\mathbf{M}_{i,j} = 1\), otherwise \(\bar{\mathbf{X}}_{i,j} = \tilde{{\mathbf{X}}}_{i,j}\).

Notations. We adopt the following notations. We use superscripts to denote the domain of the data, and subscripts to denote denoising steps, e.g., \(\mathbf{x}_k^f\) refers to the data at diffusion step \(k\) in the frequency domain and \(\boldsymbol{\epsilon}_k^t\) refers to the noise at diffusion step \(k\) in the time domain. \(\beta\) is a predefined scheduler over \(T\) denoising steps and \(\alpha_{k} = 1-\beta_k,k\in\{1,2,\cdots,T\}\). We denote \(\bar{\alpha}_{i:j}=\Pi_{k=i}^j\alpha_k\) and if \(i=1\), we omit the lower index \(i\) and use \(\bar{\alpha}_j\) for short. Specially, we define \(\bar{\alpha}_{0:j}=1\) for arbitrary \(j\). Besides, for time series data, we following the same notations as CSDI, i.e., \(\mathbf{X}^{\text{ta}},\mathbf{X}^{\text{obs}}\) denotes the imputation targets and ground truth values and \(\mathbf{M}^{\text{cond}}\) denotes the condition matrix.

3 Methodology↩︎

3.1 Hybrid Time-Frequency Diffusion Model↩︎

3.1.1 Forward Process of HyFAD↩︎

The forward process of HyFAD consists of two successive parts in a single step, i.e., noise is first injected in the frequency domain and subsequently in the time domain. By iteratively performing these two successive noise injection procedures, the original data distribution \(p(x_0)\) is transformed into a target Gaussian distribution. Given frequency transform \(\mathcal{F}\), we first map the data into corresponding domain via \(\mathcal{F}\), and then inject random gaussian noise \(\boldsymbol{\epsilon}_k^f\sim\mathcal{N}(0,\mathbf{I})\) into this frequency domain representation following a DDPM-style forward diffusion process, which is formulated as: \[\mathbf{x}_k^f = \sqrt{\alpha_k^f}\mathcal{F}(\mathbf{x}_{k-1}^t) + \sqrt{\beta_k^f}\sqrt{1-\lambda}(\Lambda\boldsymbol{\epsilon}_k^f), \label{eq:addnoisef1}\tag{2}\] where \(\lambda\) is the coefficient for balancing variance between frequency and time domain noise and \(\boldsymbol{\epsilon}_k^f\sim\mathcal{N}(0,\mathbf{I})\). Then the frequency-domain noisy sample \(\mathbf{x}_{k}^f\) in Eq.@eq:eq:addnoisef1 is transformed back to time domain: \[\mathcal{F}^{-1}(\mathbf{x}_k^f) = \sqrt{\alpha_k^f}\mathbf{x}_{k-1}^t + \sqrt{\beta_k^f}\sqrt{1-\lambda}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}_k^f), \label{eq:addnoisef2}\tag{3}\] where \(\boldsymbol{\epsilon}_k^f\sim\mathcal{N}(0,\mathbf{I})\) and we can obtain the posterior \(q(\mathcal{F}^{-1}(\mathbf{x}_{k}^f)\vert \mathbf{x}_{k-1}^t)\mathrel{\vcenter{:}}=\mathcal{N}(\sqrt{\alpha_k^f}\mathbf{x}_{k-1}^t,\beta_k^f(1-\lambda)\mathcal{F}^{-1}\Lambda\Lambda^T(\mathcal{F}^{-1})^{T})\).

Then another independent gaussian noise \(\boldsymbol{\epsilon}_k^t\sim\mathcal{N}(0,\mathbf{I})\) is injected to \(\mathcal{F}^{-1}(\mathbf{x}_k^f)\), i.e., noise is injected in the time domain: \[\mathbf{x}_k^t = \sqrt{\alpha_k^t}\mathcal{F}^{-1}(\mathbf{x}_k^f) + \sqrt{\beta_k^t}\sqrt{\lambda}\boldsymbol{\epsilon}_k^t. \label{eq:addnoiset}\tag{4}\] Thus, the posterior can be written as \(q(\mathbf{x}_k^t\vert \mathcal{F}^{-1}(\mathbf{x}_{k}^f))\mathrel{\vcenter{:}}=\mathcal{N}(\sqrt{\alpha_k^t}\mathcal{F}^{-1}(\mathbf{x}_{k}^f),\beta_k^t\lambda\mathbf{I})\). With this, the noise injection loop \(k\) is completed. We can formulate the relationship between \(\mathbf{x}_{k}^t\) and \(\mathbf{x}_{k-1}^t\) as: \[\mathbf{x}_k^t = \sqrt{\alpha_k^t\alpha_k^f}\mathbf{x}_{k-1}^t + \sqrt{\alpha_k^t(1-\alpha_k^f)}\sqrt{1-\lambda}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}_k^f)+\sqrt{1-\alpha_k^t}\sqrt{\lambda}\boldsymbol{\epsilon}_k^t. \label{eq:addnoiseall}\tag{5}\] By iteratively preforming the noise injection process according to Eq.@eq:eq:addnoisef2 , Eq.@eq:eq:addnoiset and Eq.@eq:eq:addnoiseall , the closed-form relationship between \(\mathbf{x}_k^t\) and \(\mathbf{x}_0^t\) can be derived as: \[\mathbf{x}_k^t = \sqrt{\bar{\alpha}_k^t\bar{\alpha}_k^f}\mathbf{x}_0^t + \sqrt{1-\lambda}\sum_{s=1}^k\sqrt{\beta_{s}^f}\sqrt{\frac{\bar{\alpha}_k^t\bar{\alpha}_k^f}{\bar{\alpha}_{s-1}^t\bar{\alpha}_s^f}}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}_s^f) + \sqrt{\lambda}\sum_{s=1}^k\sqrt{\beta_s^t}\sqrt{\frac{\bar{\alpha}_k^t}{\bar{\alpha}_s^t}}\boldsymbol{\epsilon}_s^t. \label{eq:addnoisetotal}\tag{6}\] According to Eq.@eq:eq:addnoisetotal , for any diffusion step \(k\), \(\mathbf{x}_k^t\) can be expressed as a linear transformation of the original signal \(\mathbf{x}_0^t\) with two groups of independent gaussian noise \(\boldsymbol{\epsilon}_s^t\) and \(\boldsymbol{\epsilon}_s^f\), \(s\in\{1,2,\cdots,k\}\). Therefore, the posterior \(q(\mathbf{x}_k^t\vert\mathbf{x}_0^t)\) remains gaussian, i.e., \[q(\mathbf{x}_k^t\vert\mathbf{x}_0^t) \mathrel{\vcenter{:}}= \mathcal{N}(\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k), \label{eq:posterior}\tag{7}\] where \(\boldsymbol{\mu}_k = \sqrt{\bar{\alpha}_k^t\bar{\alpha}_k^f}\mathbf{x}_0^t \label{eq:postmean}\) and \(\boldsymbol{\Sigma}_k = (1-\lambda)\sum_{s=1}^k\beta_s^f\frac{\bar{\alpha}_k^t}{\bar{\alpha}_{s-1}^t}\frac{\bar{\alpha}_k^f}{\bar{\alpha}_s^f}\left(\mathcal{F}^{-1}\Lambda\Lambda^T(\mathcal{F}^{-1})^T\right) +\lambda\sum_{s=1}^k\beta_s^t\frac{\bar{\alpha}_k^t}{\bar{\alpha}_s^t}\mathbf{I}. \label{eq:postvariance}\)

3.1.2 Reverse Process of HyFAD↩︎

To progressively reconstruct the original data from the diffused samples, we construct a reverse Markov chain that is consistent with the forward noise injection process. Following the DDPM formulation, we parameterize the reverse distribution \(p_{\theta}(\mathbf{x}_{k-1}^t\vert\mathbf{x}_{k}^t)\) using a Gaussian approximation, aiming to model the true posterior \(q(\mathbf{x}_{k-1}^t\vert \mathbf{x}_{k}^t,\mathbf{x}_0^t)\). This enables closed-form stepwise sampling throughout the denoising trajectory. Since our method introduces a combination of time-domain noise and frequency-domain noise (Eq.@eq:eq:addnoiseall ) at each diffusion step, the reverse process must explicitly account for the structure of this composite noise to ensure faithful reconstruction of original data.

To improve sampling stability, we adopt a DDIM-style sampling strategy [20], which follows a deterministic trajectory in the reverse process. We initialize from a Gaussian prior that is consistent with the terminal distribution of the forward diffusion. As in Eq.@eq:eq:posterior , the sampling process starts from \[\mathbf{x}_{T}^t=\sqrt{\lambda}\boldsymbol{\epsilon}^t+\sqrt{1-\lambda}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}^f),\boldsymbol{\epsilon}^t,\boldsymbol{\epsilon}^f\sim\mathcal{N}(0,\mathbf{I}).\] At each step \(k\), the model predicts the time domain noise \(\hat{\boldsymbol{\epsilon}}_{k}^t\) and frequency domain noise \(\hat{\boldsymbol{\epsilon}}_{k}^f\). Reverse denoising is then performed in a step-wise, coupled manner within the same diffusion loop. First, the time-domain noise in \(\mathbf{x}_{k}^t\) is removed according to: \[\mathcal{F}^{-1}(\mathbf{x}_{k-1}^f) = \frac{1}{\sqrt{\alpha_{k}^t}}\left(\mathbf{x}_{k}^t-\frac{\beta_{k}^t}{\sqrt{1-\bar{\alpha}_{k}^t}}\hat{\boldsymbol{\epsilon}}_{k}^t\right). \label{eq:timedenoising}\tag{8}\] Then the frequency-domain noise in \(\mathcal{F}^{-1}(\mathbf{x}_{k-1}^f)\) is removed for \(\mathbf{x}_{k-1}^t\): \[\mathbf{x}_{k-1}^t = \frac{1}{\sqrt{\alpha_{k}^f}}\left(\mathcal{F}^{-1}(\mathbf{x}_{k-1}^f)-\sqrt{1-\lambda}\sqrt{\beta_{k}^f}\mathcal{F}^{-1}(\Lambda\boldsymbol{\hat{\epsilon}}_k^f)\right). \label{eq:freqdenosing}\tag{9}\]

Eq.@eq:eq:timedenoising and Eq.@eq:eq:freqdenosing define distinct update rules for time and frequency domain denoising process. To summarize, the step-wise denoising strategy i.e., time-to-frequency denoising, serves as the mirror counterpart of the forward diffusion process (frequency-to-time noise injection). Empirically, the time-domain branch is effective in stabilizing global structures and low-frequency trends, while the frequency-domain branch focuses on high-frequency component, together forming a coarse-to-fine denoising mechanism. The forward and reverse process is illustrated in Fig.2.

3.1.3 Loss Function↩︎

As stated in Sec.3.1.1 and Sec.3.1.2, the overall noise injection process consists of two successive sub-processes, i.e., time and frequency domain noise injection. We adopt the commonly used noise estimation approach in diffusion models. The loss function consists of three components: (i) time domain noise estimation loss, i.e., the \(l_2\) distance between the estimated time domain noise and forward time domain noise, (ii) frequency domain noise estimation loss, i.e., the \(l_2\) distance between the estimated frequency domain noise and forward frequency domain noise. (iii) consistency loss, i.e., the \(l_2\) distance between total noise from step \(k\) to step \(k-1\). This consistency term mitigates potential mismatches arising from estimating the two branches separately and discourages arbitrary residual reallocation across branches. Accordingly, the overall loss is defined as: \[\mathcal{L}^{d} = \mathbb{E}_{\mathbf{x}_0\sim q(\mathbf{x}_0),k}(\Vert \boldsymbol{\epsilon}_k^t-\boldsymbol{\hat{\epsilon}}_k^t\Vert_2^2 +\Vert \boldsymbol{\epsilon}_k^f-\mathcal{F}^{-1}(\Lambda\boldsymbol{\hat{\epsilon}}_k^f)\Vert_2^2 +\omega\Vert(\boldsymbol{\epsilon}_k^t+\boldsymbol{\epsilon}_k^f)-(\hat{\boldsymbol{\epsilon}}_k^t+\mathcal{F}^{-1}(\Lambda\hat{\boldsymbol{\epsilon}}_k^f))\Vert_2^2), \label{eq:loss}\tag{10}\] where \(\omega\) is the weight for consistency loss. In addition, since our model is specifically designed for time series imputation task, we focus solely on the reconstruction error over masked regions during training, i.e., \(\mathcal{L}= \mathbf{M}^t\odot\mathcal{L}^{d}\). The detailed training and sampling algorithm is presented in Sec.7.7.

3.2 Frequency-aware Diffusion Embedding↩︎

Considering the relationship between denoising steps and the frequency components of the data, we present the following proposition [21]–[23]:

Motivated by Prop. 1, we propose a frequency-aware diffusion embedding to capture the stage-dependent recoverability of different frequency components during denoising. As the noise level changes across diffusion steps, different frequency bands should be emphasized differently. To this end, the proposed embedding integrates three components: noise-aware frequency gate, stage-aware frequency schedule, and missing-aware band reliability.

Noise-aware Frequency Gate. Given the accumulated diffusion noise, whose expected spectral energy satisfies \(\mathbb{E}[\vert \hat{\epsilon}_t(\omega)\vert^2] = \int_0^t \vert g(s)\vert^2\mathrm{d}s\), the diffusion step \(t\) corresponds to a frequency-agnostic noise scale \(N_t=\int_0^t \vert g(s)\vert^2\mathrm{d}s\), which decreases monotonically in the reverse process. Since different bands have different signal-to-noise ratios, treating all frequency components uniformly is suboptimal. We therefore introduce a noise-aware gate, with calibration parameters \(\gamma\) and \(s\) controlling the gating threshold and its scale: \[\hat{g}_b(t) = \sigma\left(\frac{\log P_b^{\mathrm{sig}} -\log (\gamma N_ts)}{\tau}\right),\] where \(P_b^{\mathrm{sig}}\) is the spectral energy of frequency band \(b\) computed from the signal proxy detailed in Sec. 3.3.1, and \(\tau\) controls the sharpness of the gate. To avoid numerical issues, we set a minimum value \(g_{\min}\): \[g_b(t) = g_{\min} + (1-g_{\min})\hat{g}_{b}(t).\]

Stage-aware Frequency Schedule. Following Prop. 1, we further encode the progressive low-to-high frequency recovery process into the diffusion embedding. For each band \(b\), the stage-dependent scheduling weight \(S_b(t)\) measures its relative importance at diffusion step \(t\): \[S_b(t) = \exp\left(-\left(\frac{\omega_b}{c(t)}\right)^p\right),\] where \(\omega_b\) is the normalized frequency coordinate, \(c(t)\) gradually relaxes the effective bandwidth as denoising progresses, and \(p>0\) is a constant. Unlike the noise-aware gate, \(S_b(t)\) depends only on \(t\) and \(\omega_b\), providing a stage-dependent frequency prior independent of actual spectral energy or noise magnitude.

Band Reliability. For time series imputation, missing values can introduce spectral artifacts and weaken the reliability of frequency-domain modeling. We therefore define a missing-aware reliability weight for each frequency band: \[R_b = \frac{P_b^{\mathrm{sig}}}{P_b^{\mathrm{sig}}+\kappa P_b^{\mathrm{mask}}},\] where \(P_b^{\mathrm{mask}}\) denotes the spectral energy induced by the missing pattern in band \(b\), and \(\kappa\) is a constant.

Finally, we combine the noise-aware gate, stage-aware schedule, and band reliability through element-wise multiplication, and project the result into the embedding space: \[\mathbf{g}(t) = \mathrm{Proj}\left(g(t)\odot S(t)\odot R,\, \frac{d}{2}\right), \label{eq:proj}\tag{11}\] where \(\mathrm{Proj}(\cdot,\frac{d}{2})\) resamples \(\lfloor\frac{L}{2}\rfloor + 1\) frequency bins to \(\frac{d}{2}\) embedding dimensions. The resulting gating vector modulates sinusoidal basis functions to form the frequency-aware diffusion embedding: \[\mathbf{e}(t) = \left[ \sin(t\,\mathbf{f}_d)\odot \mathbf{g}(t), \; \cos(t\,\mathbf{f}_d)\odot \mathbf{g}(t) \right], \label{eq:grid}\tag{12}\] where \(\mathbf{f}_d\) denotes a linearly spaced frequency grid and \(t\) is the diffusion step. Details are provided in Appendix 7.5.

3.3 Model Details↩︎

3.3.1 Missing-aware Frequency Modeling↩︎

Direct frequency transformation of the masked series \(\tilde{\mathbf{X}}=\mathbf{M}\odot \mathbf{X}\) often induces spectral leakage, distorting the true spectral morphology. To mitigate this, we introduce a Signal Proxy mechanism that serves as a diffusion-synchronized spectral reference. Designed to ensure observational consistency and reliable estimation, the proxy evolves with the diffusion process to reflect the changing signal-to-noise ratio. Formally, we define the signal proxy \(\tilde{\mathbf{x}}_{\mathrm{proxy}}\) at step \(t\) as: \[\tilde{\mathbf{x}}_{\mathrm{proxy}} = \mathbf{M}\odot \mathbf{x}_{\mathrm{obs}} + (1-\mathbf{M})\odot \mathbf{x}_t\] where \(\mathbf{x}_{\mathrm{obs}}\) denotes the observed values and \(\mathbf{x}_t\) is the temporal latent variable. This formulation preserves the authentic signal structure at observed locations while populating missing intervals with state estimates consistent with the current diffusion phase. This diffusion-aware mechanism effectively minimizes spectral interference, enabling robust frequency representation learning.

3.3.2 Model Architecture↩︎

We present the details of the denoising network.We adopt a cascaded dual-branch network as the core denoising architecture, where both branches share the same network structure. Considering the intra-channel and inter-channel dependencies in multivariate time series, the temporal-feature transformer architecture from CSDI [4] is adopted as the diffusion backbone.

Specifically, our model takes the noisy data, conditional information and step embedding as input. The input is first processed by a temporal attention module and a feature attention module for feature extraction, followed by a linear projection to estimate time-domain noise \(\hat{\boldsymbol{\epsilon}}_{t}\) To ensure that the frequency-domain branch focuses exclusively on estimating frequency-domain noise, we transform the output of the time-domain branch into the frequency domain via a given transformation \(\mathcal{F}\) and use it as the input to the frequency-domain branch. The frequency-domain branch is symmetric to the time-domain branch, except that the step embedding is replaced with the frequency-aware diffusion embedding introduced in Sec. 3.2. Finally, the frequency-domain output is mapped back to the time domain via \(\mathcal{F}^{-1}\). Overall, the model produces two outputs: the estimated time-domain noise \(\hat{\boldsymbol{\epsilon}}_{t}\) and frequency-domain noise \(\hat{\boldsymbol{\epsilon}}_{f}\). The details of our model architecture is in the Appendix 7.6.

4 Experiments↩︎

Experimental Settings. All experiments are implemented in Python 3.12 with Pytorch [24] on a single Nvidia RTX 3090 GPU. We use Adam [25] with an initial learning rate of \(1\times10^{-3}\) and a multi-step scheduler that decays the learning rate by 10 at 75% and 90% of the total epochs. Quadratic noise schedules are adopted for both time- and frequency-domain diffusion, and the frequency transformation \(\mathcal{F}\) is set to rDFT as described in Sec. 2.1.

Datasets and Evaluation Metrics. We evaluate HyFAD on two widely used imputation benchmarks: PhysioNet [26] and Air Quality [27]. We report MAE and RMSE between the imputed values and ground truth at missing positions. More details are provided in Appendix 7.8.

Baselines. We compare HyFAD with classical methods (Mean and Lerp), deep learning approaches (BRITS [28], SAITS [29], TimesNet [30]), deep generative models (GP-VAE [31], SSGAN [32], CSDI [4], LSCD [33]), and a time-series foundation model (ModernTCN [34]).

4.1 Time Series Imputation Results↩︎

We evaluate HyFAD under different missing-data scenarios on two real-world datasets, i.e., PhysioNet and Air Quality. Table 1 reports the quantitative results, with visualizations provided in Appendix 7.13.

As shown in Table 1, HyFAD achieves the best performance across all missing rates on PhysioNet. Compared with the second-best baseline, it reduces MAE by 2.8%, 5.6%, and 8.9% under 10%, 50%, and 90% missingness, respectively, while also consistently lowering RMSE. On Air Quality, HyFAD also obtains the best MAE and RMSE, demonstrating its effectiveness across different real-world time-series scenarios. Please refer to Appendix.7.13.3 for more detailed comparison.

4.2 Ablation Studies and Parameter Analysis↩︎

To validate the effectiveness of the proposed frequency-aware diffusion embedding, we replace it with the diffusion embedding in [4]. We further analyze several key hyperparameters, including the terminal values of the noise schedules in time- and frequency-domain diffusion, i.e., \(\beta_{\mathrm{end}}^t\) and \(\beta_{\mathrm{end}}^f\), the balancing parameter \(\lambda\) that controls the relative strengths of the two diffusion processes, and \(\omega\), the weight of the consistency loss.

As shown in Table 2, the proposed Frequency-Aware Diffusion Embedding consistently outperforms the standard diffusion embedding across all missing rates on PhysioNet. In addition, since CSDI can be viewed as the time-domain-only counterpart of HyFAD, the consistent improvements over CSDI in Table 1 further validate the necessity of hybrid time-frequency modeling. These results support the core design of HyFAD: the time branch preserves global structure, while the frequency branch complements local dynamics and high-frequency details.

For noise allocation, the best performance is achieved when the temporal process is dominant (\(\lambda=0.75\)). The performance drops at \(\lambda=0.25\) and \(0.5\) indicate that frequency-domain diffusion is more suitable for refining local dynamics that are not fully captured in the time domain, rather than serving as the primary noise carrier. For diffusion schedules, the optimal temporal noise level appears at \(\beta_\mathrm{end}^t=0.5\), while either excessive (\(0.75\)) or insufficient (\(0.25\)) noise weakens the model’s ability to learn the underlying data distribution. Regarding \(\beta_{\mathrm{end}}^f\), small values provide insufficient spectral information for refinement, whereas overly large values over-regularize the temporal branch and lead to degraded performance.

Finally, \(\omega=0.4\) achieves the best overall trade-off, yielding the lowest RMSE across all missing rates while maintaining competitive MAE. Smaller weights lead to insufficient coupling, whereas larger weights may slightly improve some MAE values but generally degrade RMSE, suggesting over-constrained cross-domain fitting.

5 Related Works↩︎

Time series imputation aims to recover complete sequences from partially observed data. Early methods such as mean and linear interpolation (Lerp) rely on simple statistical assumptions and often fail to capture complex temporal dynamics. Deep learning-based approaches, including BRITS [28], SAITS [29], and TimesNet [30], improve imputation by modeling complex temporal dependencies.

Diffusion models have recently been introduced to time series imputation, motivated by their success in other domains [6], [35], [36]. CSDI [4] formulates imputation as conditional score matching with observed values as conditions, while SSSD [5] incorporates state space models for temporal modeling. Recent works also explore frequency-domain information: FGTI [11] introduces frequency-space conditions for multi-domain representation learning, and FourierDiffusion [13] transforms signals into the spectral domain and conducts diffusion entirely in frequency space.

Frequency-domain modeling has also been widely studied in time-series analysis, as spectral characteristics provide complementary information beyond the time domain. Fedformer [37] combines decomposition with a frequency-domain Transformer and frequency mode selection for efficient long-term forecasting. FiLM [38] uses Legendre memory with Fourier- and low-rank-based frequency-enhanced layers to preserve dominant low-frequency patterns and suppress noise. TFAD [39] combines time- and frequency-domain analysis with sequence decomposition, context comparison, and data augmentation for anomaly detection, where the frequency branch mainly provides a discriminative representation space.

HyFAD differs from prior works by explicitly coupling time- and frequency-domain diffusion for time-series imputation. Rather than using frequency information only for representation learning, HyFAD performs forward and reverse diffusion directly in both domains. It further introduces frequency-aware diffusion embedding and step-aware spectral guidance to support progressive denoising and reconstruction, improving recovery of both global trends and fine-grained temporal variations.

6 Conclusion↩︎

In this paper, we propose HyFAD, a hybrid time-frequency diffusion model with frequency-aware diffusion embedding. By integrating time- and frequency-domain processes within a unified diffusion framework, HyFAD enables coarse-to-fine reconstruction. The proposed frequency-aware diffusion-step embedding adaptively emphasizes appropriate spectral components during denoising, enhancing reconstruction of high-frequency components. Experiments demonstrate consistently strong performance across multiple datasets and varying missing rates.

7 Appendix↩︎

7.1 Choice of \(\mathcal{F}\) and \(\Lambda\)↩︎

Notably, our framework imposes no strict restrictions on the transform \(\mathcal{F}\), as any invertible linear frequency-domain transformation satisfies our requirements. On the other hand, the \(\Lambda\) matrix is intrinsically coupled with the chosen transform \(\mathcal{F}\). Any modification to the transform \(\mathcal{F}\) necessitates a corresponding adjustment to the \(\Lambda\) matrix to ensure energy conservation between the time and frequency domains.

7.2 Details in Forward Process↩︎

Posterior of Eq.@eq:eq:addnoisef2 : Eq.@eq:eq:addnoisef2 presents \[\mathcal{F}^{-1}(\mathbf{x}_k^f) = \sqrt{\alpha_k^f}\mathbf{x}_{k-1}^t + \sqrt{\beta_k^f}\sqrt{1-\lambda}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}_k^f),\qquad\boldsymbol{\epsilon}_k^f\sim\mathcal{N}(0,\mathbf{I}).\]

\(\mathcal{F}^{-1}(\mathbf{x}_{k-1}^t)\) is a linear transformation of \(\mathbf{x}_k^t\) along with a Gaussian term, Therefore, the posterior \(q(\mathcal{F}^{-1}(\mathbf{x}_{k}^f)\vert \mathbf{x}_{k-1}^t)\) remains gaussian. The posterior mean is \(\sqrt{\alpha_k^f}\mathbf{x}_{k-1}^t\) and the standard deviation is \[(\sqrt{\beta_k^f}\sqrt{1-\lambda}\mathcal{F}^{-1}\Lambda)(\sqrt{\beta_k^f}\sqrt{1-\lambda}\mathcal{F}^{-1}\Lambda)^T=\beta_k^f(1-\lambda)\mathcal{F}^{-1}\Lambda\Lambda^T(\mathcal{F}^{-1})^T\] Here \(\mathcal{F}^{-1}\) and \((\mathcal{F}^{-1})^T\) denotes a matrix \(\mathbf{U}\) due to the linearity of transformation \(\mathcal{F}\), i.e., \(\mathcal{F}(\mathbf{x}) = \mathbf{Ux}\).

Posterior of Eq.@eq:eq:addnoiset : Eq.@eq:eq:addnoiset presents \(\mathbf{x}_k^t = \sqrt{\alpha_k^t}\mathcal{F}^{-1}(\mathbf{x}_k^f) + \sqrt{\beta_k^t}\sqrt{\lambda}\boldsymbol{\epsilon}_k^t\), \(\boldsymbol{\epsilon}_k^t\sim\mathcal{N}(0,\mathbf{I})\).\(\mathbf{x}_k^t\) is a linear transformation of \(\mathcal{F}^{-1}(\mathbf{x}_k^f)\) with a Gaussian term \(\sqrt{\beta_k^t}\sqrt{\lambda}\boldsymbol{\epsilon}_k^t\). Therefore, the posterior is still Gaussian with mean \(\sqrt{\alpha_{k}^t}\mathcal{F}^{-1}(\mathbf{x}_k^f)\) and standard deviation \(\beta_k^t\lambda\mathbf{I}\).

Proof of Eq.@eq:eq:addnoisetotal  (the relationship between \(\mathbf{x}_k^t\) and \(\mathbf{x}_k^0\) in the forward process):

From Eq.@eq:eq:addnoiseall : \[\mathbf{x}_k^t = \sqrt{\alpha_k^t\alpha_k^f}\mathbf{x}_{k-1}^t + \sqrt{\alpha_k^t(1-\alpha_k^f)}\sqrt{1-\lambda}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}_k^f)+\sqrt{1-\alpha_k^t}\sqrt{\lambda}\boldsymbol{\epsilon}_k^t \label{eq:appen1}\tag{13}\] we can get: \[\mathbf{x}_k^t = \sqrt{\bar{\alpha}_k^t\bar{\alpha}_k^f}\mathbf{x}_0^t + \sqrt{1-\lambda}\sum_{s=1}^k\sqrt{\beta_{s}^f}\sqrt{\frac{\bar{\alpha}_k^t}{\bar{\alpha}_{s-1}^t}}\sqrt{\frac{\bar{\alpha}_k^f}{\bar{\alpha}_s^f}}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}_s^f) + \sqrt{\lambda}\sum_{s=1}^k\sqrt{\beta_s^t}\sqrt{\frac{\bar{\alpha}_k^t}{\bar{\alpha}_s^t}}\boldsymbol{\epsilon}_s^t \label{eq:appen2}\tag{14}\]

Posterior of Eq.@eq:eq:addnoisetotal : Eq.@eq:eq:addnoisetotal indicates \(\mathbf{x}_k^t\) is linear transformation of \(\mathbf{x}_0^t\) with two groups of independent gaussian noise. Therefore, \(q(\mathbf{x}_k^t\vert\mathbf{x}_0^t)\) is still gaussian with the standard deviation of the sum of two the two groups of gaussian noise.

7.3 Details in the reverse process↩︎

Noise prior in the reverse process. At the end of the forward process, \(\bar{\alpha}_k^t,\bar{\alpha}_k^f\to 0\), therefore, the mean of \(\mathbf{x}_k^t\) is \(0\). For the standard deviation term, it is the linear combination of two independent gaussian noises, so the reverse process starts from \(\mathbf{x}_{T}=\sqrt{\lambda}\boldsymbol{\epsilon}^t+\sqrt{1-\lambda}\mathcal{F}^{-1}(\Lambda\boldsymbol{\epsilon}^f),\boldsymbol{\epsilon}^t,\boldsymbol{\epsilon}^f\sim\mathcal{N}(0,\mathbf{I})\).

7.4 Proof of Proposition 1↩︎

7.5 Details of Frequency-aware Diffusion Embedding↩︎

Factor \(s\) in Noise-aware Frequency Gate. Due to substantial variations in spectral energy scales across different datasets and variables, we introduce a global calibration factor \(s\) to align the magnitude of band-wise spectral energy, which ensures the thresholding operation in the frequency-aware diffusion embedding is performed on a consistent relative scale.

During training, the factor \(s\) is updated per batch using a global exponential moving average (EMA). Specifically, given the band-wise spectral energy \(P_b^s\) of the current batch, \(s_0\) is initialized as the median of \(P_b^s\): \[s_0 = \mathrm{median}(P_b^s),\] \(s\) is then updated as: \[s_k \leftarrow \alpha s_{k-1} + (1-\alpha)\hat{s}_{k},\] where \(\alpha \in (0,1)\) denotes the EMA decay rate and we set \(\alpha=0.99\) in our implementation, \(s_{k-1}\) is the factor after batch \(k-1\) and \(\hat{s}_k\) is the observations from batch \(k\), i.e., \(\hat{s}_k=\mathrm{median}(P_k^s)\). During inference, the factor \(s\) is fixed to the value accumulated during training and is no longer updated. This design enforces a strict separation between training and testing statistics, thereby preventing potential data leakage and ensuring consistent gating behavior across test samples.

\(\omega_b\) and \(c(t)\) in Stage-aware Frequency Schedule. Let \(b\in\{0,1,\cdots,B-1\}\) index the \(B=\lfloor \frac{L}{2}\rfloor + 1\) frequency bands. The relative position of the \(b\)-th band is: \[\tilde{\omega}_b = \frac{b}{B-1},\] To align the relative frequency coordinate with the numerical scale of the stage-dependent cutoff, we further multiply a max frequency \(\omega_{\max}\): \[\omega_b = \tilde{\omega}_b\cdot \omega_{\max}\]

The function \(c(t)\) is a stage-dependent frequency scaling function that controls the effective frequency range emphasized by the model at diffusion step \(t\), which is defined as: \[c(t) = c_{\min} + (c_{\max} - c_{\min})\left(1-\frac{t}{T}\right)^q, \label{eq:ct}\tag{19}\] where \(t\) is the current diffusion step and \(T\) is the total number of diffusion steps. \(c_{\min}\) and \(c_{\max}=\omega_{\max}\) are the minimum and maximum cutoff frequency, \(q>0\) is a transition exponent. \(c(t)\) remains close to \(c_{\min}\) during the early stages of the reverse diffusion process, thereby emphasizing low-frequency components. As the denoising process progresses, \(c(t)\) gradually approaches \(c_{\max}\), allowing the model to progressively incorporate higher-frequency details.

Embedding Details. The projection operator \(\mathrm{Proj}(\cdot,\tfrac{d}{2})\) in Eq.@eq:eq:proj is implemented via one-dimensional linear interpolation along the frequency axis using torch.nn.functional.interpolate (mode=’linear’), which resamples the frequency-wise gating vector from \(\lfloor L/2\rfloor+1\) frequency bins to \(\tfrac{d}{2}\) embedding dimensions.

The frequency grid \(\mathbf{f}_d= [f_1,f_2,\cdots,f_{\frac{d}{2}}]\in \mathbb{R}^{\frac{d}{2}}\) in Eq.@eq:eq:grid is constructed by uniform sampling from \([0,f_{\max}]\): \[f_i = \frac{i-1}{\frac{d}{2}-1}f_{\max}\]

7.6 Architecture Details↩︎

Fig.3 presents the diffusion backbone of our denoising block. The block adopts a time-frequency coupled dual-branch residual The input sequence is first projected using a \(1\times1\) convolution followed by a ReLU activation. Time and Frequency diffusion step embeddings are generated via fully connected layers with SiLU activations and then concatenated with external side information.

During residual modeling, the network first processes the representation in the time-domain branch. After an initial \(1\times1\) convolution and ReLU projection, the features are passed through a temporal Transformer layer and a feature Transformer layer to capture inter-channel and intra-channel dependencies. The outputs are fused using a gated activation unit (GAU) and further refined by stacked \(1\times1\) convolutions. Residual connections feed the results back into the input stream to form the input for the next layer.

Then the time-domain representation is then fed into the frequency-domain branch with a similar structure to get the frequency domain output,. By processing in a time-to-frequency order, the architecture establishes a complementary relationship between global trend modeling and fine-grained frequency refinement.

7.7 Training and Sampling Details↩︎

The detailed training and sampling algorithm is presented in Alg.4 and 5.

7.8 Dataset Details↩︎

The PhysioNet dataset [26] consists of 4000 irregularly-sampled medical time series data including 35 variables (e.g., Albumin and heart rate) collected from ICU with a total length of 48 hours and has been made available by their authors under the terms of the Creative Commons Attribution License 3.0 (CCAL). Consistent with previous studies [4], the dataset is processed hourly to get 48 timesteps. The processed dataset contains nearly 80% originally missing values without ground truth. In our experiments, we random select 10/50/90% of the observed values as the imputation targets (i.e., ground truth of test dataset).

The air quality dataset [40] contains PM2.5 data collected from 36 monitor stations in Beijing and is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license. All the air quality data are collected hourly for 12 months. The original dataset contains 13.3% missing values with artificial ground-truth with a non-random missing pattern.

The details of the datasets are presented in Tab.3.

7.9 Evaluation Metrics↩︎

In this section, we present the details of evaluation metrics in our experiments. \(y, \hat{y}\in\mathbb{R}^{K\times L}\) denote the ground truth and output of our model and \(M\) is the indicator matrix.

Mean Absolute Error (MAE) calculates the average \(L_1\) error between the imputed samples and the ground truth of the time series: \[\mathbf{MAE}(y,\hat{y}) = \frac{\sum_{i=1}^K\sum_{j=1}^LM_{ij}\vert y_{ij}-\hat{y}_{ij}\vert_1}{\sum_{i=1}^K\sum_{j=1}^LM_{ij}}\] MAE reflects the overall deviation across all points, emphasizes the overall accuracy of the model outputs. Compared with RMSE, MAE is more robust to outliers in the data.

Rooted Mean Square Error (RMSE) calculates the average \(L_2\) error between the imputed samples and the ground truth of time series: \[\mathbf{RMSE}(y,\hat{y}) = \sqrt{\frac{\sum_{i=1}^K\sum_{j=1}^LM_{ij}\Vert y_{ij}-\hat{y}_{ij}\Vert_2^2}{\sum_{i=1}^K\sum_{j=1}^LM_{ij}}}\] RMSE highlights potential large errors, so a few big deviations may dominate the score. Compared with MAE, RMSE is less robust to outliers but better captures worst-case performance.

7.10 Limitations↩︎

Here we provide an in-depth discussion regarding the limitations of our current work. First, our exploration of frequency-domain representations is primarily confined to the Discrete Fourier Transform (DFT). Future research could extend this framework to alternative transforms, such as the Discrete Cosine Transform (DCT) or Discrete Wavelet Transform (DWT), to better capture different kinds of features from different types of time series data. Second, this study focuses on deterministic time series imputation using DDIM for efficient sampling. A promising future direction involves adapting DDPM-based sampling schemes to our hybrid time-frequency architecture, as well as evaluating the model’s performance on probabilistic imputation tasks to better quantify uncertainty. Another limitation is that the frequency-domain sampling process can be less numerically stable than sampling in the time domain, particularly when reconstructing signals from sparse spectral representations. We intend to investigate specialized stabilization methods for frequency-domain diffusion in future works.

7.11 Implementation Details↩︎

Global factor \(s_k\) is updated as \(s_k = \alpha s_{k-1} + (1-\alpha)\hat{s}_k\), where \(\hat{s}_k\) is computed as the median spectral energy of the current training batch. It depends only on the observed values of the training samples and their diffusion states without any validation or test information. Moreover, this EMA is updated only during training and kept fixed at inference time, so it does not introduce data leakage.

The noise scheduler does not rely on validation data. The scheduler is predefined, and the noise scale \(N_t\) in the frequency branch (Eq. 18) is computed analytically from the diffusion coefficients rather than estimated from validation or test statistics. Therefore, neither the construction nor the use of the scheduler involves validation information.

All baselines are evaluated under the same setting. We strictly follow the CSDI setting for data preprocessing, splitting, and evaluation on both PhysioNet and Air Quality, ensuring a fair comparison across methods. For validation, we follow the same setting as CSDI: the validation set is only used to select the best checkpoint, and we do not use early stopping in our training.

7.12 Computational Cost↩︎

Compared with CSDI, the extra cost of HyFAD mainly comes from frequency-domain branch, time–frequency transforms and the joint optimization objective. Since the backbone is Transformer-based, the overall complexity is \(O(L^2)\) with transform cost \(O(L\log L)\). We also provide training and inference speed on PhysioNet (batch size = 16) in Table.4.

7.13 Experiment details↩︎

7.13.1 Hyperparameters↩︎

The hyperparameters in our experiments is detailed in Table.5 and the hyperparameters in frequency-aware diffusion embedding in presented in Table.6

7.13.2 Ablation on explicit frequency modeling↩︎

To empirically assess the contribution of explicit frequency modeling, we compare HyFAD with its pure time-domain counterpart. Because the time-domain and frequency-domain diffusion processes in HyFAD are tightly coupled throughout both the forward and reverse diffusion procedures, simply removing the frequency branch from the denoising network would lead to an incomplete and potentially unfair ablation. Instead, we remove the frequency-domain process from both diffusion stages, which reduces the framework to a purely time-domain diffusion model. This counterpart is equivalent to CSDI, an established generative baseline for time-series imputation. Therefore, the comparison between HyFAD and CSDI serves both as a controlled study of the effect of frequency modeling and as a direct comparison with a representative time-domain generative imputation method.

To further examine whether the imputed series preserves periodic structures across different spectral ranges, we adopt Spectral-MAE and partition the spectrum into three sub-intervals and the results are in Tab.7. The Spectral-MAE is defined as the normalized PSD estimates of the ground-truth series \(P_{\mathrm{GT}}(\omega)\) against those of the reconstruction \(P_{\mathrm{Pred}}(\omega)\): \[\mathrm{S}\text{-}\mathrm{MAE} = \frac{1}{\vert\Omega\vert}\sum_{\omega\in\Omega}\lvert\frac{P_{\mathrm{GT}}(\omega)}{\sum_{\omega^\prime}P_{\mathrm{GT}}(\omega^\prime)}-\frac{P_{\mathrm{pred}}(\omega)}{\sum_{\omega^\prime}P_{\mathrm{pred}}(\omega^\prime)}\rvert\] where \(\Omega\) is the set of evaluated frequencies. The results show that CSDI exhibits systematic weaknesses in capturing both low- and high-frequency patterns, highlighting the importance of explicit frequency-domain modeling. Since Spectral-MAE can be dominated by low-frequency components, the performance improvement is not significant. We additionally compute the error in the logarithmic spectral domain and report band-wise Log-SMAE. Moreover, we report the Leading Frequency Error (LFE), which evaluates whether the dominant frequency component of the original signal is correctly recovered by the imputed sequence. LFE complements Log-SMAE by focusing on the recovery of the dominant periodic structure: a large LFE indicates that the dominant frequency is not well recovered, whereas Log-SMAE measures the fidelity of the full spectral distribution. The log-SMAE and LFE are formulated as: \[\log\mathrm{S}\text{-}\mathrm{MAE} = \frac{1}{\vert\Omega\vert}\sum_{\omega\in\Omega}\lvert\log\frac{P_{\mathrm{GT}}(\omega)}{\sum_{\omega^\prime}P_{\mathrm{GT}}(\omega^\prime)}-\log\frac{P_{\mathrm{pred}}(\omega)}{\sum_{\omega^\prime}P_{\mathrm{pred}}(\omega^\prime)}\rvert\] \[\mathrm{LFE} = \frac{1}{N}\sum_{i=1}^N\vert f_{\mathrm{GT},i}^{\star}-f_{\mathrm{pred,i}}^{\star}\vert,\qquad f_{i}^{\star} = \frac{1}{2\pi}\mathop{\arg\max}_{\omega\in\Omega}P_i(\omega)\]

The additional spectral evaluation results provide further empirical evidence for the effect of explicit frequency modeling. As shown in Table.8, HyFAD achieves lower Log-SMAE than the pure time-domain counterpart across all three frequency bands and all missing ratios. In the low-frequency band, the relative reductions are 24.4%, 19.7%, and 10.8% under 10%, 50%, and 90% missingness, respectively. In the middle-frequency band, the reductions are 22.8%, 17.9%, and 11.2%, while in the high-frequency band, the reductions are 21.3%, 18.1%, and 6.8%.

HyFAD also consistently reduces LFE in the middle- and high-frequency bands across all missing ratios. Specifically, the relative reductions in the middle-frequency band are 25.0%, 16.7%, and 12.2%, and those in the high-frequency band are 20.8%, 19.6%, and 9.8%. In contrast, the low-frequency LFE values of both methods are already very small, indicating that the dominant low-frequency periodicity can be relatively well recovered even by the pure time-domain baseline. This suggests that the main advantage of HyFAD lies in recovering the more challenging middle- and high-frequency periodic structures.

Overall, these results show that the benefit of explicit frequency modeling is not limited to aggregate improvements in standard time-domain metrics. Instead, HyFAD brings consistent improvements in spectral reconstruction, particularly in the middle- and high-frequency regions where the pure time-domain baseline exhibits larger errors. This further validates the effectiveness of incorporating frequency-domain diffusion for modeling heterogeneous spectral components in time-series imputation.

7.13.3 Comparison with LSCD and CSDI↩︎

To better demonstrate the superiority of HyFAD over the two next-best baselines, CSDI and LSCD, We present the results of HyFAD, LSCD and CSDI on 5 independent runs and report the corresponding 95% CI and \(p\)-values in Table.9.

The detailed comparison in Table 9 further confirms the superiority and stability of HyFAD over the two strongest baselines, LSCD and CSDI. HyFAD achieves the lowest mean RMSE and MAE across all datasets and missingness rates. On PhysioNet, HyFAD consistently outperforms LSCD and CSDI under 10%, 50%, and 90% missingness, with especially clear gains in RMSE. For example, compared with LSCD, HyFAD reduces RMSE by 7.4%, 5.5%, and 7.5% under the three missing ratios, respectively. Similar improvements are observed over CSDI, showing that HyFAD remains robust as the missing ratio increases.

The confidence intervals also indicate that HyFAD maintains stable performance across repeated runs, and its intervals are generally lower than those of LSCD and CSDI. Moreover, all reported p-values are below 0.05, demonstrating that the improvements of HyFAD over both LSCD and CSDI are statistically significant. On the Air Quality dataset, HyFAD also achieves the best RMSE and MAE, although the MAE improvement is relatively smaller. Overall, these results show that HyFAD not only improves average imputation accuracy but also provides statistically reliable gains over competitive diffusion-based baselines.

7.13.4 Visualization Results↩︎

The imputation results on Air Quality dataset and PhysioNet dataset with missing ratio 10%, 50% and 90% are presented in Fig.6,7,8,9.

7.14 Social impacts↩︎

This work may have positive societal impacts by improving time series imputation in domains such as healthcare monitoring, air-quality assessment, industrial sensing, and scientific data analysis, where missing observations are common and better reconstruction can support more reliable downstream analysis. However, the method may also pose risks if imputed values are treated as ground-truth measurements, especially in high-stakes settings such as clinical decision-making, public safety, or critical infrastructure. Incorrect or biased imputations could lead to misleading conclusions or harmful decisions, particularly for underrepresented populations, regions, or sensor conditions. Therefore, imputed values should be clearly distinguished from observed data, and deployment in sensitive domains should include domain-specific validation, uncertainty assessment, privacy protection, and expert review.

References↩︎

---

1. Corresponding Author↩︎