75、非平稳窄带声学干扰的主动抑制

非平稳窄带声学干扰的主动抑制

在声学系统中，主动噪声控制（ANC）技术对于减少窄带干扰具有重要意义。为了实现有效的噪声抑制，需要解决系统参数识别和控制器设计等问题。本文将探讨相关的方法和策略。

1. 二次路径识别问题及应用领域

在ANC系统中，二次路径的识别是一个关键问题。目前，最流行的二次路径识别方法是辅助噪声法。该方法通过向控制回路引入白噪声形式的附加信号来激励系统。然而，这种方法存在明显的缺点：
- 系统中的干扰总体水平会增加，这与控制目标相冲突。
- 基本形式的辅助噪声法收敛速度较慢，而改进方法会显著增加控制器的复杂度和计算成本。

ANC系统有许多重要的商业应用：
- 主动耳机 ：由于体积小，通常采用反馈控制，高频噪声通过被动方式消除。
- 主动降噪摩托车头盔 ：能显著降低噪声。
- ANC抗打鼾系统 ：由于系统组件放置自由度受限以及睡眠者头部移动导致二次路径变化，该应用具有挑战性。
- 医疗设备噪声消除 ：包括改善患者舒适度的婴儿培养箱ANC系统，以及减少噪声对医护人员长期影响的MRI噪声主动控制。

2. 问题描述

为了简化数学分析，假设所有信号（如控制信号、参考信号、窄带干扰和测量噪声）均为复数值。考虑离散时间系统输出处窄带干扰的减少问题，该系统由以下方程描述：
[y(t) = K_p(q^{-1})u(t - 1) + d(t) + v(t)]
其中：
- (t \in Z = {…, -1, 0, 1, …}) 表示归一化离散时间。
- (y(t)) 表示受干扰的复数值输出信号。
- (q^{-1}) 是后向移位算子。
- (d(t) = a(t)e^{j\varphi(t)})，(\varphi(t) = \sum_{i = 1}^{t}\omega(i)) 是具有缓慢变化幅度 (a(t)) 和缓慢变化角频率 (\omega(t) \in (-\pi, \pi]) 的复数值正弦干扰。
- (v(t)) 是满足假设 (A1) 的复数值宽带噪声，即 ({v(t)}) 是均值为零、方差为 (\sigma_v^2) 的独立复数值随机变量序列，且 (v(t)) 的实部和虚部相互独立。
- (K_p(q^{-1})) 是满足条件 (A2) 的稳定、线性、单输入单输出系统的传递函数，即 (K_p(e^{-j\omega}) \neq 0)，(\forall\omega \in (-\pi, \pi])。

目标是设计一个最小方差控制器，即找到使ANC系统输出信号的均方值最小的控制规则：
[u(t) : E[|y(t)|^2] \to \min]

将考虑纯反馈控制策略和结合反馈控制与前馈控制（如果参考信号可用）的混合策略。

3. 反馈控制策略

3.1 系统和干扰完全先验已知的情况

假设系统的传递函数 (K_p(q^{-1})) 已知，且干扰形式为 (d(t) = a_0e^{j\omega_0t})（即 (a(t) \equiv a_0) 且 (\omega(t) \equiv \omega_0)）且可测量。为了补偿这种形式的正弦干扰，应生成一个正弦控制信号 (u(t))，使得系统输出与干扰幅度相同但极性相反。根据线性系统对正弦输入的缩放和移位特性，可得到稳态近似：
[K_p(q^{-1})u(t - 1) = k_pu(t - 1)]
其中 (k_p = K_p(e^{-j\omega_0}) \in C) 是系统在频率 (\omega_0) 处的复数值增益/衰减。结合上述关系，可得控制信号的形式为：
[u(t) = -\frac{d(t + 1)}{k_p}]
在该控制规则下，开环系统的稳态输出信号为 (y(t) = v(t))，其方差达到最小值 (\lim_{t \to \infty}E[|y(t)|^2] = \sigma_v^2)。

3.2 系统完全先验已知，干扰部分先验已知（已知频率）的情况

