Этот пример показывает, как применить адаптивные фильтры к ослаблению акустического шума через активное управление шумом.
При активном регулировании шума пытаются уменьшить объем нежелательного шума, распространяющегося по воздуху, с помощью электроакустической системы, использующей измерение датчики, такие как микрофоны и выходные приводы, такие как громкоговорители. Сигнал шума обычно поступает от некоторого устройства, такого как вращающаяся машина, так что можно измерить шум около его источника. Цель активной системы управления шумом состоит в том, чтобы создать «антишум», который ослабляет нежелательный шум в желаемой тихой области, используя адаптивный фильтр. Эта задача отличается от традиционного адаптивного шумоподавления тем, что: - Нужный сигнал отклика не может быть непосредственно измерен; доступен только ослабленный сигнал. - Активная система управления шумом должна учитывать вторичный путь ошибки «громкоговоритель - микрофон» при его адаптации.
Для получения дополнительной детали реализации задач активного контроля шума см. С.М. Куо и 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'});
Здесь мы моделируем активную систему контроля шума. Чтобы подчеркнуть различие, мы запускаем систему без активного контроля шума для первых 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
