Загрузка интересующего сигнала

Сначала мы загрузим интересующий сигнал, короткий клип из хора Handel's Hallelujah.

audioReader = dsp.AudioFileReader('handel.ogg','SamplesPerFrame',44100);
timeScope = timescope('SampleRate',audioReader.SampleRate,...
    'YLimits',[-1,1],'TimeSpanSource','property','TimeSpan',1);

while ~isDone(audioReader)
    x = audioReader() / 2;
    timeScope(x);
end

Прослушивание звукового клипа

Интересующий сигнал можно прослушать с помощью устройства записи аудиоустройств.

release(audioReader);
audioWriter = audioDeviceWriter;
while ~isDone(audioReader)
    x = audioReader() / 2;
    audioWriter(x);
end

Формирование шумового сигнала

Давайте сделаем сигнал периодического шума - синусоиду с частотой 1000 Гц.

sine = dsp.SineWave('Amplitude',0.5,'Frequency',1000,...
    'SampleRate',audioReader.SampleRate,...
    'SamplesPerFrame',audioReader.SamplesPerFrame);

Теперь нарисуем 10 мс этой синусоиды выше. Он показывает 10 периодов в 10 мс, как и должно.

eta = sine();
Fs  = sine.SampleRate;
plot(1/Fs:1/Fs:0.01,eta(1:floor(0.01*Fs)));
xlabel('Time [sec]');
ylabel('Amplitude');
title('Noise Signal, eta(n)');

Прослушивание шума

Периодический шум - это чистый тон. Следующий код воспроизводит одну секунду (один кадр из 44100 выборок) шумового сигнала.

eta = sine();
release(audioWriter);
audioWriter(eta);

Измеренный сигнал

Сигнал, который мы на самом деле измеряем, является суммой этих двух сигналов, и мы называем этот сигнал s (n). График s (n) показывает, что огибающая музыкального сигнала в значительной степени скрыта. Слушая 3-секундный клип от измеренного сигнала, шум явно заметен...

release(audioReader);
release(timeScope);
release(audioWriter);
count = 1;
while count < 4
    s = (audioReader() / 2) + sine();
    timeScope(s);
    audioWriter(s);
    count = count + 1;
end

Адаптивная конфигурация фильтра

Адаптивный усилитель линии (ALE) основан на простой концепции линейного прогнозирования. Почти периодический сигнал можно точно предсказать, используя линейные комбинации его прошлых выборок, тогда как непериодический сигнал не может. Таким образом, задержанная версия измеренного сигнала s (n-D) используется в качестве опорного входного сигнала x (n) для адаптивного фильтра, и требуемый ответный сигнал d (n) делается равным s (n). Параметрами для выбора в такой системе являются задержка D сигнала и длина L фильтра, используемые в адаптивной линейной оценке. Величина задержки зависит от величины корреляции в интересующем сигнале. Так как у нас нет этого сигнала (если бы мы это сделали, нам бы не нужен ALE!), мы просто выберем значение D = 100 и изменим его позже. Такой выбор предполагает, что образцы Хора Аллилуйи некоррелированы, если они находятся на расстоянии более 12 мс друг от друга. Кроме того, мы выберем значение L = 32 для адаптивного фильтра, хотя это тоже может быть изменено.

D = 100;
delay = dsp.Delay(D);

Наконец, мы будем использовать некоторые блочные адаптивные алгоритмы, которые требуют, чтобы длины векторов для x (n) и d (n) были целыми кратными длине блока. Мы выберем длину блока N = 49, с которой начать.

Блок LMS

Первым алгоритмом, который мы рассмотрим, является алгоритм Block LMS. Этот алгоритм аналогичен хорошо известному алгоритму наименьшего квадрата (LMS), за исключением того, что он использует обновления блочных коэффициентов вместо обновлений коэффициентов выборки за выборкой. Алгоритм блочной LMS нуждается в длине фильтра, длине N блока и значении mu размера шага. Как нам выбрать мю? Начнем со значения mu = 0,0001 и уточним его в ближайшее время.

L  = 32;
N  = 49;
mu = 0.0001;
blockLMSFilter = ...
    dsp.BlockLMSFilter('Length',L,'StepSize',mu,'BlockSize',N);

Запуск фильтра

Выходной сигнал y (n) должен в значительной степени содержать периодическую синусоиду, тогда как сигнал ошибки e (n) должен содержать музыкальную информацию, если мы сделали все правильно. Поскольку у нас есть исходный музыкальный сигнал v (n), мы можем построить график e (n) против v (n) на том же графике, который показан выше, вместе с остаточным сигналом e (n) -v (n). Похоже, что система сходится примерно через 5 секунд адаптации с таким размером шага. Реальное доказательство, однако, получается путём прослушивания;

release(audioReader);
release(timeScope);
release(audioWriter);
while ~isDone(audioReader)
    x = audioReader() / 2;
    s = x + sine();
    d = delay(s);
    [y,e] = blockLMSFilter(s,d);
    timeScope(e);
    audioWriter(e);
end

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

Источник шума FM

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

eta = 0.5 * sin(2*pi*1000/Fs*(0:396899)' + 10*sin(2*pi/Fs*(0:396899)'));
signalSource = dsp.SignalSource(eta,...
    'SamplesPerFrame',audioReader.SamplesPerFrame,...
    'SignalEndAction','Cyclic repetition');
release(audioReader);
release(timeScope);
release(audioWriter);
while ~isDone(audioReader)
    x = audioReader() / 2;
    s = x + signalSource();
    timeScope(s);
    audioWriter(s);
end

«Бородавка» в сигнале отчетливо слышна. Режекторный фильтр с фиксированным коэффициентом не удаляет ЧМ-модулированную синусоиду. Посмотрим, сможет ли ALE на основе Block LMS. Мы увеличим значение размера шага до mu = 0,005, чтобы помочь ALE отслеживать изменения шумового сигнала.

mu = 0.005;
release(blockLMSFilter);
blockLMSFilter.StepSize = mu;

Запуск адаптивного фильтра

Теперь мы фильтруем шумный музыкальный сигнал адаптивным фильтром и сравниваем ошибку с шумным музыкальным сигналом.

release(audioReader);
release(timeScope);
release(audioWriter);
while ~isDone(audioReader)
    x = audioReader() / 2;
    s = x + signalSource();
    d = delay(s);
    [y,e] = blockLMSFilter(s,d);
    timeScope([x,e]);
    audioWriter(e);
end

release(audioReader);
release(timeScope);
release(audioWriter);