有限新息率采样中Fast Cadzow去噪方法研究

Research on the Fast Cadzow Denoising Method in Finite Rate of Innovation Sampling

摘要: 本文研究了有限新息率(Finite Rate of Innovation, FRI)信号采样中的Fast Cadzow迭代去噪方法。针对传统Cadzow算法在高维数据降噪中计算复杂度高、收敛速度慢的问题，使用了一种基于子空间投影优化的Fast Cadzow算法，通过减少奇异值分解(SVD)的计算量，提升了算法的效率。本文详细分析了Fast Cadzow算法的数学原理和实现过程，并通过数值仿真实验对其在不同信噪比(SNR)条件下的降噪效果和计算效率进行了验证。实验结果表明，与传统Cadzow算法相比，Fast Cadzow算法在PSNR、计算效率上均表现出显著优势，尤其在低信噪比环境中展现出更强的抗噪能力。本文的研究成果为FRI信号重构时的高效降噪提供了一种新的解决方案，对实际信号处理中的稀疏信号重构具有重要的应用价值。

Abstract: This paper investigates the Fast Cadzow iterative denoising method in Finite Rate of Innovation (FRI) signal sampling. Addressing the limitations of the traditional Cadzow algorithm, which exhibits high computational complexity and slow convergence in high-dimensional data denoising, we employ an optimized Fast Cadzow algorithm based on subspace projection. This approach significantly reduces the computational load of singular value decomposition (SVD) operations, improving the algorithm’s efficiency. We thoroughly analyze the mathematical principles and implementation of the Fast Cadzow algorithm and validate its denoising performance and computational efficiency through numerical simulations under different signal-to-noise ratio (SNR) conditions. Experimental results demonstrate that, compared to the traditional Cadzow algorithm, the Fast Cadzow algorithm achieves significant improvements in PSNR and computational efficiency, with stronger noise resistance, especially in low SNR environments. The findings of this study offer an effective solution for efficient denoising in FRI signal reconstruction and hold substantial application value for sparse signal reconstruction in practical signal processing scenarios.

1. 引言

有限新息率(Finite Rate of Innovation, FRI)信号处理是现代信号处理领域中的重要技术，广泛应用于需要从稀疏采样数据中恢复信号的诸多场景，例如雷达成像、医学成像、地震勘探等。在这些应用中，信号往往具有稀疏特性，其新息率(即信号中的非零信息)远低于经典采样定理所要求的采样频率。这类信号采样和重构问题的独特性在于，信号的恢复不依赖于其频带宽度，而是基于有限的新息率，即信号每周期由有限的参数。Vetterli等人[1]首次提出了针对这类信号的采样和重构方法，在有限新息率信号的降噪处理中，噪声的干扰往往会显著影响信号的重建精度。为此，研究者提出了多种降噪算法，其中Cadzow算法凭借其有效的低秩逼近特性成为一种经典方法。Cadzow算法通过将信号转换为Hankel矩阵形式，然后进行奇异值分解(Singular Value Decomposition, SVD) [2]，对矩阵进行低秩逼近，从而去除噪声。然而，传统的Cadzow算法在处理高维数据时计算复杂度较高，尤其在FRI信号的迭代处理中，其SVD的计算成本较大，导致算法收敛速度缓慢，限制了实际应用的效率和范围。因此，提升Cadzow算法在FRI信号降噪中的效率成为当前研究的关键问题。

为了解决传统Cadzow算法的计算瓶颈，近年来提出的Fast Cadzow算法[3]引入了基于子空间投影的优化策略，通过减少每次迭代中的计算量来加速收敛。Fast Cadzow算法在经典Cadzow算法的基础上，采用了额外的矩阵子空间投影，限制了SVD操作的计算维度，从而显著提升了迭代效率。在信号重建质量不降低的前提下，Fast Cadzow算法的计算复杂度较传统Cadzow算法有显著降低，使其在高维信号的快速降噪和重构中具有独特优势[4]。

