В этом примере показано, как применять адаптивные фильтры к ослаблению акустического шума посредством активного управления шумом.
При активном управлении шумом предпринимается попытка уменьшить объем нежелательного шума, распространяющегося по воздуху, с использованием электроакустической системы с использованием измерительных датчиков, таких как микрофоны и выходные приводы, такие как громкоговорители. Сигнал шума обычно поступает от какого-либо устройства, такого как вращающаяся машина, так что можно измерить шум вблизи своего источника. Целью системы управления активным шумом является создание «противошумного», который ослабляет нежелательный шум в желаемой тихой области с использованием адаптивного фильтра. Эта проблема отличается от традиционного адаптивного шумоподавления тем, что: - Требуемый ответный сигнал не может быть непосредственно измерен; доступен только ослабленный сигнал. - Активная система шуморегулирования должна учитывать в своей адаптации вторичный тракт ошибки «громкоговоритель-микрофон».
Для получения более подробной информации о задачах активного управления шумом см. S.M. Kuo и D.R. Morgan, «Active Noise Control Systems: Algorithms and DSP Implementations», Wiley-Interscience, Нью-Йорк, 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 мс. В то время как любой адаптивный алгоритм фильтрации 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 Гц и имеет длину фильтра 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');
Типичные активные приложения управления шумом включают звуки вращающегося оборудования из-за их раздражающих характеристик. Здесь мы синтетически генерируем шум, который может исходить от типичного электродвигателя.
Наиболее популярным адаптивным алгоритмом для активного управления шумом является алгоритм filtered-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 = dsp.SpectrumAnalyzer('SampleRate',Fs,'OverlapPercent',80, ... 'SpectralAverages',20,'PlotAsTwoSidedSpectrum',false, ... 'ShowLegend',true, ... 'ChannelNames', {'Original noisy signal', 'Attenuated noise'});
Здесь мы моделируем активную систему управления шумом. Чтобы подчеркнуть разницу, мы запускаем систему без активного управления шумом для первых 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
