自相关滤波详解

自相关滤波详解

自相关滤波是信号处理领域中一种重要的技术，它利用信号与其自身的相关性来提取有用信息，滤除噪声，增强特定信号特征。

自相关的基本概念

自相关是描述信号在不同时间点之间相似度的一个重要指标。对于一个确定性连续时间信号 $x(t)$，其自相关函数 $R_{xx}(\tau )$ 定义为：

$$R_{xx}(\tau )=\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_{-T}^{T}x(t)x(t+\tau )dt$$

对于随机过程，自相关函数定义为该信号与其时间延迟版本的乘积期望值：

$$R_{xx}(\tau )=\mathbb{E}[x(t)x(t+\tau )]$$

其中 $\mathbb{E}[\cdot ]$ 表示期望操作，$\tau$ 是时间延迟。对于离散时间信号 $x[n]$，自相关函数可以表示为：

$$R_{xx}[k]=\mathbb{E}[x[n]x[n+k]]$$

若信号是平稳随机过程，则其自相关函数只与时间延迟 $k$ 有关，而与具体时间 $n$ 无关。在实际应用中，由于我们只能获取有限长度的信号，自相关函数通常通过时间平均的方式来估计。对于长度为 $N$ 的离散信号 $x[n],n=0,1,...,N-1$，有偏自相关估计为：

$${\hat{R}}_{xx}[k]=\frac{1}{N}\sum _{n=0}^{N-|k|-1}x[n]x[n+|k|],|k|<N$$

无偏自相关估计为：

$${\hat{R}}_{xx}[k]=\frac{1}{N-|k|}\sum _{n=0}^{N-|k|-1}x[n]x[n+|k|],|k|<N$$

自相关函数具有几个重要的数学性质：

1. 对称性：$R_{xx}(-\tau )=R_{xx}(\tau )$ 对于实值信号
2. 最大值原理：$|R_{xx}(\tau )|\le R_{xx}(0)$ 对于任意 $\tau$
3. 自相关函数的傅里叶变换等于信号功率谱密度： S x x ( f ) = F { R x x ( τ ) } S_{xx}(f) = \mathcal{F}\{R_{xx}(\tau)\} Sxx​(f)=F{Rxx​(τ)}

最后一条性质是维纳-辛钦定理（Wiener-Khinchin Theorem），是频谱分析的理论基础之一。对于能量有限的确定性信号，功率谱密度可以表示为：

$$S_{xx}(f)=|X(f){|}^2$$

其中 $X(f)$ 是信号 $x(t)$ 的傅里叶变换。

自相关滤波的基本原理

自相关滤波的核心思想是利用信号自相关特性来区分有用信号和噪声。在许多实际应用中，有用信号往往具有强的相关性，而噪声（特别是白噪声）则具有弱的相关性。考虑一个包含噪声的连续时间信号模型：

$$y(t)=s(t)+w(t)$$

其中 $s(t)$ 是有用信号，$w(t)$ 是加性噪声。假设噪声是零均值的，并且与信号不相关，则该混合信号的自相关函数为：

$$R_{yy}(\tau )=R_{ss}(\tau )+R_{ww}(\tau )+R_{sw}(\tau )+R_{ws}(\tau )$$

如果噪声与信号不相关，即 $R_{sw}(\tau )=R_{ws}(\tau )=0$，则：

$$R_{yy}(\tau )=R_{ss}(\tau )+R_{ww}(\tau )$$

如果噪声是白噪声，其自相关函数为：

$$R_{ww}(\tau )={\sigma }_w^2\delta (\tau )$$

其中 ${\sigma }_w^2$ 是噪声方差，$\delta (\tau )$ 是狄拉克德尔塔函数。对于离散情况：

$$R_{ww}[k]={\sigma }_w^2\delta [k]$$

因此，当 $\tau \ne 0$ 或 $k\ne 0$ 时，$R_{yy}(\tau )=R_{ss}(\tau )$ 或 $R_{yy}[k]=R_{ss}[k]$，这意味着非零延迟的自相关值主要由信号决定，噪声的贡献被大大减小。自相关滤波就是基于这一原理，通过对信号进行自相关处理，然后应用适当的滤波器来提取所需的信息。

