Matlab：使用 Filtered-x LMS FIR 自适应滤波器实现有源噪声控制

• 🏠 Matlab 官方手册： Active Noise Control Using a Filtered-X LMS FIR Adaptive Filter - MATLAB & Simulink - MathWorks 中国。

一、有源噪声控制 (Active Noise Control)

• 有源噪声控制通过电声系统，使用监测传感器（比如麦克风）与输出制动器（比如扬声器）来减少不希望的噪声音量在空气中传播。有源噪声控制系统的目的是生成一种“抗噪声”，使用自适应滤波器在期望的安静区域中衰减不需要的噪声，原理图如下所示：

二、次级传播路径的实现代码 (Secondary Propagation Path)

🚀 该问题与传统自适应噪声消除不同之处在于：
我们无法直接测量期望的响应信号，只有噪声源信号可以被监测到。有源噪声控制系统在自适应过程中必须考虑次级通道(Secondary Path)中的传播。

次级传播路径是抗噪声从输出扬声器到安静区内的误差麦克风的路径。下面的Matlab指令可以生成一个扬声器到误差麦克风的脉冲响应，即次级通道的脉冲响应。其带宽限制在160-2000Hz的范围内，且滤波器长度为0.1秒，另外，采样频率为8000Hz。

🦄 通过Matlab生成一个固定频带宽度的次级通道的具体步骤如下所示：

• 1、首先设置次级通道相关超参数： 采样频率(Fs)、滤波器长度(N)、频率上 / 下边沿(Fhigh / Flow)、信号传播延迟(delayS)、阻带衰减(Ast)、滤波器阶数(Nfilt)；
• 2、然后根据次级通道相关超参数设计带通滤波器， 主要包括次级通道带通滤波器设计与次级通道脉冲响应获取两个流程：
  • （1）次级通道带通滤波器设计： 将滤波器阶数(Nfilt)、频率上 / 下边沿(Fhigh / Flow)、阻带衰减(Ast)、采样频率(Fs)输入Matlab带通滤波器设计类方法fdesign.bandpass，得到次级通道的带通滤波器模板filtSpecs；使用design类方法将滤波器模板filtSpecs设计成一个具体的次级通道带通滤波器bandpass。其中bandpass的设计方法为Chebyshev II(cheby2)型滤波器，滤波器结构为直接II型 SOS 二阶块(df2tsos, Transposed Direct From II, Second-Order Secitons)。
  • （2）次级通道脉冲响应获取： 构造输入次级通道的参考信号，将其输入次级通道带通滤波器得到次级通道的脉冲响应，并进行正则化处理。

%%%%%%%%%%%%%%%% 生成并绘制第二通道的脉冲响应 %%%%%%%%%%%%%%%%

% 超参数设置
Fs       = 8e3;     % 采样频率为 8kHz
N         = 800;    % 采样点为800个，对应于时域中的 0.1s
Flow    = 160;    	% 频带下边沿160Hz
Fhigh  = 2000;  	% 频带上边沿2000Hz
delayS = 7;
Ast      = 20;      % 阻带衰减 (Stopband Attenuation) 为 20dB
Nfilt     = 8;      % 滤波器阶数

% 使用带通滤波器方法生成具有固定带宽次级通道的脉冲响应
% 得到次级通道的带通滤波器bandpass
filtSpecs = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
bandpass = design(filtSpecs,'cheby2','FilterStructure','df2tsos', ...
    'SystemObject',true);