假设干扰的频率 (\omega_0) 恒定且已知，但其幅度根据随机游走模型缓慢变化：
[a(t + 1) = a(t) + \eta(t + 1)]
其中 ({\eta(t)}) 是满足假设 (A3) 的独立随机变量序列，即均值为零、方差为 (\sigma_{\eta}^2)，且与 ({v(t)}) 独立。在这些假设下，干扰信号可表示为递归形式：
[d(t + 1) = e^{j\omega_0}d(t) + \tilde{\eta}(t + 1)]
其中 ({\tilde{\eta}(t)})，(\tilde{\eta}(t) = e^{j\omega_0}\eta(t)) 是均值为零、方差为 (\tilde{\sigma} {\eta}^2 = \sigma {\eta}^2) 的独立复数值随机变量序列。

使用近似 (K_p(q^{-1})u(t - 1) \approx k_pu(t - 1))，可得到最优最小方差控制规则的稳态形式：
[\hat{d}(t + 1|t) = e^{j\omega_0}[\hat{d}(t|t - 1) + \mu_ay(t)]]
[u(t) = -\frac{\hat{d}(t + 1|t)}{k_p}]
其中 (\hat{d}(t + 1|t)) 是干扰的预测值，(\mu_a \in R) 是适应增益，(\mu_a = -\frac{\xi}{2} + \sqrt{\frac{\xi^2}{4} + \xi})，(\xi = \frac{\sigma_{\eta}^2}{\sigma_v^2})。在该控制规则下，系统输出为 (y(t) \approx c(t) + v(t))，其中 (c(t)) 是抵消误差：
[c(t) = d(t) - \hat{d}(t|t - 1) = e^{j\omega_0}(1 - \mu_a)c(t - 1) + \tilde{\eta}(t) - \mu_ae^{j\omega_0}v(t + 1)]
稳态输出信号的方差为 (\lim_{t \to \infty}var[y(t)] = \sigma_c^2 + \sigma_v^2)，其中 (\sigma_c^2 = \lim_{t \to \infty}var[c(t)] = \frac{\sigma_{\eta}^2 + \sqrt{\sigma_{\eta}^4 + 4\sigma_{\eta}^2\sigma_v^2}}{2})。

3.3 系统和干扰部分先验已知（已知频率）的情况

在实际情况中，系统增益 (k_p) 和适应增益 (\mu_a) 不仅未知，而且可能随时间变化。考虑将真实系统增益 (k_p) 替换为“标称”增益 (k_n \neq k_p)，最优适应增益 (\mu_a) 替换为任意选择的增益 (\mu \neq \mu_a) 的情况，得到控制器：
[\hat{d}(t + 1|t) = e^{j\omega_0}[\hat{d}(t|t - 1) + \mu y(t)]]
[u(t) = -\frac{\hat{d}(t + 1|t)}{k_n}]
定义系统增益建模误差 (\beta \in C)，(\beta = \frac{k_p}{k_n} \neq 1)。输出信号可表示为 (y(t) = c(t) + v(t))，其中抵消误差为：
[c(t) = d(t) - \beta\hat{d}(t|t - 1) = e^{j\omega_0}(1 - \mu\beta)c(t - 1) + \tilde{\eta}(t) - \mu\beta e^{j\omega_0}v(t + 1)]
比较可知，当适应增益系数 (\mu) 满足 (\mu\beta = \mu_a) 时，抵消误差与最优控制器产生的误差相同。这意味着通过适当选择 (\mu) 的值，可以“补偿”由于采用不适当的系统增益 (k_n) 而导致的建模误差。

为了在线调整复数值适应增益 (\mu)，设计了自适应算法SONIC（Self-Optimizing Narrowband Interference Canceller）：
[z(t) = e^{j\omega_0}[(1 - c_{\mu})z(t - 1) - \frac{c_{\mu}}{\hat{\mu}(t - 1)}y(t - 1)]]
[r(t) = \rho r(t - 1) + |z(t)|^2]
[\hat{\mu}(t) = \hat{\mu}(t - 1) - \frac{y(t)z^ (t)}{r(t)}]
[\hat{d}(t + 1|t) = e^{j\omega_0}[\hat{d}(t|t - 1) + \hat{\mu}(t)y(t)]]
[u(t) = -\frac{\hat{d}(t + 1|t)}{k_n}]
其中 (c_{\mu} \in [0.005, 0.05]) 是一个小的正常数，(z^ (t)) 是 (z(t)) 的复共轭。变量 (z(t)) 由 (z(t) = \frac{\partial y(t; \mu)}{\partial \mu}) 给出，这里使用了Wirtinger导数。