本文聚焦于有限新息率采样中的Fast Cadzow迭代去噪方法，详细探讨该算法的数学原理和实现过程，并通过数值仿真实验验证其在稀疏信号重建中的有效性。本文将针对有限新息率采样信号的特点，引入Fast Cadzow算法，并详细分析其在降噪处理中的理论基础。再通过数值仿真实验验证Fast Cadzow算法在迭代速度和计算效率上的显著提升，同时分析不同噪声水平下算法的重建效果。对比传统Cadzow算法和Fast Cadzow算法在有限新息率采样中的表现，阐明Fast Cadzow算法在信号处理效率和降噪精度方面的优势。通过本文的研究，希望能够为有限新息率信号处理中的降噪问题提供一种高效、准确的解决方案，为FRI信号的降噪与重建应用提供理论支持和实践参考。

2. 有限新息率

2.1. 新息率与FRI信号

Vetterli等人提出了一种技术，针对一类既不受带宽限制也不属于移位不变空间的信号而进行采样和重构。这类信号被称为有限采样率(FRI)信号，完全由每个周期的有限个参数来确定。Vetterli等人表明，狄拉克脉冲流、非均匀样条和非均匀多项式样条属于FRI信号。给定一组已知的脉冲基函数$h_k\left(t\right)$ ，任意的幅值$a_k$ 与时延$t_k$ ，就能够组成以下的脉冲信号$x\left(t\right)$ ：

$x\left(t\right)={\displaystyle \sum _{k=1}^K{a}_k{h}_k}\left(t-t_k\right),t\in \left[0,T\right)$ (1)

其中K为脉冲的个数，T为脉冲串的时间持续长度。

在(1)中，脉冲基函数$h_k\left(t\right)$ 是确定的，那么只需要将幅值$a_k$ 与时延$t_k$ 估计出来，就能够将原始信号$x\left(t\right)$ 重构出来。

在这里引入计数函数$C_x\left(t_a,t_b\right)$ 来计算$\left[t_a,t_b\right]$ τ时间内脉冲信号$x\left(t\right)$ 的自由度的个数，也就是未知参量的个数，则新息率$\rho$ 的定义为：

$\rho =\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{C}_x\left(-\frac{\tau }{2},\frac{\tau }{2}\right)$ (2)

而新息率$\rho$ 如果是有限的信号就被称作有限新息率信号(FRI信号)如图1所示。在有限新息率理论中，脉冲基函数$h_k\left(t\right)$ 一般为狄拉克脉冲、高斯脉冲和洛伦兹脉冲等[5]-[7]，即：

$h_k\left(t\right)=\delta \left(t\right)$ (3)

$h_k\left(t\right)={\text{e}}^{-\frac{t^2}{2{\sigma }^2}}$ (4)

$h_k\left(t\right)=\frac{1}{1+{\left(\frac{t}{\theta }\right)}^2}$ (5)

其中$\delta \left(t\right)$ 为狄拉克脉冲，$\sigma$ 为高斯脉冲的标准差，决定着高斯脉冲持续时间以及频谱宽度，$\theta$ 控制洛伦兹脉冲的宽度[8]。每种FRI信号都是一定的新息率，Dirac脉冲和高斯脉冲均有幅值$a_k$ 和时延$t_k$ 所以其新息率为$\rho =2K/T$ ，而洛伦兹脉冲除此之外还有一个新的自由度$\theta$ 所以新息率$\rho =3K/T$ 。

Figure 1. Common FRI signals. (a) Dirac pulse, (b) Gaussian pulse, (c) Lorentzian pulse

图1. 常见的FRI信号。(a) 狄拉克脉冲，(b) 高斯脉冲，(c) 洛伦兹脉冲

2.2. FRI采样结构与信号重构

经典FRI采样框架是一个单独的采样通道，如图2：

Figure 2. FRI sampling framework

图2. FRI采样框架