自相关滤波器的数学模型

一个一般化的自相关滤波器可以用如下操作序列描述：

1. 计算输入信号 $y(t)$ 的自相关函数 $R_{yy}(\tau )$
2. 设计和应用一个滤波器 $h(\tau )$ 到自相关函数
3. 从滤波后的自相关函数中重构出增强的信号

在数学上，滤波过程可以表示为自相关函数与滤波器核的卷积：

$${\hat{R}}_{ss}(\tau )=(R_{yy}\ast h)(\tau )={\int }_{-\infty }^{\infty }{R}_{yy}(\alpha )h(\tau -\alpha )d\alpha$$

对于离散信号：

$${\hat{R}}_{ss}[k]=(R_{yy}\ast h)[k]=\sum _{m=-\infty }^{\infty }{R}_{yy}[m]h[k-m]$$

在频域中，这相当于：

$${\hat{S}}_{ss}(f)=S_{yy}(f)\cdot H(f)$$

其中 ${\hat{S}}_{ss}(f)$, $S_{yy}(f)$ 和 $H(f)$ 分别是 ${\hat{R}}_{ss}(\tau )$, $R_{yy}(\tau )$ 和 $h(\tau )$ 的傅里叶变换。滤波器 $h(\tau )$ 或 $H(f)$ 的设计取决于信号和噪声的特性以及具体应用目标。一种理论上最优的滤波器是维纳滤波器，它基于最小均方误差准则：

$$H_{wiener}(f)=\frac{S_{ss}(f)}{S_{ss}(f)+S_{ww}(f)}$$

其中 $S_{ss}(f)$ 和 $S_{ww}(f)$ 分别是信号和噪声的功率谱密度。从滤波后的自相关函数重构信号是一个挑战性问题。对于随机信号，完全精确的重构通常是不可能的，因为自相关处理丢失了相位信息。然而，对于某些特定类型的信号，如最小相位信号，可以通过求解同态方程来实现重构：

$$\log {\hat{S}}_{ss}(f)=\log |\hat{X}(f){|}^2$$

然后通过复倒谱分析计算：

x ^ ( t ) = F − 1 { e 1 2 log ⁡ S ^ s s ( f ) + j ϕ m i n ( f ) } \hat{x}(t) = \mathcal{F}^{-1}\{e^{\frac{1}{2}\log \hat{S}_{ss}(f) + j\phi_{min}(f)}\} x^(t)=F−1{e21​logS^ss​(f)+jϕmin​(f)}

其中 ${\varphi }_{min}(f)$ 是与 $\log {\hat{S}}_{ss}(f)$ 通过希尔伯特变换关联的最小相位函数：

ϕ m i n ( f ) = − H { log ⁡ S ^ s s ( f ) } \phi_{min}(f) = -\mathcal{H}\{\log \hat{S}_{ss}(f)\} ϕmin​(f)=−H{logS^ss​(f)}

其中 H { ⋅ } \mathcal{H}\{\cdot\} H{⋅} 表示希尔伯特变换：

H { g ( f ) } = 1 π ∫ − ∞ ∞ g ( ξ ) f − ξ d ξ \mathcal{H}\{g(f)\} = \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{g(\xi)}{f - \xi} d\xi H{g(f)}=π1​∫−∞∞​f−ξg(ξ)​dξ

这种方法在语音处理和地震信号分析中有应用。

自相关滤波器的频域解释

自相关滤波在频域中有清晰的解释。根据维纳-辛钦定理，信号的功率谱密度（PSD）是其自相关函数的傅里叶变换：

S y y ( f ) = F { R y y ( τ ) } = ∫ − ∞ ∞ R y y ( τ ) e − j 2 π f τ d τ S_{yy}(f) = \mathcal{F}\{R_{yy}(\tau)\} = \int_{-\infty}^{\infty} R_{yy}(\tau) e^{-j2\pi f\tau} d\tau Syy​(f)=F{Ryy​(τ)}=∫−∞∞​Ryy​(τ)e−j2πfτdτ

对于离散信号：

$$S_{yy}(f)=\sum _{k=-\infty }^{\infty }{R}_{yy}[k]e^{-j2\pi fk}$$

