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

Сначала мы загрузим сигнал интереса, короткий клип из хора Генделя «Аллилуйя».

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

Теперь удаление чистой синусоиды из синусоиды плюс музыкальный сигнал не является особенно сложным, если известна частота нарушающей синусоиды. Простая двухполюсная система с двумя нулевыми узкополосными фильтрами может выполнить эту задачу. Итак, давайте сделаем задачу немного сложнее, добавив модулированный 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

«Варбл» в сигнале хорошо слышен. Узкополосный фильтр с фиксированным коэффициентом не удаляет синусоиду с FM-модуляцией. Посмотрим, сможет ли 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);