基本原理就是原始信号$x\left(t\right)$ 通过采样核$g\left(t\right)$ 进行滤波，采样核的本质就是一个低通滤波器，然后再对低通滤波后的待采信号$y\left(t\right)$ 按照一定的采样率$T_s$ 进行均匀采样[9]，得到采样结果$y\left[n\right]$ ($n=0,1,\cdots ,L-1$ )。则采样结果可表示为：

$y\left[n\right]={\displaystyle {\int }_{-\infty }^{\infty }x\left(t\right)g\left(t-n{T}_s\right)\text{d}t}=\langle g\left(t-n{T}_s\right),x\left(t\right)\rangle n=0,1,\cdots ,L-1$ (6)

则经典FRI采样步骤如下。

通过采样值来计算采样值的傅里叶级数系数。采样结果$y\left[n\right]$ 如下式子：

$y\left[n\right]={\displaystyle {\int }_{-\infty }^{\infty }x\left(t\right)g\left(t-n{T}_s\right)\text{d}t}n=0,1,\cdots ,L-1$ (7)

其中L表示采样点数。将$y\left[n\right]$ 中的$x\left(t\right)$ 用如下傅里叶级数形式来替换，

$x\left(t\right)={\displaystyle \sum _{m\in \mathbb{Z}}X\left[m\right]{\text{e}}^{j\frac{2\pi }{T}mt}}t\in \left[0,T\right)$ (8)

就能得到其傅里叶系数$X\left[m\right]$ ：

$\begin{array}{c}X\left[m\right]=\frac{1}{T}{\displaystyle {\int }_0^{T}x\left(t\right){\text{e}}^{-j\frac{2\pi }{T}mt}\text{d}t}\\ =\frac{1}{T}{\displaystyle \sum _{k=1}^K{a}_k}{\displaystyle {\int }_0^T{h}_k\left(t-t_k\right){\text{e}}^{-j\frac{2\pi }{T}mt}\text{d}t}\\ =\frac{1}{T}H\left(\frac{2\pi m}{T}\right){\displaystyle \sum _{k=1}^K{a}_k}{\text{e}}^{-j\frac{2\pi }{T}m{t}_k}\end{array}$ (9)

其中的$H\left(\omega \right)$ 就是公式(2)中的脉冲基函数$h_k\left(t\right)$ 的连续傅里叶变换，假设直接令$H\left[m\right]=\frac{1}{T}H\left(\frac{2\pi m}{T}\right)$ ，则采样样本$Y\left[m\right]$ 如下：

$Y\left[m\right]=X\left[m\right]/H\left[m\right]={\displaystyle \sum _{k=1}^K{a}_k}{\text{e}}^{-j\frac{2\pi }{T}m{t}_k}$ (10)

再做一次代换，联立

${\omega }_k=\frac{2\pi }{T}{t}_k,u_k={\text{e}}^{-j{\omega }_k}$ (11)

则公式(10)可变换为：

$Y\left[m\right]={\displaystyle \sum _{k=1}^K{a}_k}{u}_k^m$ (12)

通过将$Y\left[m\right]$ 中的$a_k,u_k$ 求解出来，即是从频率谱中估计出原始信号的角频率${\omega }_k$ 以及幅度$a_k$ ，就将重构问题转化为了谱估计问题。

零化滤波器法[10]便是解决这类谱估计问题的方法，核心思想是构造一个特定的滤波器，该滤波器对信号中的每个非零频率成分(如谐波)进行零化，从而更好的对样本进行估计求解。经典FRI采样理论的重构算法也是使用该方法，步骤如下。

首先构造一个多项式滤波器$H\left(z\right)$ ，使得其零点位置对应于信号中的每个频率成分。

$H\left(z\right)={\displaystyle \prod _{k=1}^K\left(1-u_k{z}^{-1}\right)}={\displaystyle \sum _{k=0}^K{h}_k{z}^{-k}}$ (13)

对式13进行Z变换：