在频域中，自相关滤波可以表示为：

$${\hat{S}}_{ss}(f)=S_{yy}(f)\cdot H(f)$$

对于白噪声污染的信号，其功率谱密度为：

$$S_{yy}(f)=S_{ss}(f)+{\sigma }_w^2$$

通过设计适当的频域滤波器 $H(f)$，我们可以增强特定频率范围内的信号分量，同时抑制噪声。对于维纳滤波器：

$$H_{wiener}(f)=\frac{S_{ss}(f)}{S_{ss}(f)+S_{ww}(f)}$$

当信噪比（SNR）很高 $S_{ss}(f)\gg S_{ww}(f)$ 时，$H_{wiener}(f)\approx 1$；当SNR很低 $S_{ss}(f)\ll S_{ww}(f)$ 时，$H_{wiener}(f)\approx 0$。维纳滤波器的输出信号功率谱密度为：

$${\hat{S}}_{ss}(f)=\frac{S_{ss}^2(f)}{S_{ss}(f)+S_{ww}(f)}$$

维纳滤波器的均方误差可以表示为：

$$MSE={\int }_{-\infty }^{\infty }\frac{S_{ss}(f)S_{ww}(f)}{S_{ss}(f)+S_{ww}(f)}df$$

对于有色噪声情况，预白化滤波可以提高性能。预白化滤波器的频率响应为：

$$W(f)=\frac{1}{\sqrt{S_{ww}(f)}}$$

预白化后的信号谱为：

$$S_{y^{\prime }}(f)=W^2(f)S_{yy}(f)=\frac{S_{ss}(f)}{S_{ww}(f)}+1$$

在预白化域中应用维纳滤波：

$$H_{wiener}^{\prime }(f)=\frac{\frac{S_{ss}(f)}{S_{ww}(f)}}{\frac{S_{ss}(f)}{S_{ww}(f)}+1}=\frac{S_{ss}(f)}{S_{ss}(f)+S_{ww}(f)}$$

总体滤波器为：

$$H_{total}(f)=W(f)\cdot H_{wiener}^{\prime }(f)=\frac{S_{ss}(f)}{S_{ww}(f)(S_{ss}(f)+S_{ww}(f))}$$

自相关滤波的时域实现

在时域中实现自相关滤波通常涉及以下数学步骤：

1. 估计输入信号的自相关函数：
   对于长度为N的离散信号 $y[n]$，其自相关估计为：

   $${\hat{R}}_{yy}[k]=\frac{1}{N-|k|}\sum _{n=0}^{N-|k|-1}y[n]y[n+|k|],|k|<N$$

2. 设计滤波器时域响应：
   对于给定的频域目标响应 $H(f)$，可以通过逆傅里叶变换计算时域响应：

   h [ n ] = F − 1 { H ( f ) } h[n] = \mathcal{F}^{-1}\{H(f)\} h[n]=F−1{H(f)}

   实际中，通常会应用窗口函数 $w[n]$ 以减小频谱泄漏：

   $$h_w[n]=h[n]\cdot w[n]$$

3. 对自相关函数应用滤波器：

   $${\hat{R}}_{ss}[k]=\sum _{m=-M}^M{\hat{R}}_{yy}[k-m]h_w[m]$$

   其中 $M$ 是滤波器长度的一半。

4. 从滤波后的自相关函数重构信号：
   对于自回归（AR）过程，可以通过求解Yule-Walker方程组获得AR参数：

   $$\sum _{i=1}^p{a}_i{\hat{R}}_{ss}[k-i]=-{\hat{R}}_{ss}[k],k=1,2,...,p$$

   或矩阵形式：

   $$\left[\begin{array}{cccc}{\hat{R}}_{ss}[0]& {\hat{R}}_{ss}[1]& \cdots & {\hat{R}}_{ss}[p-1]\\ {\hat{R}}_{ss}[1]& {\hat{R}}_{ss}[0]& \cdots & {\hat{R}}_{ss}[p-2]\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{R}}_{ss}[p-1]& {\hat{R}}_{ss}[p-2]& \cdots & {\hat{R}}_{ss}[0]\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}-{\hat{R}}_{ss}[1]\\ -{\hat{R}}_{ss}[2]\\ ⋮\\ -{\hat{R}}_{ss}[p]\end{array}\right]$$

   求解得到AR参数 { a i } \{a_i\} {ai​} 后，可以计算激励方差：

   $${\sigma }_e^2={\hat{R}}_{ss}[0]+\sum _{i=1}^p{a}_i{\hat{R}}_{ss}[i]$$

   然后生成白噪声序列 $e[n]$，方差为 ${\sigma }_e^2$，通过AR模型方程重构信号：

   $$\hat{s}[n]=e[n]-\sum _{i=1}^p{a}_i\hat{s}[n-i]$$

