Активное управление шумом

При активном регулировании шума пытаются уменьшить объем нежелательного шума, распространяющегося по воздуху, с помощью электроакустической системы, использующей измерение датчики, такие как микрофоны и выходные приводы, такие как громкоговорители. Сигнал шума обычно поступает от некоторого устройства, такого как вращающаяся машина, так что можно измерить шум около его источника. Цель активной системы управления шумом состоит в том, чтобы создать «антишум», который ослабляет нежелательный шум в желаемой тихой области, используя адаптивный фильтр. Эта задача отличается от традиционного адаптивного шумоподавления тем, что: - Нужный сигнал отклика не может быть непосредственно измерен; доступен только ослабленный сигнал. - Активная система управления шумом должна учитывать вторичный путь ошибки «громкоговоритель - микрофон» при его адаптации.

Для получения дополнительной детали реализации задач активного контроля шума см. С.М. Куо и D.R. Morgan, «Active Noise Control Systems: Algorithms and DSP Implementations», Wiley-Interscience, New York, 1996.

Вторичный путь распространения

Вторичный путь распространения является путем, которым противошумовые взятия от выхода громкоговорителя до микрофона ошибки в тихой зоне. Следующие команды генерируют импульсную характеристику микрофона от громкоговорителя до ошибки, которая ограничена диапазоном 160-2000 Гц и длиной фильтра 0,1 секунды. Для этой задачи активного регулирования шума мы будем использовать частоту дискретизации 8000 Гц.

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 мс. Хотя для этой цели могли бы использоваться любые адаптивные конечные импульсные характеристики алгоритма фильтрации, нормализованный 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 Гц и имеет длину фильтра 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');

Шум, который будет отменен

Типичные активные приложения управления шумом включают звуки вращающегося машинного оборудования из-за их раздражающих характеристик. Здесь мы синтетически генерируем шум, который может исходить от типичного электродвигателя.

Инициализация активного контроля шума

Самым популярным адаптивным алгоритмом для активного контроля шума является LMS-алгоритм filtered-X. Этот алгоритм использует оценку вторичного пути, чтобы вычислить сигнал выхода, чей вклад в датчик ошибки разрушительно мешает нежелательному шуму. Опорный сигнал является шумной версией нежелательного звука, измеренного около его источника. Мы будем использовать фильтр контроллера длину около 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 = dsp.SpectrumAnalyzer('SampleRate',Fs,'OverlapPercent',80, ...
    'SpectralAverages',20,'PlotAsTwoSidedSpectrum',false, ...
    'ShowLegend',true, ...
    'ChannelNames', {'Original noisy signal', 'Attenuated noise'});

Симуляция активного контроля шума с использованием фильтрованного LMS-алгоритма X

Здесь мы моделируем активную систему контроля шума. Чтобы подчеркнуть различие, мы запускаем систему без активного контроля шума для первых 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