Активное шумоподавление с помощью адаптивного КИХ-фильтра с дополнительной фильтрацией выходного сигнала

Данный пример демонстрирует применение адаптивных фильтров для ослабления акустического шума в системах активного шумоподавления.

Оглавление:

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

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

Более детально задачи активного шумоподавления рассмотрены в книге S.M. Kuo и D.R. Morgan, "Active Noise Control Systems: Algorithms and DSP Implementations", Wiley- Interscience, New York, 1996.

Путь вторичного распространения

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

Fs     = 8e3; % 8 КГц
N      = 800; % 800 отсчетов на 8 КГц = 0.1 секунды
Flow   = 160; % нижняя частота среза: 160 Гц
Fhigh  = 2000; % верхняя частота среза: 2000 Гц
delayS = 7;
Ast    = 20; % подавление 20 дБ
Nfilt  = 8; % порядок фильтра

% Создание полосового фильтра для имитации канала с ограниченной полосой
% пропускания
Fd = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
Hd = design(Fd,'cheby2','FilterStructure','df2tsos',...
    'SystemObject',true);

% Фильтрация шума для получения импульсной характеристики канала
H = step(Hd,[zeros(delayS,1); log(0.99*rand(N-delayS,1)+0.01).* ...
    sign(randn(N-delayS,1)).*exp(-0.01*(1:N-delayS)')]);
H = H/norm(H);

t = (1:N)/Fs;
plot(t,H,'b');
xlabel('Время, с');
ylabel('Значения коэффициентов');
title('Импульсная характеристика вторичного пути распространения сигнала');

Определение вторичного пути распространения

Первой задачей системы активного шумоподавления является определение импульсной характеристики пути вторичного распространения. Этот шаг обычно выполняется перед шумоподавлением с помощью синтезированного случайного сигнала, проигрываемого динамиками, при отсутствии шума. Нижеприведенные команды генерируют случайный сигнал длительностью 3.75 с, а также измеренный микрофоном сигнал с ошибкой.

ntrS = 30000;
s = randn(ntrS,1); % синтез случайного сигнал
Hfir = dsp.FIRFilter('Numerator',H.');
dS = step(Hfir,s) + ... % случайный сигнал прошедший через вторичный канал
    0.01*randn(ntrS,1); % шум микрофона

Создание фильтра для оценки вторичного пути распространения

В большинстве случаев для адекватного управления алгоритмом длительность отклика фильтра, оценивающего вторичный путь распространения, должна быть короче самого вторичного пути. Мы будем использовать фильтр 250 порядка, что соответствует импульсной характеристике длиной 31 мс. Для этой цели подходит любой алгоритм адаптивной КИХ-фильтрации, но обычно используют нормализованный алгоритм нахождения минимальной среднеквадратической ошибки (normalized LMS-алгоритм) ввиду его простоты и устойчивости.

M = 250;
muS = 0.1;
hNLMS = dsp.LMSFilter('Method','Normalized LMS','StepSize', muS,...
    'Length', M);
[yS,eS,Hhat] = step(hNLMS,s,dS);

n = 1:ntrS;
plot(n,dS,n,yS,n,eS);
xlabel('Число итераций');
ylabel('Уровень сигнала');
title('Идентификация вторичного пути распространения с NLMS-алгоритма');
legend('Ожидаемый сигнал','Сигнал на выходе','Сигнал ошибки');

Точность полученной оценки

Как точно оценивается импульсная характеристика вторичного пути? Этот график показывает коэффициенты настоящего пути и пути, рассчитанного алгоритмом. Только конец полученной импульсной характеристики имеет неточности. Эта остаточная ошибка не навредит производительности системы активного шумоподавления во время ее работы над выбранной задачей.

plot(t,H,t(1:M),Hhat,t,[H(1:M)-Hhat(1:M); H(M+1:N)]);
xlabel('Время, с');
ylabel('Значения коэффициентов');
title('Определение импульсной характеристики вторичного пути распространения');
legend('Действительная','Оцененная','Ошибка');

Основной путь распространения сигнала

Путь распространения шума, который должен быть подавлен, может быть также описан с помощью линейного фильтра. Следующие команды генерируют импульсную характеристику пути источник шума-микрофон с ограниченной полосой 200-800 Гц и имеет длительность отклика, равную 0.1 с.

delayW = 15;
Flow   = 200; % нижняя частота среза: 200 Hz
Fhigh  = 800; % верхняя частота среза: 800 Hz
Ast    = 20;  % подавление 20 дБ
Nfilt  = 10;  % порядок фильтра

% Создание полосового фильтра для имитации импульсного отклика с
% ограниченной полосой
Fd2 = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
Hd2 = design(Fd2,'cheby2','FilterStructure','df2tsos',...
    'SystemObject',true);

% Фильтрация шума для получения импульсной характеристики
G = step(Hd2,[zeros(delayW,1); log(0.99*rand(N-delayW,1)+0.01).*...
    sign(randn(N-delayW,1)).*exp(-0.01*(1:N-delayW)')]);
G = G/norm(G);

plot(t,G,'b');
xlabel('Время, с');
ylabel('Значения коэффициентов');
title('Импульсная характеристика первичного пути распространения');

Подавляемый шум

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

Инициализация системы

Самым распространенным алгоритмом для систем активного шумоподавления является LMS- алгоритм с дополнительной фильтрацией выходного сигнала фильтра перед формированием сигнала ошибки (Filtered-x LMS algorithm). Этот алгоритм использует оценку вторичного пути распространения для расчета выходного сигнала, который разрушительно влияет на нежелательный шум в области датчика измерения ошибки. Опорным сигналом является зашумленная версия нежелательного звука, измеренная вблизи его источника. Мы будем использовать управляемый фильтр с длительностью отклика около 44 мс и шагом подстройки, равным 0.0001.

% КИХ фильтр используемый для моделирования первичного пути распространения
Hfir = dsp.FIRFilter('Numerator',G.');

% Адаптивный фильтр реализующий алгоритм Filtered-X LMS
L = 350;
muW = 0.0001;
Hfx = dsp.FilteredXLMSFilter('Length',L,'StepSize',muW,...
    'SecondaryPathCoefficients',Hhat);

% Синтез шума с помощью синусоид
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); % случайная начальная фаза
Hsin = dsp.SineWave('Amplitude',A,'Frequency',F,'PhaseOffset',phase,...
    'SamplesPerFrame',512,'SampleRate',Fs);

% Проигрыватель аудио для воспроизведения результатов работы алгоритма
Hpa = dsp.AudioPlayer('SampleRate',Fs,'QueueDuration',2);

% Анализотор спектра
Hsa = dsp.SpectrumAnalyzer('SampleRate',Fs,'OverlapPercent',80,...
    'SpectralAverages',20,'PlotAsTwoSidedSpectrum',false,...
    'ShowLegend',true);

Симуляция разработанной системы активного шумоподавления

Здесь мы сымитируем работу системы активного шумоподавления. Чтобы подчеркнуть разницу, первые 200 итераций шумоподавление будет отключено. Звук на микрофоне до подавления представляет характерный «вой» промышленных моторов.

Результирующий алгоритм сходится примерно через 5 с (имитационных) после включения адаптивного фильтра. Сравнивая спектры сигнала остаточной ошибки и исходного зашумленного сигнала, можно наблюдать, что большая часть периодичных компонент была успешно подавлена. Однако эффективность стационарного шумоподавления может быть неравномерна по всем частотам. Такое часто бывает в реальных системах, применяемых для задач активной борьбы с шумом. При прослушивании сигнала ошибки раздражающий «вой» значительно снижается.

for m = 1:400
    s = step(Hsin); % генерация синусоид со случайной фазой
    x = sum(s,2);   % генерация шума сложением всех синусоид
    d = step(Hfir,x) + ... % распространение шума через первичный канал
        0.1*randn(size(x)); % добавление шума, сопроводающего процесс измерения
    if m <= 200
        % отключение шумоподавления на первые 200 итераций
        e = d;
    else
        % включение алгоритма шумоподавления
        xhat = x + 0.1*randn(size(x));
        [y,e] = step(Hfx,xhat,d);
    end
    step(Hpa,e);     % воспроизведение сигнала на выходе
    step(Hsa,[d,e]); % спектр исходного (канал 1) и ослабленного (канал 2) сигналов
end
release(Hpa); % отключение динамиков
release(Hsa); % отключение спектроанализатора
Warning: The queue has underrun by 3456 samples. Try increasing queue duration,
buffer size, or throughput rate.