对于更复杂的ARMA过程，可以使用迭代方法，如迭代预白化或最大似然估计来识别模型参数。

自相关滤波的统计特性分析

自相关滤波器的性能分析需要考虑统计特性。对于输入信号 $y[n]=s[n]+w[n]$，自相关估计的均值为：

$$\mathbb{E}[{\hat{R}}_{yy}[k]]=R_{ss}[k]+R_{ww}[k]$$

对于白噪声 $w[n]$，$R_{ww}[k]={\sigma }_w^2\delta [k]$，因此：

$$\mathbb{E}[{\hat{R}}_{yy}[k]]=R_{ss}[k]+{\sigma }_w^2\delta [k]$$

自相关估计的方差与信号的四阶矩（四次累积量）有关。对于高斯过程，四阶矩可以用二阶矩表示：

$$\text{Var}[{\hat{R}}_{yy}[k]]\approx \frac{1}{N}\sum _{m=-\infty }^{\infty }(R_{yy}[m+k]R_{yy}[m-k]+R_{yy}[m]R_{yy}[m+2k])$$

对于长度为 $N$ 的信号，自相关估计的方差近似为：

$$\text{Var}[{\hat{R}}_{yy}[k]]\approx \frac{1}{N}\sum _{m=-(N-1)}^{N-1}(1-\frac{|m|}{N})(R_{yy}[m+k]R_{yy}[m-k])$$

功率谱密度估计的方差分析更为复杂。对于有偏谱估计，在频率 $f$ 处的方差为：

$$\text{Var}[{\hat{S}}_{yy}(f)]\approx \frac{1}{N}\sum _{k=-(N-1)}^{N-1}(1-\frac{|k|}{N})w[k{]}^2\left(S_{yy}^2(f)+\frac{1}{N}\sum _{m=-(N-1)}^{N-1}{S}_{yy}(f-\frac{m}{N})S_{yy}(f+\frac{m}{N})\right)$$

其中 $w[k]$ 是窗口函数。对于维纳滤波器，输出信噪比可以表示为：

$$SN{R}_{out}(f)=\frac{|H_{wiener}(f){|}^2{S}_{ss}(f)}{|1-H_{wiener}(f){|}^2{S}_{ss}(f)+|H_{wiener}(f){|}^2{S}_{ww}(f)}$$

代入维纳滤波器表达式，可得：

$$SN{R}_{out}(f)=\frac{S_{ss}^2(f)}{S_{ww}(f)\cdot S_{ss}(f)}=\frac{S_{ss}(f)}{S_{ww}(f)}=SN{R}_{in}(f)$$

这表明维纳滤波器不会提高单个频率分量的信噪比，而是通过优化各频率分量的权重来最小化总体均方误差。滤波后信号的自相关函数与原始信号自相关函数的关系为：

R ^ s s [ k ] = F − 1 { S ^ s s ( f ) } = F − 1 { S y y ( f ) ⋅ H ( f ) } \hat{R}_{ss}[k] = \mathcal{F}^{-1}\{\hat{S}_{ss}(f)\} = \mathcal{F}^{-1}\{S_{yy}(f) \cdot H(f)\} R^ss​[k]=F−1{S^ss​(f)}=F−1{Syy​(f)⋅H(f)}

$${\hat{R}}_{ss}[k]=(R_{yy}\ast r_h)[k]$$

其中 $r_h[k]$ 是滤波器 $h[n]$ 的自相关函数。

自相关滤波的数学框架