${\displaystyle \sum _{i=0}^K{h}_i}Y\left[m-i\right]={\displaystyle \sum _{i=0}^K{h}_i}{\displaystyle \sum _{k=0}^{K-1}{a}_k}{u}_k^{m-i}={\displaystyle \sum _{k=0}^{K-1}{a}_k}{u}_k^m{\displaystyle \sum _{i=0}^K{h}_i}{u}_k^{-i}=0$ (14)

其中${\mu }_k$ 为零化滤波器的根，$h_k$ 为滤波器的系数，零化的含义就是该滤波器与$Y\left[m\right]$ 的内积为0，令$h_0=1$ ，则上式矩阵表达则为：

$\left(\begin{array}{cccc}Y\left[-1\right]& Y\left[-2\right]& \cdots & Y\left[-K\right]\\ Y\left[0\right]& Y\left[-1\right]& \dots & Y\left[-K+1\right]\\ ⋮& ⋮& \ddots & ⋮\\ Y\left[K-1\right]& Y\left[K-2\right]& \cdots & Y\left[-1\right]\end{array}\right)\left(\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_K\end{array}\right)=-\left(\begin{array}{c}Y\left[0\right]\\ Y\left[1\right]\\ ⋮\\ Y\left[K\right]\end{array}\right)$ (15)

根据矩阵推导出只需要$2K+1$ 个傅里叶级数系数就可以计算零化滤波器的系数$h_k$ ，再将$h_k$ 代入到式(12)，便能求出零化滤波器的根${\mu }_k$ ，再将得到根直接代入式(11)，则可恢复出时延$t_k$ ：

$t_k=-\frac{\tau \cdot \angle u_k}{2\pi }$ (16)

根据${\mu }_k$ 构造一个K列的范德蒙德矩阵：

$\left(\begin{array}{c}Y\left[0\right]\\ Y\left[1\right]\\ ⋮\\ Y\left[K-1\right]\end{array}\right)=\frac{1}{\tau }\left(\begin{array}{cccc}1& 1& \cdots & 1\\ u_0& u_1& \cdots & u_{K-1}\\ ⋮& ⋮& \ddots & ⋮\\ u_0^{K-1}& u_1^{K-1}& \cdots & u_{K-1}{}^{K-1}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{K-1}\end{array}\right)$ (17)

利用最小二乘法即可估计得到幅值$a_k$ 。至此便重构得到了时延$t_k$ 与幅值$a_k$ 。

零化滤波器法虽然从理论上给出了FRI信号的临界采样率[11]，但在高新息率信号的处理上存在局限性。随着信号新息率的增加，零化滤波器法在多项式求根时的计算量急剧增加[12]，导致重构效率显著下降。此外，在噪声环境下，该方法的重构效果较差，因为噪声通常具有无限的新息率，严重影响了重构精度。因此，需要探索和采用其他改进算法，以提高其在复杂环境中的重构性能。

3. Fast Cadzow算法

Cadzow算法是一种经典的低秩矩阵逼近方法，通过构造Hankel矩阵并进行SVD，可以有效去除噪声并恢复信号的稀疏结构。为了解决处理FRI信号时，传统Cadzow算法存在计算复杂度高和收敛速度慢的问题，为了解这个问题，本文使用了一种优化的Fast Cadzow算法[13]，通过引入快速矩阵投影，显著提升了算法的效率和收敛速度。使用了Fast Cadzow降噪算法的FRI采样框架如图3所示：

Figure 3. FRI sampling framework with Fast Cadzow noise reduction algorithm

图3. 包含Fast Cadzow降噪算法的FRI采样框架

3.1. Fast Cadzow算法的核心思想

Fast Cadzow算法的核心在于通过引入一个矩阵子空间投影(subspace projection)步骤，减少了每次SVD分解的维数，从而降低计算复杂度。具体而言，在每次迭代中，Fast Cadzow算法首先将Hankel矩阵投影到一个低维子空间中，再进行SVD分解和低秩逼近。该子空间投影操作使得矩阵在SVD时的计算量显著减少，从而提高了算法的整体效率。

