使用FX-LMS FIR 自适应滤波器的主动噪声控制

远日点613

已于 2025-02-24 14:38:29 修改

阅读量145

点赞数

文章标签： matlab 信号处理

于 2025-02-23 22:20:32 首次发布

原文链接：https://ww2.mathworks.cn/help/audio/ug/active-noise-control-using-a-filtered-x-lms-fir-adaptive-filter.html

版权

翻译来自MATLAB文档：Active Noise Control Using a Filtered-X LMS FIR Adaptive Filter

若官方有中文版了可以提醒删除。

主动噪声控制

在主动噪声控制中，目标是通过使用电声系统减少空气中传播的不需要的噪声音量，该系统利用测量传感器（如麦克风）和输出执行器（如扬声器）。噪声信号通常来自某种设备，例如旋转机械，因此可以在其源附近测量噪声。主动噪声控制系统的目标是产生一种“反噪声”，以在所需的安静区域内衰减不需要的噪声，使用自适应滤波器。这个问题与传统的自适应噪声消除有所不同，主要体现在以下几点：所需的响应信号无法直接测量；只有衰减后的信号可用；主动噪声控制系统在自适应过程中必须考虑次级扬声器到麦克风的误差路径。

有关主动噪声控制任务的更多实施细节，请参见 S.M. Kuo 和 D.R. Morgan 的《主动噪声控制系统：算法与DSP实现》，Wiley-Interscience，纽约，1996。

次级传播路径

次级传播路径是反噪声从输出扬声器到安静区域内误差麦克风的路径。以下命令生成一个频带限制在160 - 2000 Hz范围内、过滤长度为0.1秒的扬声器到误差麦克风的脉冲响应。对于这个主动噪声控制任务，我们将使用8000 Hz的采样频率。

Fs     = 8e3;  % 8 kHz
N      = 800;  % 800 samples@8 kHz = 0.1 seconds
Flow   = 160;  % Lower band-edge: 160 Hz
Fhigh  = 2000; % Upper band-edge: 2000 Hz
delayS = 7;
Ast    = 20;   % 20 dB stopband attenuation
Nfilt  = 8;    % Filter order

% Design bandpass filter to generate bandlimited impulse response
filtSpecs = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
bandpass = design(filtSpecs,'cheby2','FilterStructure','df2tsos', ...
    'SystemObject',true);

% Filter noise to generate impulse response
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,'b');
xlabel('Time [sec]');
ylabel('Coefficient value');
title('True Secondary Path Impulse Response');

估计次级传播路径

主动噪声控制的第一步是估计次级传播路径的脉冲响应。此步骤通常在噪声控制之前进行，使用合成随机信号通过输出扬声器播放，同时不允许出现不需要的噪声。以下命令生成3.75秒的随机噪声以及在误差麦克风处测得的信号。

ntrS = 30000;
randomSignal = randn(ntrS,1); % Synthetic random signal to be played
secondaryPathGenerator = dsp.FIRFilter('Numerator',secondaryPathCoeffsActual.');
secondaryPathMeasured = secondaryPathGenerator(randomSignal) + ... % random signal propagated through secondary path
    0.01*randn(ntrS,1); % measurement noise at the microphone

设计次级传播路径估计

通常，次级路径滤波器估计的长度不需要像实际次级路径那样长，在大多数情况下，较短的长度也足以实现有效控制。我们将使用250个抽头的次级路径滤波器长度，对应于31毫秒的脉冲响应长度。虽然可以使用任何自适应FIR滤波算法来实现这一目的，但通常使用归一化LMS算法，因为它简单且稳健。输出信号和误差信号的图示显示，该算法在约10000次迭代后收敛。

M = 250;
muS = 0.1;
secondaryPathEstimator = dsp.LMSFilter('Method','Normalized LMS','StepSize', muS, ...
    'Length', M);
[yS,eS,SecondaryPathCoeffsEst] = secondaryPathEstimator(randomSignal,secondaryPathMeasured);

n = 1:ntrS;
figure, plot(n,secondaryPathMeasured,n,yS,n,eS);
xlabel('Number of iterations');
ylabel('Signal value');
title('Secondary Identification Using the NLMS Adaptive Filter');
legend('Desired Signal','Output Signal','Error Signal');

次级路径估计的准确性

次级路径脉冲响应估计的准确性如何？该图显示了真实路径和估计路径的系数。只有真实脉冲响应的尾部没有被准确估计。这一残余误差在所选任务中并不会显著影响主动噪声控制系统的性能。

figure, plot(t,secondaryPathCoeffsActual, ...
    t(1:M),SecondaryPathCoeffsEst, ...
    t,[secondaryPathCoeffsActual(1:M)-SecondaryPathCoeffsEst(1:M); secondaryPathCoeffsActual(M+1:N)]);
xlabel('Time [sec]');
ylabel('Coefficient value');
title('Secondary Path Impulse Response Estimation');
legend('True','Estimated','Error');

主要传播路径

要消除的噪声传播路径也可以通过线性滤波器来表征。以下命令生成一个输入到误差麦克风的脉冲响应，该脉冲响应带限于200 - 800 Hz范围，滤波器长度为0.1秒。

delayW = 15;
Flow   = 200; % Lower band-edge: 200 Hz
Fhigh  = 800; % Upper band-edge: 800 Hz
Ast    = 20;  % 20 dB stopband attenuation
Nfilt  = 10;  % Filter order

% Design bandpass filter to generate bandlimited impulse response
filtSpecs2 = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
bandpass2 = design(filtSpecs2,'cheby2','FilterStructure','df2tsos', ...
    'SystemObject',true);

% Filter noise to generate impulse response
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,'b');
xlabel('Time [sec]');
ylabel('Coefficient value');
title('Primary Path Impulse Response');