SONIC算法的自优化循环不仅可以补偿系统建模误差，还可以优化控制器。在假设 (A1) - (A3) 下，(\lim_{t \to \infty}E[\hat{\mu}(t)] \approx \frac{\mu_a}{\beta})，这意味着SONIC控制器在均值意义上收敛到最优控制器。此外，由于使用递归预测误差（RPE）方法优化 (\mu)，该控制器在假设 (A1) - (A3) 不成立的情况下也能尝试最小化抵消误差，对系统动态变化和信噪比变化等具有较好的响应。

通过一个简单的模拟实验可以证明SONIC的鲁棒性。实验中系统动态发生切换，窄带干扰和宽带噪声的参数设置为 (\sigma_v = 0.1)，(\sigma_e = 0.001)，(\omega_0 = 0.1)，(d(0) = 1)。SONIC算法的参数设置为 (c_{\mu} = 0.005)，(\rho = 0.9995)，(k_n = e^{j\omega_0})，并配备了安全保护机制。实验结果表明，算法能够很好地应对初始收敛问题和系统动态的突然变化，具有自稳定特性。

系统切换时间表和相应的建模误差如下表所示：
| 时间间隔 | 系统 | (|\beta|) | (Arg\beta[^{\circ}]) |
| ---- | ---- | ---- | ---- |
| (0 < t < 15000) | (K_1(z) = \frac{0.0952}{1 - 0.9048z^{-1}}) | (0.708) | (-47.9) |
| (15000 \leq t < 30000) | (K_2(z) = \frac{0.0238}{1 - 0.9762z^{-1}}) | (0.234) | (-79.3) |
| (30000 \leq t < 45000) | (K_3(z) = \frac{0.2}{1 - 0.8z^{-1}}) | (0.913) | (-27.1) |
| (45000 \leq t \leq 60000) | (K_4(z) = \frac{+0.1 - 0.14z^{-1}}{1 - 1.8391z^{-1} + 0.8649z^{-2}}) | (1.960) | (121.1) |

以下是SONIC算法的流程图：

graph TD;
    A[初始化: \(\hat{d}(0), r(0), z(0), \hat{\mu}(0)\)] --> B[计算 \(z(t)\)];
    B --> C[计算 \(r(t)\)];
    C --> D[计算 \(\hat{\mu}(t)\)];
    D --> E[计算 \(\hat{d}(t + 1|t)\)];
    E --> F[计算 \(u(t)\)];
    F --> G[更新系统状态];
    G --> B;

3.4 系统和干扰部分先验已知（未知频率）的情况

假设窄带干扰的频率 (\omega(t)) 未知且随时间缓慢变化。将已知频率 (\omega_0) 替换为瞬时频率的当前估计值，得到修改后的算法：
[\hat{d}(t + 1|t) = e^{j\hat{\omega}(t|t - 1)}[\hat{d}(t|t - 1) + \mu y(t)]]
[\hat{\omega}(t + 1|t) = g[Y(t)]]
[u(t) = -\frac{\hat{d}(t + 1|t)}{k_n}]
其中 (\hat{\omega}(t + 1|t)) 是基于截至时刻 (t) 的观测值 (Y(t) = {y(i), i \leq t}) 对瞬时频率 (\omega(t + 1)) 的一步预测。

瞬时频率的估计可以使用简单的梯度算法：
[\hat{\omega}(t + 1|t) = \hat{\omega}(t|t - 1) + \gamma Arg\left(\frac{\hat{d}(t + 1|t)e^{-j\hat{\omega}(t|t - 1)}}{\hat{d}(t|t - 1)}\right)]
在缓慢时变干扰信号的情况下，可近似为：
[\hat{\omega}(t + 1|t) = \hat{\omega}(t|t - 1) + \gamma Im\left(\frac{\mu y(t)}{\hat{d}(t|t - 1)}\right)]
得到控制算法：
[\hat{d}(t + 1|t) = e^{j\hat{\omega}(t|t - 1)}[\hat{d}(t|t - 1) + \mu y(t)]]
[\hat{\omega}(t + 1|t) = \hat{\omega}(t|t - 1) + \gamma Im\left(\frac{\mu y(t)}{\hat{d}(t|t - 1)}\right)]
[u(t) = -\frac{\hat{d}(t + 1|t)}{k_n}]