Fast Cadzow算法的迭代过程包括以下步骤：

1) Hankel化：将初始信号向量$z_0=y$ 转换为Hankel矩阵$H_{Z_k}$ ，构造如下：

$H{z}_k=\left[\begin{array}{cccc}{z}_k\left[0\right]& z_k\left[1\right]& \cdots & z_k\left[K-1\right]\\ z_k\left[1\right]& z_k\left[2\right]& \cdots & z_k\left[K\right]\\ ⋮& ⋮& \ddots & ⋮\\ z_k\left[L-1\right]& z_k\left[L\right]& \cdots & z_k\left[N-1\right]\end{array}\right]$ (18)

其中，N是信号长度，L和K分别为Hankel矩阵的行数和列数。

2) 子空间投影：根据当前迭代结果，构造投影子空间$T_k$ 并将$H_{Z_k}$ 投影到该子空间中。在第一次迭代时，设投影子空间$T_0$ 为整个空间。在随后的迭代中，基于上一次迭代的低秩近似矩阵$L_k$ 来确定子空间$T_k$ 。假设$L_k$ 的SVD分解为$L_k=U_k\Sigma V_k^{*}$ ，则$T_k$ 是由$U_k$ 和$V_k$ 的列空间和行空间所生成的直和空间，投影子空间$T_k$ 可以表示为：

$T_k=\left\{U_k{B}^{*}+C{V}_k^{*}|B\in {ℂ}^{K\times r},C\in {ℂ}^{L\times r}\right\}$ (19)

在该子空间上进行投影，得到：

$W_k=P_{T_k}\left(H_{Z_k}\right)$ (20)

3) 低秩逼近：对投影后的矩阵$W_k$ 进行SVD分解，保留前r个奇异值得到低秩近似矩阵$L_{K+1}$ ：

$L_{K+1}=T_r\left(W_k\right)$ (21)

其中$T_r$ 表示取矩阵的前r个奇异值的操作，从而实现矩阵的低秩逼近。

4) 逆Hankel化：将低秩矩阵$L_{k+1}$ 通过逆Hankel化操作转换回一维信号$z_{k+1}=H^{†}{L}_{k+1}$ ，用于下一次迭代。

5) 收敛判断：计算当前迭代结果$z_{k+1}$ 与上一轮结果$z_k$ 的差值。如果满足收敛条件，即$\Vert z_{k+1}-z_k\Vert <\text{tol}$ (其中tol为预设的阈值)，则算法终止；否则继续迭代。

通过上述步骤，Fast Cadzow算法将SVD操作限制在低维子空间中，极大地减少了每次迭代的计算复杂度，使得算法的收敛速度和计算效率大幅提升。

3.2. Fast Cadzow算法的优势

计算效率的提升：传统Cadzow算法在每次迭代中都需要对Hankel矩阵进行SVD，这一步骤的计算复杂度为$O\left(N^3\right)$ (其中N为信号的长度)，在处理大规模信号时计算代价显著。而Fast Cadzow算法通过在每次迭代中引入子空间投影操作，将SVD操作限制在一个低维子空间内，从而有效地降低了每次迭代的计算复杂度。具体来说，Fast Cadzow算法的复杂度降低为$O\left(N{r}^2\right)$ (其中r为近似的低秩)，显著减少了计算成本。因此，Fast Cadzow算法特别适合处理大规模或高维度的信号。

收敛速度提高：在每次迭代中，Fast Cadzow算法通过将Hankel矩阵投影到子空间内，使得每次迭代后的低秩逼近结果更加接近原始信号的特征，从而加快了算法的收敛速度。相较于传统Cadzow算法在每轮迭代中重复全局SVD计算的方式，Fast Cadzow的收敛速度更快，能够在较少的迭代次数内逼近目标信号。这一改进使得算法在同等去噪效果下能够节省时间成本，特别是在噪声较大的环境中，算法的快速收敛特性尤为明显。

4. 实验结果