% 将参考信号输入带通滤波器得到次级通道的脉冲响应secondaryPathCoeffsActual
secondaryPathCoeffsActual = bandpass([zeros(delayS,1); ...
                       log(0.99*rand(N-delayS,1)+0.01).* ...
                       sign(randn(N-delayS,1)).*exp(-0.01*(1:N-delayS)')]);
secondaryPathCoeffsActual = ...
    secondaryPathCoeffsActual/norm(secondaryPathCoeffsActual);

% 绘图
t = (1:N)/Fs;
plot(t, secondaryPathCoeffsActual, ...      % 需要绘制的x,y数据
    'Color', [0.4940 0.1840 0.5560], ...    % 线条颜色
    'LineStyle', '--', ...     % 线条类型
    'LineWidth', 2, ...        % 线条宽度
    'Marker', 'square', ...    % 数据点形状
    'MarkerIndices', 1:3:length(secondaryPathCoeffsActual), ...  % 每隔3个点绘制一个数据点
    'MarkerEdgeColor', [0.1, 0.1, 0.1], ...    % 数据点边缘颜色
    'MarkerFaceColor', [0.8, 0.8, 0.8], ...    % 数据点面颜色
    'MarkerSize', 5, ...      % 数据点大小
    'LineJoin', 'round')      % 线条转折点的形状
    
xlabel('Time [sec]', ...
    'FontSize', 16, ...
    'FontWeight', 'bold');
ylabel('Coefficient value', ...
    'FontSize', 16, ...
    'FontWeight','bold');
title('True Secondary Path Impulse Response', ...
    'FontSize', 18, ...
    'FontWeight', 'bold');

set(gca,'FontSize',14);     % 设置坐标刻度字体大小

三、次级通道 (Secondary Path) 传播路径的估计

🚀 有源噪声控制的第一任务：估计次级传播路径的脉冲响应。
对次级通道传播路径的脉冲响应通常是在噪声达到误差麦克风之前进行的，评估次级通道传播路径脉冲响应的思路为：

• (1) 通过上面得到的次级通道脉冲响应secondaryPathCoeffsActual作为分子(Numerator)，通过dsp.FIRFilter方法得到次级通道信号生成器FIR滤波器secondaryPathGenerator；
• (2) 构造一个长度为30000的随机信号randomSignal。将randomSignal加上麦克风处的随机干扰信号0.01*randn(30000, 1)构成合成信号输入次级通道FIR滤波器secondaryPathGenerator到次，最终得到次级通道的估计输出信号secondaryPathMeasured。

下面的代码产生了3.75秒的随机噪声以及误差麦克风的测量信号：

%%%%%%%%%%%%%%%% 次级通道传播的脉冲估计 %%%%%%%%%%%%%%%%

ntrS = 30000;       % 信号长度，由于采样频率为800，所以信号的时间长度3.75s
randomSignal = randn(ntrS, 1);    % 合成的随机信号

% 构造次级通道的滤波器secondaryPathGenerator
secondaryPathGenerator = dsp.FIRFilter('Numerator', secondaryPathCoeffsActual.');
% 将合成信号收入次级通道滤波器secondaryPathGenerator，得到次级通道
% 的估计输出信号secondaryPathMeasured
secondaryPathMeasured = secondaryPathGenerator(randomSignal) + ... % 次级通道中传播的随机噪声
    0.01*randn(ntrS, 1);  % 误差麦克风处监测到的噪声

四、基于NLMS自适应滤波器的次级通道估计

实际情况下，次级通道的滤波器长度可以比实际通过其噪声信号的长度短，并且在通常情况下，并不需要对整个噪声信号进行自适应控制。这里使用250个信号点(taps)长度作为次级通道滤波器的长度，其对应于31ms的脉冲响应长度。

虽然，很多自适应FIR滤波器都可以实现有源噪声控制，但是由于归一化最小均方差(NLMS, Normalized Least Mean Square)算法简单且具有良好的鲁棒性，本文使用NLMS进行有源噪声控制。

自适应滤波器是利用前一时刻滤波器参数，自动调整当前时刻滤波器参数，以适应信号随时间变化的统计特征，从而实现最优滤波，其原理如下图所示：

• $x(n)$为输入信号；
• $d(n)$为期望信号；
• $y(n)$为输入信号通过参数可调数字滤波器（自适应滤波器）的输出信号；
• $e(n)$为误差信号。

🚀 基于NLMS自适应滤波器的次级通道估计的具体步骤为：

• （1）通过dsp.LMSFilter方法构造次级通道估计器secondaryPathEstimator。
• （2）输入随机信号randomSignal与次级通道的估计输出信号secondaryPathMeasured输入次级通道估计器secondaryPathEstimator得到期望信号(用于消除噪声的信号)secondaryPathCoeffsEst、NLMS自适应滤波器的输出信号yS与误差信号eS。

%%%%%%%%%%%%%%%% 三、基于NLMS的次级通道估计 %%%%%%%%%%%%%%%%

M = 250;      % 滤波器长度
muS = 0.1;   % 步长

% 构造基于NLMS滤波器的次级通道估计器secondaryPathEstimator
secondaryPathEstimator = dsp.LMSFilter('Method', 'Normalized LMS', ...  % 使用NLMS方法
    'StepSize', muS, ...    % 步长
    'Length', M);             % 滤波器长度

% 将输入信号randomSignal与上面得到的次级通道输出secondaryPathMeasured
% 输入NLMS自适应滤波器得到：
% （1）secondaryPathCoeffsEst：NLMS自适应滤波器权重系数
% （2）yS：输入信号通过NLMS自适应滤波器的次级通道输出
% （3）eS：误差麦克风的信号
[yS, eS, SecondaryPathCoeffsEst] = secondaryPathEstimator(randomSignal, secondaryPathMeasured);

% 绘制结果
n = 1:ntrS;
figure, plot(n, secondaryPathMeasured, n, yS, n, eS);
xlabel("Number of iterations", "FontSize", 16, "FontWeight", "bold");
ylabel("Signal value", "FontSize", 16, "FontWeight", "bold");
title("Secondary Identification Using the NLMS Adaptive Filter", "FontSize", 18, "FontWeight", "bold")
set(gca,'FontSize',14);     % 设置坐标刻度字体大小

由上图可以看出，NLMS算法在大约10000次迭代后收敛。

五、基于NLMS的次级通道估计的准确性

下面的代码实现展示次级路径脉冲响应估计的准确度情况。

%%%%%%%%%%%%%%%% 基于NLMS的次级通道估计的准确性 %%%%%%%%%%%%%%%%

figure;
% 次级通道的真实响应
plot(t, secondaryPathCoeffsActual, "LineStyle","--", "LineWidth", 1.5, "Color", "k"), hold on;
% 次级通道的估计响应
plot(t(1:M), SecondaryPathCoeffsEst, "LineStyle","-.", "LineWidth", 1.5), hold on;
% 真实响应与估计响应的误差
plot(t, [secondaryPathCoeffsActual(1:M)-SecondaryPathCoeffsEst(1:M); ...
    secondaryPathCoeffsActual(M+1:N)], "LineStyle","-", "LineWidth", 1.5)

xlabel("Time [sec]", "FontSize", 16, "FontWeight", "bold")
ylabel("Coefficient value", "FontSize", 16, "FontWeight", "bold")
title("Secondary Path Impulse Response Estimation", "FontSize", 18, "FontWeight", "bold")
legend('True', 'Estimated', 'Error');
set(gca,'FontSize',14); 

六、主传播通道(Primary Propagation Path)

将要被消除的噪声在主传播通道中的传播也可以由一个线性滤波器刻画。下面的指令生成一个从输入到误差麦克风的一个频率介于200-800Hz并且滤波器长度为0.1s的脉冲响应。

%%%%%%%%%%%%%%%% 主传播通道 %%%%%%%%%%%%%%%%
delayW = 15;        % 主通道的信号延迟
Flow     = 200;     % 频带下边界
Fhigh   = 800;      % 频带上边界
Ast       = 20;        % 20 dB阻带衰减
Nfilt      = 10;       % 滤波器阶数

% 设计带通滤波器bandpass2，以产生200-800频带内的脉冲响应
filtSpecs2 = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
bandpass2 = design(filtSpecs2,'cheby2','FilterStructure','df2tsos', ...
    'SystemObject',true);

% 将噪声输入带通滤波器bandpass2以生成脉冲响应
primaryPathCoeffs = bandpass2([zeros(delayW, 1); ...
    log(0.99*rand(N-delayW, 1)+0.01).*sign(randn(N-delayW, 1)).*exp(-0.01*(1:N-delayW)')]);
primaryPathCoeffs = primaryPathCoeffs / norm(primaryPathCoeffs);

% 绘制主通道的脉冲响应
figure, plot(t, primaryPathCoeffs, "Color", 'b', 'LineStyle', '-.', 'LineWidth', 1.5, ...
    'Marker', 'o', 'MarkerIndices', 1:3:length(primaryPathCoeffs), ...
    'MarkerEdgeColor', [0.1, 0.1, 0.1], 'MarkerFaceColor', [0.8, 0.8, 0.8], 'MarkerSize', 5)
    
xlabel('Time [sec]', ...
    'FontSize', 16, ...
    'FontWeight', 'bold');
ylabel('Coefficient value', ...
    'FontSize', 16, ...
    'FontWeight','bold');
title('Primary Path Impulse Response', ...
    'FontSize', 18, ...
    'FontWeight', 'bold');

set(gca,'FontSize',14);     % 设置坐标刻度字体大小

七、有源噪声控制(ANC, Active Noise Control)

7.1 有源噪声的初始化

在有源噪声控制中最流行的自适应算法为filterd-X LMS算法。它使用次级通道的估计结果计算一个输出信号，该输出在误差传感器处与不期望的噪声相互抵消，以达到噪声消除的目的。参考信号(reference signal)是早噪声源附近进行监测得到的。下面的代码使用大约44ms及步长为0.0001的滤波器控制器：

% 根据主通道脉的冲响应结果primaryPathCoeffs，构造FIR滤波器以便构造主通道传播模型（传递函数）
primaryPathGenerator = dsp.FIRFilter('Numerator', primaryPathCoeffs.');

% 使用Filtered-X LMS自适应滤波器对噪声进行控制
L = 350;
muW = 0.0001;
noiseController = dsp.FilteredXLMSFilter('Length',L,'StepSize',muW, ...
    'SecondaryPathCoefficients',SecondaryPathCoeffsEst);

% 构造Sine波形的构造器合成地生成噪声
A  = [.01 .01 .02 .2 .3 .4 .3 .2 .1 .07 .02 .01];
La = length(A);
F0 = 60;
k  = 1:La;
F  = F0*k;
phase = rand(1, La); % 随机初始化相位
sine = audioOscillator('NumTones', La, 'Amplitude',A,'Frequency',F, ...
    'PhaseOffset',phase,'SamplesPerFrame',512,'SampleRate',Fs);

% 使用音频播放器分别播放消去噪声前后的输入噪声信号
player = audioDeviceWriter('SampleRate', Fs);

% 使用频谱分析器展示原始信号与消去噪声的信号
scope = spectrumAnalyzer("SampleRate", Fs, "OverlapPercent", 80, ...
    "PlotAsTwoSidedSpectrum", false, ...
    "ShowLegend", true, ...
    "ChannelNames", {'Original noise signal', 'Attenuated noise'});

7.2 基于Filtered-X LMS算法的有源噪声的仿真

下面进行有源噪声控制系统的设计，为了突出有源噪声控制系统的影响，在开始200步迭代中不使用有源噪声控制，并在消噪之前监听误差麦克风处的信号。

一旦自适应滤波器开始工作，大概经过6-10s的仿真时间算法开始收敛，即噪声被有效消除。通过对比残差信号的频谱，我们可以看到输入噪声在很大程度上得到了衰减。

%%%%%%%%%%%%%%%% 基于Filtered-X LMS算法的有源噪声的仿真 %%%%%%%%%%%%%%%%

for m = 1:400
    % 通过上面的sin()生成带有随机相位的合成噪声
    x = sine();
    d = primaryPathGenerator(x) + ... % 输入噪声在主通道中的传播
        0.1*randn(size(x)); % 添加测量噪声

    % 开始迭代
    if m <= 200
        % 200步迭代之前，不进行有源噪声控制
        e = d;
    else
        % 200步迭代之后，添加有源噪声控制
        xhat = x + 0.1*randn(size(x));
        [y, e] = noiseController(xhat, d);
    end
    
    player(e);      % 播放噪声信号
    scope([d, e])   % 展示原始信号的频率(Channel 1)与衰减的噪声(通道2)
end

release(player);    % 释放播放器对象
release(scope);     % 释放频谱显示对象

代码执行结果如下图所示：