使用局部平均法和近似线性滤波（ALF）技术分析该算法的跟踪特性。在一定假设下，系统输出可近似表示为 (y(t) \approx d(t) - \beta\hat{d}(t|t - 1) + v(t))，其中 (\beta = \frac{k_T}{k_n})，(k_T) 是系统在时间区间 (T) 内的平均增益。通过ALF技术得到线性关系：
[\Delta\hat{x}(t) = (1 - \mu\beta)\Delta\hat{x}(t - 1) + ja_0^2\Delta\hat{\omega}(t - 1)\mu\beta s(t - 1)]
[\Delta\hat{\omega}(t + 1) = \Delta\hat{\omega}(t) - \frac{\gamma}{a_0^2}Im[\mu\beta\Delta\hat{x}(t)] - \frac{\gamma}{a_0^2}Im[\mu\beta z(t)] + w(t + 1)]
其中 (\Delta x(t) = c(t)d^ (t))，(s(t) = v(t)d^ (t))。

当瞬时频率根据随机游走模型演化时，假设真实系统增益已知（(\beta = 1)），可得到频率跟踪误差的稳态方差表达式：
[E[(\Delta\hat{\omega}(t))^2] \approx \frac{\gamma^2\mu}{4a_0^2}\sigma_v^2 + \left(\frac{1}{2\gamma} + \frac{1}{2\mu}\right)\sigma_w^2]
通过计算可得到使均方误差最小的系数 (\mu_{\omega}) 和 (\gamma_{\omega})，且最优调整后的算法在统计上是有效的，其性能与高分辨率但计算复杂度高的ESPRIT算法相当。

结合前面描述的所有自适应机制，得到扩展版的SONIC算法：
- 自优化 ：
[z(t) = e^{j\hat{\omega}(t|t - 1)}[(1 - c_{\mu})z(t - 1) - \frac{c_{\mu}}{\hat{\mu}(t - 1)}y(t - 1)]]
[p(t) = \rho r(t - 1) + |z(t)|^2]
[\hat{\mu}(t) = \hat{\mu}(t - 1) - \frac{y(t)z^ (t)}{p(t)}]
- 预测控制 ：
[\hat{d}(t + 1|t) = e^{j\hat{\omega}(t|t - 1)}[\hat{d}(t|t - 1) + \hat{\mu}(t)y(t)]]
[u(t) = -\frac{\hat{d}(t + 1|t)}{k_n[\hat{\omega}(t|t - 1)]}]
- 频率估计 *：
[\hat{\omega}(t + 1|t) = \hat{\omega}(t|t - 1) + \gamma Im\left(\frac{\hat{\mu}(t)y(t)}{\hat{d}(t|t - 1)}\right)]

在一个实际的实值系统实验中，扩展版SONIC算法在修改后在普通PC上实现。采样率为1 kHz，左扬声器产生干扰，右扬声器产生反声，误差麦克风放置在距右扬声器15 cm处。干扰的瞬时频率在241 Hz到250 Hz之间以20秒为周期缓慢正弦变化。系统的标称增益设置为 (k_n = 1)，SONIC的其他设置为 (c_{\mu} = 0.01)，(\mu_{max} = 0.05)，(\gamma = 0.01)。应用扩展版SONIC算法可使谐波噪声衰减15 - 20 dB。

综上所述，通过不同的反馈控制策略和自适应算法，可以有效地抑制非平稳窄带声学干扰，提高ANC系统的性能。这些算法在不同的先验知识条件下都能发挥作用，并且具有一定的鲁棒性和自适应性。在实际应用中，可以根据具体情况选择合适的算法和参数设置。

非平稳窄带声学干扰的主动抑制

4. 算法对比与优势分析

在主动噪声控制（ANC）领域，不同的算法有着各自的特点和适用场景。下面将SONIC算法与其他常见算法进行对比，以突出其优势。

算法	优点	缺点	适用场景
辅助噪声法	原理简单，能对二次路径进行识别	增加系统干扰总体水平，收敛慢，改进后复杂度和计算成本高	对干扰增加不太敏感，对计算资源要求不高的场景
ESPRIT算法	频率估计分辨率高	计算复杂度高	对频率估计精度要求极高，且有足够计算资源的场景
SONIC算法	能补偿系统建模误差，具有自稳定特性，统计效率高，对系统动态和信噪比变化响应好	-	非平稳窄带声学干扰场景，尤其是系统参数和干扰特性变化的情况