从更一般的角度看，自相关滤波可以在函数空间中表述为一个投影问题。令 $\mathcal{H}$ 表示希尔伯特空间，包含所有可能的信号，$\mathcal{S}\subset \mathcal{H}$ 是信号子空间，$\mathcal{W}\subset \mathcal{H}$ 是噪声子空间。理想的自相关滤波是将观测信号 $y$ 投影到信号子空间 $\mathcal{S}$ 上：

$$\hat{s}=P_{\mathcal{S}}(y)$$

其中 $P_{\mathcal{S}}$ 是投影算子。

在子空间方法中，信号的自相关矩阵可以进行特征分解：

$${\mathbf{R}}_{yy}=\sum _{i=1}^N{\lambda }_i{\mathbf{u}}_i{\mathbf{u}}_i^H$$

其中 ${\lambda }_i$ 是特征值，${\mathbf{u}}_i$ 是对应的特征向量。

当信号空间维度为 $K<N$，且信噪比足够高时，自相关矩阵的 $K$ 个最大特征值对应的特征向量张成信号子空间，其余特征向量张成噪声子空间：

$${\mathbf{R}}_{yy}=\sum _{i=1}^K{\lambda }_i{\mathbf{u}}_i{\mathbf{u}}_i^H+\sum _{i=K+1}^N{\lambda }_i{\mathbf{u}}_i{\mathbf{u}}_i^H={\mathbf{U}}_s{\mathbf{\Lambda }}_s{\mathbf{U}}_s^H+{\mathbf{U}}_n{\mathbf{\Lambda }}_n{\mathbf{U}}_n^H$$

信号子空间投影算子为：

$${\mathbf{P}}_s={\mathbf{U}}_s{\mathbf{U}}_s^H$$

信号估计为：

$$\hat{\mathbf{s}}={\mathbf{P}}_s\mathbf{y}$$

在频域中，子空间方法可以通过奇异值分解（SVD）实现。对于信号的时频表示 $Y(t,f)$，其SVD为：

$$Y(t,f)=\sum _{i=1}^r{\sigma }_i{u}_i(t)v_i^{\ast }(f)$$

其中 ${\sigma }_i$ 是奇异值，$u_i(t)$ 和 $v_i(f)$ 分别是时域和频域奇异向量。通过只保留前 $K$ 个奇异值分量，可以实现低秩近似：

$$\hat{S}(t,f)=\sum _{i=1}^K{\sigma }_i{u}_i(t)v_i^{\ast }(f)$$

这种表述在多分辨率分析中特别有用。

非线性自相关滤波

传统的自相关滤波是基于二阶统计量的线性方法。对于非高斯信号或非线性系统，高阶统计量（如三阶矩、四阶矩）包含了更多信息。$n$ 阶累积量是信号统计特性的更完整描述，对于零均值随机变量 $X$，其一维情况下定义为：

$${\kappa }_n=\mathbb{E}[X^n]-\sum _{p\in P(n)}\prod _{B\in p}{\kappa }_{|B|}$$

其中 $P(n)$ 是整数 $n$ 的所有分划集合，$|B|$ 是分划块 $B$ 的元素个数。特别地，二阶累积量等于二阶矩（方差），三阶累积量等于三阶矩，四阶累积量等于四阶矩减去二阶矩的组合：

$${\kappa }_2=\mathbb{E}[X^2]-\mathbb{E}[X{]}^2=\text{Var}[X]$$
$${\kappa }_3=\mathbb{E}[X^3]-3\mathbb{E}[X]\mathbb{E}[X^2]+2\mathbb{E}[X{]}^3$$
$${\kappa }_4=\mathbb{E}[X^4]-4\mathbb{E}[X]\mathbb{E}[X^3]-3\mathbb{E}[X^2{]}^2+12\mathbb{E}[X{]}^2\mathbb{E}[X^2]-6\mathbb{E}[X{]}^4$$

对于多维随机变量，累积量可以推广为张量。$n$ 阶累积量张量的元素为：

$${\kappa }_{i_1,i_2,...,i_n}=\frac{{\partial }^n}{\partial t_{i_1}\partial t_{i_2}...\partial t_{i_n}}\log \Phi (\mathbf{t}){|}_{\mathbf{t}=0}$$