4.1. 实验设置

本实验是为了验证Fast Cadzow算法在FRI信号降噪中的性能。通过将Fast Cadzow与传统Cadzow算法进行对比，评价其在降噪效果、收敛速度、计算效率等方面的表现。

实验信号：生成长度为N = 100的稀疏信号，包含5个随机位置的狄拉克脉冲。生成信号后，加入不同信噪比的高斯白噪声。

噪声模型：信噪比设置为10 dB、20 dB和30 dB，以评估算法在不同噪声环境下的表现。

实验参数：迭代次数设为10次，以对比两种算法在不同SNR下的收敛速度。

实验环境：MATLAB。

4.2. 评价指标

1) 峰值信噪比(PSNR)：PSNR用于评价重建信号的质量，定义为：

$\text{PSNR}=10\cdot {\log }_{10}\left(\frac{\max {\left(\text{original-signal}\right)}^2}{\text{MSE}}\right)$ (22)

其中，MSE (均方误差)为重建信号和原始信号之间的误差：

$\mathrm{MSE}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(\text{ denoised\_signal}\left[i\right]-\text{original\_signal}\left[i\right]\right)}^2}$ (23)

PSNR越高表示重建信号越接近原始信号。

2) 计算效率：通过记录每种算法的平均计算时间，比较两者的计算效率。

4.3. 实验结果及分析

在不同SNR条件下，Fast Cadzow算法在PSNR方面均表现出优于传统Cadzow算法的效果。图4是不同SNR下的PSNR曲线对比图。

SNR = 10 dB：两种算法均能够恢复信号，但Fast Cadzow算法在前几次迭代中收敛更快，PSNR显著提高。

SNR = 20 dB：两种算法的PSNR差距缩小，但Fast Cadzow依旧在迭代次数和PSNR上占优。

SNR = 30 dB：在噪声较低的条件下，Fast Cadzow算法与传统Cadzow的PSNR相近，但前者的收敛速度更快。

Fast Cadzow算法在不同噪声水平下的计算时间均低于传统Cadzow算法。这是由于Fast Cadzow通过子空间投影降低了每次迭代的计算量。表1为不同SNR下的平均计算时间：

Table 1. The average calculation time under different SNRS

表1. 不同SNR下的平均计算时间

SNR (dB)	传统Cadzow算法(秒)	Fast Cadzow算法(秒)
10	5.34	2.14
20	4.76	1.95
30	3.98	1.76

Figure 4. PSNR curve comparison diagram

图4. PSNR曲线对比图

实验结果表明，Fast Cadzow算法在保持较高降噪效果的同时，提高了计算效率和收敛速度。特别是在低噪声条件下，Fast Cadzow算法的优势更加明显，这说明其通过子空间投影所实现的加速效果，使其适用于高效FRI信号的降噪任务。传统Cadzow算法在收敛速度和计算效率上均略显不足，而Fast Cadzow算法则提供了较优的降噪与重构性能。

5. 结论

本文以有限新息率(FRI)信号为研究对象，使用了一种基于子空间投影优化的Fast Cadzow迭代去噪方法，并通过数值实验分析了算法的降噪性能和计算效率。主要结论如下：

1. 构建了Fast Cadzow算法的迭代模型，通过引入子空间投影，显著降低了奇异值分解(SVD)的计算复杂度。在不同信噪比(SNR)条件下，实验结果表明该方法能够快速收敛，有效减少每次迭代的计算时间，验证了算法的高效性。

2. 实验表明，Fast Cadzow算法在低信噪比条件下展现了较强的抗噪性能。在信噪比为10 dB时，峰值信噪比(PSNR)较传统Cadzow算法提升约2 dB；在高信噪比(20 dB、30 dB)条件下，仍能保持良好的降噪效果，说明算法在稀疏信号重构中的稳定性较强。

以上研究结果为有限新息率信号处理中的高效去噪提供了一种可行的解决方案，并为复杂信号的重构和应用提供了理论依据和实践参考。