要消除的噪声

典型的主动噪声控制应用涉及旋转机械的声音，因为它们具有恼人的特性。在这里，我们合成生成可能来自典型电动机的噪声。

主动噪声控制的初始化

最流行的主动噪声控制自适应算法是滤波-X LMS 算法。该算法使用次级路径估计来计算一个输出信号，使其在误差传感器处的贡献与不希望的噪声发生相消干扰。参考信号是测量在噪声源附近的不希望声音的带噪声版本。我们将使用约 44 毫秒的控制器滤波器长度和 0.0001 的步长来处理这些信号统计。

% FIR Filter to be used to model primary propagation path
primaryPathGenerator = dsp.FIRFilter('Numerator',primaryPathCoeffs.');

% Filtered-X LMS adaptive filter to control the noise
L = 350;
muW = 0.0001;
noiseController = dsp.FilteredXLMSFilter('Length',L,'StepSize',muW, ...
    'SecondaryPathCoefficients',SecondaryPathCoeffsEst);

% Sine wave generator to synthetically create the noise
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); % Random initial phase
sine = audioOscillator('NumTones', La, 'Amplitude',A,'Frequency',F, ...
    'PhaseOffset',phase,'SamplesPerFrame',512,'SampleRate',Fs);

% Audio player to play noise before and after cancellation
player = audioDeviceWriter('SampleRate',Fs);

% Spectrum analyzer to show original and attenuated noise
scope = spectrumAnalyzer('SampleRate',Fs, ...
    'PlotAsTwoSidedSpectrum',false, ...
    'ShowLegend',true, ...
    'ChannelNames', {'Original noisy signal', 'Attenuated noise'});

使用滤波-X LMS 算法的主动噪声控制仿真

在这里，我们模拟主动噪声控制系统。为了强调差异，我们在前 200 次迭代中运行没有主动噪声控制的系统。在消除之前，听到的误差麦克风的声音具有这种电动机特有的工业“鸣叫”特征。

一旦启用自适应滤波器，结果算法在约 5 秒的适应时间后收敛。比较残余误差信号的频谱和原始噪声信号的频谱，我们发现大多数周期性成分已被显著衰减。然而，稳态消除性能可能在所有频率上并不均匀。这在应用于主动噪声控制任务的实际系统中往往是这样的。听到的误差信号中，恼人的“鸣叫”显著减少。

for m = 1:400
    % Generate synthetic noise by adding sine waves with random phase
    x = sine();
    d = primaryPathGenerator(x) + ...  % Propagate noise through primary path
        0.1*randn(size(x)); % Add measurement noise
    if m <= 200
        % No noise control for first 200 iterations
        e = d;
    else
        % Enable active noise control after 200 iterations
        xhat = x + 0.1*randn(size(x));
        [y,e] = noiseController(xhat,d);
    end
    player(e);     % Play noise signal
    scope([d,e]); % Show spectrum of original (Channel 1)
                     % and attenuated noise (Channel 2)
end
release(player); % Release audio device
release(scope); % Release spectrum analyzer