其中 $\Phi (\mathbf{t})$ 是特征函数。高阶自相关函数定义为信号的高阶矩。对于信号 $x(t)$，其 $n$ 阶自相关函数为：

$$R_n^x({\tau }_1,{\tau }_2,...,{\tau }_{n-1})=\mathbb{E}[x(t)x(t+{\tau }_1)x(t+{\tau }_2)...x(t+{\tau }_{n-1})]$$

相应的，$n$ 阶累积量函数定义为：

$$C_n^x({\tau }_1,{\tau }_2,...,{\tau }_{n-1})=\text{cum}[x(t),x(t+{\tau }_1),x(t+{\tau }_2),...,x(t+{\tau }_{n-1})]$$

对于非高斯信号，高阶累积量不为零，包含了比二阶统计量更多的信息。基于高阶统计量的自相关滤波可以捕捉信号中的非线性相关性。非线性自相关滤波的一种形式是基于Volterra级数的滤波器，其输出可以表示为输入的非线性函数：

$$y(t)=h_0+{\int }_{-\infty }^{\infty }{h}_1({\tau }_1)x(t-{\tau }_1)d{\tau }_1+{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{h}_2({\tau }_1,{\tau }_2)x(t-{\tau }_1)x(t-{\tau }_2)d{\tau }_{1}d{\tau }_2+...$$

其中 $h_n({\tau }_1,{\tau }_2,...,{\tau }_n)$ 是 $n$ 阶Volterra核。对于离散信号：

$$y[n]=h_0+\sum _{k_1=0}^{N-1}{h}_1[k_1]x[n-k_1]+\sum _{k_1=0}^{N-1}\sum _{k_2=0}^{N-1}{h}_2[k_1,k_2]x[n-k_1]x[n-k_2]+...$$

Volterra滤波器的设计通常通过最小化均方误差来实现：

$$\underset{h_n}{\min }\mathbb{E}[|s[n]-y[n]{|}^2]$$

这导致了一组复杂的多维Wiener-Hopf方程。

另一种非线性方法是通过核函数将信号映射到高维特征空间，然后在该空间中应用线性滤波：

$$y[n]=\sum _{k=0}^{N-1}{\alpha }_{k}K(x[n],x[k])$$

其中 $K(\cdot ,\cdot )$ 是核函数，如径向基函数（RBF）核：

$$K(x_i,x_j)=\exp \left(-\frac{\Vert x_i-x_j{\Vert }^2}{2{\sigma }^2}\right)$$

这种方法在处理高度非线性系统中特别有效。

自相关滤波的理论限制

自相关滤波虽然强大，但也存在理论限制。理解这些限制有助于正确应用该技术。首先，自相关函数只包含信号的二阶统计特性，对于非高斯过程，这是不完整的描述。信号的完整统计特性要通过高阶累积量或概率密度函数来表征。其次，自相关估计的精度受限于样本数量。对于长度为 $N$ 的信号，自相关估计的方差近似为：

$$\text{Var}[{\hat{R}}_{xx}[k]]\approx \frac{1}{N}\sum _{m=-(N-1)}^{N-1}(1-\frac{|m|}{N})R_{xx}^2[m]$$

这表明，即使对于白噪声，自相关估计的方差也不为零，特别是在延迟 $k$ 较大时。

另一个限制是相位信息的丢失。自相关函数的傅里叶变换是功率谱密度，只包含幅度信息，而没有相位信息。这意味着自相关滤波无法完全重构原始信号，除非有额外的约束（如最小相位假设）。对于维纳滤波器，其性能受限于信噪比。当信噪比很低时，滤波器会大幅衰减信号以抑制噪声，导致过度平滑，丢失信号细节。最后，自相关滤波假设信号和噪声是相互正交的，即：

$$\mathbb{E}[s(t)w(t+\tau )]=0,\forall \tau$$

当这一假设不成立时，滤波性能会下降。

自相关滤波器的计算复杂度

实现自相关滤波时，计算复杂度是一个重要考虑因素。对于长度为 $N$ 的信号，直接计算自相关函数的复杂度为 $O(N^2)$：

$${\hat{R}}_{xx}[k]=\frac{1}{N-|k|}\sum _{n=0}^{N-|k|-1}x[n]x[n+|k|],|k|<N$$

使用快速傅里叶变换（FFT）可以降低复杂度。根据维纳-辛钦定理，自相关函数可以通过功率谱密度的逆傅里叶变换计算：

R ^ x x [ k ] = F − 1 { ∣ F { x [ n ] } ∣ 2 } [ k ] \hat{R}_{xx}[k] = \mathcal{F}^{-1}\{|\mathcal{F}\{x[n]\}|^2\}[k] R^xx​[k]=F−1{∣F{x[n]}∣2}[k]

这种方法的复杂度为 $O(N\log N)$，但可能引入周期性假设导致的边缘效应。对于大规模数据，可以使用分块处理或流处理技术。将信号分成长度为 $L$ 的块，计算每个块的自相关函数，然后合并结果：

$${\hat{R}}_{xx}[k]=\frac{1}{M}\sum _{m=0}^{M-1}{\hat{R}}_{xx}^{(m)}[k]$$

其中 $M$ 是块的数量，${\hat{R}}_{xx}^{(m)}[k]$ 是第 $m$ 个块的自相关估计。对于多维信号，计算复杂度更高。对于 $d$ 维信号，每维长度为 $N$，直接计算自相关函数的复杂度为 $O(N^{2d})$，使用FFT可以降低到 $O(N^d\log N^d)$。在实时应用中，递归计算自相关函数可以进一步减少计算量：

$${\hat{R}}_{xx}^{(n)}[k]=\alpha {\hat{R}}_{xx}^{(n-1)}[k]+(1-\alpha )x[n]x[n-k]$$

其中 $\alpha$ 是平滑因子，$0<\alpha <1$。

频率域自相关滤波的高级分析

在频率域中，自相关滤波可以更深入地分析。对于有色噪声，预白化是一种重要的预处理技术。预白化滤波器的目标是将有色噪声转换为白噪声：

$$W(f)=\frac{1}{\sqrt{S_{ww}(f)}}$$

应用预白化滤波器后，信号的功率谱密度变为：

$$S_{y^{\prime }}(f)=|W(f){|}^2{S}_{yy}(f)=\frac{S_{ss}(f)+S_{ww}(f)}{S_{ww}(f)}=\frac{S_{ss}(f)}{S_{ww}(f)}+1$$

在预白化域中应用维纳滤波器：

$$H_{wiener}^{\prime }(f)=\frac{\frac{S_{ss}(f)}{S_{ww}(f)}}{\frac{S_{ss}(f)}{S_{ww}(f)}+1}=\frac{S_{ss}(f)}{S_{ss}(f)+S_{ww}(f)}$$

预白化和维纳滤波的组合等价于:

$$H_{total}(f)=W(f)\cdot H_{wiener}^{\prime }(f)=\frac{S_{ss}(f)}{S_{ww}(f)\cdot (S_{ss}(f)+S_{ww}(f))}$$

更一般地，我们可以定义频率域中的信噪比：

$$SNR(f)=\frac{S_{ss}(f)}{S_{ww}(f)}$$

维纳滤波器可以表示为SNR的函数：

$$H_{wiener}(f)=\frac{SNR(f)}{SNR(f)+1}$$

当 $SNR(f)\gg 1$ 时，$H_{wiener}(f)\approx 1$；当 $SNR(f)\ll 1$ 时，$H_{wiener}(f)\approx SNR(f)$。

维纳滤波器最小化的均方误差可以表示为：

$$MSE={\int }_{-\infty }^{\infty }\frac{S_{ss}(f)\cdot S_{ww}(f)}{S_{ss}(f)+S_{ww}(f)}df={\int }_{-\infty }^{\infty }\frac{S_{ss}(f)}{SNR(f)+1}df$$

这个表达式显示了滤波误差与信噪比的关系：在信噪比高的频带，误差小；在信噪比低的频带，误差接近信号功率。对于正则化维纳滤波，引入一个正则化参数 $\gamma$：

$$H_{reg}(f)=\frac{S_{ss}(f)}{S_{ss}(f)+S_{ww}(f)+\gamma }$$

这等价于向噪声谱中添加一个常数项，可以防止在低信噪比区域滤波器增益过大，提高数值稳定性。