SONIC算法的优势主要体现在以下几个方面：
- 自适应性强 ：能够根据系统动态和干扰特性的变化自动调整适应增益，有效应对系统建模误差和干扰的非平稳性。
- 计算效率高 ：相比于ESPRIT算法等复杂的频率估计方法，SONIC算法的计算复杂度较低，更适合实时应用。
- 鲁棒性好 ：在系统发生突然变化或假设条件不成立的情况下，仍能保持较好的控制性能，具有自稳定特性。

5. 实际应用中的注意事项

在实际应用ANC系统和SONIC算法时，需要注意以下几点：
1. 参数设置 ：
- (c_{\mu})：取值范围在 ([0.005, 0.05]) 之间，需要根据系统的具体情况进行调整。较小的值可以使算法更加稳定，但收敛速度可能较慢；较大的值可以加快收敛速度，但可能会导致系统不稳定。
- (\rho)：遗忘常数，取值在 ([0.999, 0.9999]) 之间，决定了局部分析窗口的有效宽度。较大的值可以使算法对历史数据的依赖更强，更能跟踪系统的缓慢变化；较小的值则更关注近期数据，对系统的快速变化响应更灵敏。
- (\gamma)：在频率估计中，适应增益 (0 < \gamma \ll 1)，需要根据干扰频率的变化速度进行调整。较大的值可以加快频率跟踪速度，但可能会引入较大的估计误差；较小的值则可以提高估计的稳定性，但跟踪速度较慢。
2. 安全保护机制 ：
- 为了防止算法出现不稳定或异常情况，需要设置安全保护机制，如 (\mu_{max})、(\mu_{min})、(\Delta\mu_{max}(t))、(r_{max}) 和 (r_{min}) 等参数。这些参数可以限制适应增益和其他变量的取值范围，保证算法的正常运行。
3. 初始条件 ：
- 合理的初始条件对于算法的初始收敛至关重要。例如，在SONIC算法中，(\hat{d}(0))、(r(0))、(z(0)) 和 (\hat{\mu}(0)) 的取值会影响算法的初始性能。通常可以根据系统的先验知识或经验来选择合适的初始值。

6. 未来发展趋势

随着科技的不断进步，ANC技术在非平稳窄带声学干扰抑制方面有着广阔的发展前景。以下是一些可能的未来发展趋势：
1. 多传感器融合 ：结合多个传感器的信息，如麦克风阵列、加速度计等，可以更全面地获取声学环境信息，提高干扰抑制的效果。
2. 深度学习应用 ：利用深度学习算法，如神经网络、卷积神经网络等，对声学信号进行特征提取和建模，从而实现更智能、更高效的干扰抑制。
3. 与其他技术的结合 ：将ANC技术与无线通信、物联网等技术相结合，实现远程控制和智能化管理，拓展其应用领域。

7. 总结

非平稳窄带声学干扰的主动抑制是一个具有挑战性但又非常重要的研究领域。通过对二次路径识别问题的研究和反馈控制策略的设计，我们可以有效地减少声学系统中的窄带干扰。SONIC算法作为一种自适应控制算法，在补偿系统建模误差、应对干扰的非平稳性和提高控制性能方面表现出色。

在实际应用中，需要根据具体情况选择合适的算法和参数设置，并注意参数调整、安全保护机制和初始条件等问题。未来，ANC技术有望与多传感器融合、深度学习等技术相结合，进一步提高干扰抑制的效果和应用范围。

以下是SONIC算法在实际应用中的流程图：

graph TD;
    A[系统初始化] --> B[采集声学信号];
    B --> C[计算瞬时频率估计 \(\hat{\omega}(t|t - 1)\)];
    C --> D[自优化计算 \(z(t)\)、\(p(t)\)、\(\hat{\mu}(t)\)];
    D --> E[预测控制计算 \(\hat{d}(t + 1|t)\)、\(u(t)\)];
    E --> F[生成控制信号并输出到扬声器];
    F --> G[采集输出信号 \(y(t)\)];
    G --> B;

通过以上的分析和总结，我们可以更好地理解和应用ANC技术和SONIC算法，为解决非平稳窄带声学干扰问题提供有效的方法和思路。