Улучшение сигнала Используя LMS и алгоритмы NLMS

Используя наименьшее количество среднего квадратичного (LMS) и нормированные LMS-алгоритмы, извлеките желаемый сигнал из поврежденного шумом сигнала путем отфильтровывания шума. Оба этих алгоритма доступны с dsp.LMSFilter Система object™.

Создайте сигналы для адаптации

Желаемый сигнал (выход от процесса) является синусоидой с 1 000 выборок на систему координат.

sine = dsp.SineWave('Frequency',375,'SampleRate',8000,'SamplesPerFrame',1000)
sine = 
  dsp.SineWave with properties:

          Amplitude: 1
          Frequency: 375
        PhaseOffset: 0
      ComplexOutput: false
             Method: 'Trigonometric function'
         SampleRate: 8000
    SamplesPerFrame: 1000
     OutputDataType: 'double'

s = sine();

Чтобы выполнить адаптацию, фильтр требует двух сигналов:

  • Опорный сигнал

  • Сигнал с шумом, который содержит и желаемый сигнал и добавленный шумовой компонент

Сгенерируйте шумовой сигнал

Создайте шумовой сигнал с авторегрессивным шумом (заданный как v1). В авторегрессивном шуме шум во время t зависит только от предыдущих значений и случайного воздействия.

v = 0.8*randn(sine.SamplesPerFrame,1); % Random noise part.
ar = [1,1/2];          % Autoregression coefficients.
ARfilt = dsp.IIRFilter('Numerator',1,'Denominator',ar)
ARfilt = 
  dsp.IIRFilter with properties:

            Structure: 'Direct form II transposed'
            Numerator: 1
          Denominator: [1 0.5000]
    InitialConditions: 0

  Show all properties

v1 = ARfilt(v);

Повредите желаемый сигнал создать сигнал с шумом

Чтобы сгенерировать сигнал с шумом, который содержит и желаемый сигнал и шум, добавьте шумовой v1 сигнала к желаемому s сигнала. Поврежденная шумом синусоида x :

x = s + v1;

Адаптивная обработка фильтра стремится восстановить s от x путем удаления v1. Завершать сигналы должно было выполнить адаптивную фильтрацию, процесс адаптации требует опорного сигнала.

Создайте опорный сигнал

Задайте сигнал скользящего среднего значения v2 это коррелируется с v1. v2 сигнала опорный сигнал для этого примера.

ma = [1, -0.8, 0.4, -0.2];
MAfilt = dsp.FIRFilter('Numerator',ma)
MAfilt = 
  dsp.FIRFilter with properties:

            Structure: 'Direct form'
      NumeratorSource: 'Property'
            Numerator: [1 -0.8000 0.4000 -0.2000]
    InitialConditions: 0

  Show all properties

v2 = MAfilt(v);

Создайте два адаптивных фильтра

Два подобных, шестой порядок адаптивные фильтры — LMS и NLMS — формируют базис этого примера. Установите порядок как переменную в MATLAB™ и создайте фильтры.

L = 7;
lms = dsp.LMSFilter(L,'Method','LMS')
lms = 
  dsp.LMSFilter with properties:

                   Method: 'LMS'
                   Length: 7
           StepSizeSource: 'Property'
                 StepSize: 0.1000
            LeakageFactor: 1
        InitialConditions: 0
           AdaptInputPort: false
    WeightsResetInputPort: false
            WeightsOutput: 'Last'

  Show all properties

nlms = dsp.LMSFilter(L,'Method','Normalized LMS')
nlms = 
  dsp.LMSFilter with properties:

                   Method: 'Normalized LMS'
                   Length: 7
           StepSizeSource: 'Property'
                 StepSize: 0.1000
            LeakageFactor: 1
        InitialConditions: 0
           AdaptInputPort: false
    WeightsResetInputPort: false
            WeightsOutput: 'Last'

  Show all properties

Выберите размер шага

Подобный LMS-алгоритмам имеют размер шага, который определяет объем коррекции, примененной, когда фильтр адаптируется от одной итерации до следующего. Размер шага, который является слишком маленькими увеличениями время для фильтра, чтобы сходиться на наборе коэффициентов. Размер шага, который является слишком большой силой, заставляет адаптирующийся фильтр отличать и никогда не достигать сходимости. В этом случае получившийся фильтр не может быть устойчивым.

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

maxstep функция dsp.LMSFilter объект определяет максимальный размер шага, подходящий для каждого адаптивного алгоритма фильтра LMS, который гарантирует, что фильтр сходится к решению. Часто, обозначение для размера шага является µ.

[mumaxlms,mumaxmselms] = maxstep(lms,x)
mumaxlms = 0.2127
mumaxmselms = 0.1312
[mumaxnlms,mumaxmsenlms] = maxstep(nlms,x)
mumaxnlms = 2
mumaxmsenlms = 2

Установите адаптирующийся размер шага фильтра

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

lms.StepSize  = mumaxmselms/30 
lms = 
  dsp.LMSFilter with properties:

                   Method: 'LMS'
                   Length: 7
           StepSizeSource: 'Property'
                 StepSize: 0.0044
            LeakageFactor: 1
        InitialConditions: 0
           AdaptInputPort: false
    WeightsResetInputPort: false
            WeightsOutput: 'Last'

  Show all properties

nlms.StepSize = mumaxmsenlms/20 
nlms = 
  dsp.LMSFilter with properties:

                   Method: 'Normalized LMS'
                   Length: 7
           StepSizeSource: 'Property'
                 StepSize: 0.1000
            LeakageFactor: 1
        InitialConditions: 0
           AdaptInputPort: false
    WeightsResetInputPort: false
            WeightsOutput: 'Last'

  Show all properties

Отфильтруйте с адаптивными фильтрами

Вы настроили параметры адаптивных фильтров и теперь готовы отфильтровать сигнал с шумом. Опорный сигнал v2 является входом к адаптивным фильтрам. x является желаемым сигналом в этой настройке.

Посредством адаптации y, выхода фильтров, пытается эмулировать x максимально тесно.

Поскольку v2 коррелируется только с шумовым компонентом v1 x, это может только действительно эмулировать v1. Сигнал ошибки (желаемый x), минус фактический выход y, составляет оценку части x, который не коррелируется с v2 — s, сигнал извлечь из x.

[~,elms,wlms] = lms(v2,x);
[~,enlms,wnlms] = nlms(v2,x);

Вычислите оптимальное решение

Для сравнения вычислите оптимального КИХ Винеровский фильтр.

reset(MAfilt);
bw = firwiener(L-1,v2,x); % Optimal FIR Wiener filter
MAfilt = dsp.FIRFilter('Numerator',bw)
MAfilt = 
  dsp.FIRFilter with properties:

            Structure: 'Direct form'
      NumeratorSource: 'Property'
            Numerator: [1.0001 0.3060 0.1050 0.0482 0.1360 0.0959 0.0477]
    InitialConditions: 0

  Show all properties

yw = MAfilt(v2); % Estimate of x using Wiener filter
ew = x - yw; % Estimate of actual sinusoid

Постройте график результатов

Постройте получившуюся denoised синусоиду для каждого фильтра — Винеровского фильтра, адаптивного фильтра LMS, и адаптивного фильтра NLMS — чтобы сравнить эффективность различных методов.

n = (1:1000)';
plot(n(900:end),[ew(900:end), elms(900:end),enlms(900:end)])
legend('Wiener filter denoised sinusoid',...
    'LMS denoised sinusoid','NLMS denoised sinusoid')
xlabel('Time index (n)')
ylabel('Amplitude')

Как контрольная точка, включайте сигнал с шумом как пунктирную линию в графике.

hold on
plot(n(900:end),x(900:end),'k:')
xlabel('Time index (n)')
ylabel('Amplitude')
hold off

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent Wiener filter denoised sinusoid, LMS denoised sinusoid, NLMS denoised sinusoid.

Сравните итоговые коэффициенты

Наконец, сравните Винеровские коэффициенты фильтра с коэффициентами адаптивных фильтров. При адаптации адаптивные фильтры пытаются сходиться к Винеровским коэффициентам.

[bw.' wlms wnlms]
ans = 7×3

    1.0001    0.8644    0.9690
    0.3060    0.1198    0.2661
    0.1050   -0.0020    0.1226
    0.0482   -0.0046    0.1074
    0.1360    0.0680    0.2210
    0.0959    0.0214    0.1940
    0.0477    0.0292    0.1127

Сбросьте фильтр перед фильтрацией

Можно сбросить внутренние состояния фильтра в любое время путем вызова reset функция на объекте фильтра.

Например, эти последовательные вызовы производят тот же выход после сброса объекта.

[ylms,elms,wlms] = lms(v2,x);
[ynlms,enlms,wnlms] = nlms(v2,x);

Если вы не сбрасываете объект фильтра, фильтр использует конечные состояния и коэффициенты от предыдущего запуска как начальные условия и набор данных для следующего запуска.

Исследуйте сходимость через кривые обучения

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

Эта демонстрация использует 25 демонстрационной реализации шумных синусоид.

reset(ARfilt)
reset(sine);
release(sine);
n = (1:5000)';
sine.SamplesPerFrame = 5000
sine = 
  dsp.SineWave with properties:

          Amplitude: 1
          Frequency: 375
        PhaseOffset: 0
      ComplexOutput: false
             Method: 'Trigonometric function'
         SampleRate: 8000
    SamplesPerFrame: 5000
     OutputDataType: 'double'

s = sine();
nr = 25;
v = 0.8*randn(sine.SamplesPerFrame,nr);
ARfilt = dsp.IIRFilter('Numerator',1,'Denominator',ar)
ARfilt = 
  dsp.IIRFilter with properties:

            Structure: 'Direct form II transposed'
            Numerator: 1
          Denominator: [1 0.5000]
    InitialConditions: 0

  Show all properties

v1 = ARfilt(v);
x = repmat(s,1,nr) + v1;
reset(MAfilt);
MAfilt = dsp.FIRFilter('Numerator',ma)
MAfilt = 
  dsp.FIRFilter with properties:

            Structure: 'Direct form'
      NumeratorSource: 'Property'
            Numerator: [1 -0.8000 0.4000 -0.2000]
    InitialConditions: 0

  Show all properties

v2 = MAfilt(v);

Вычислите кривые обучения

Теперь вычислите среднеквадратическую ошибку. Чтобы ускорить вещи, вычислите ошибку каждые 10 выборок.

Во-первых, сбросьте адаптивные фильтры, чтобы избегать использования коэффициентов, которые это уже вычислило и состояния, которые это сохранило. Затем постройте кривые обучения для LMS и адаптивных фильтров NLMS.

reset(lms);
reset(nlms);
M = 10; % Decimation factor
mselms = msesim(lms,v2,x,M);
msenlms = msesim(nlms,v2,x,M);
plot(1:M:n(end),mselms,'b',1:M:n(end),msenlms,'g')
legend('LMS learning curve','NLMS learning curve')
xlabel('Time index (n)')
ylabel('MSE')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent LMS learning curve, NLMS learning curve.

В этом графике вы видите расчетные кривые обучения для LMS и адаптивных фильтров NLMS.

Вычислите теоретические кривые обучения

Для LMS и алгоритмов NLMS, функции в тулбоксе помогают вам вычислить теоретические кривые обучения, наряду с минимальной среднеквадратической ошибкой (MMSE), избыточной среднеквадратической ошибкой (EMSE) и средним значением коэффициентов.

MATLAB может занять время, чтобы вычислить кривые. Фигура, показанная после кода, строит предсказанные и фактические кривые LMS.

reset(lms);
[mmselms,emselms,meanwlms,pmselms] = msepred(lms,v2,x,M);
x = 1:M:n(end);
y1 = mmselms*ones(500,1);
y2 = emselms*ones(500,1);
y3 = pmselms;
y4 = mselms;
plot(x,y1,'m',x,y2,'b',x,y3,'k',x,y4,'g')
legend('MMSE','EMSE','Predicted LMS learning curve',...
    'LMS learning curve')
xlabel('Time index (n)')
ylabel('MSE')

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent MMSE, EMSE, Predicted LMS learning curve, LMS learning curve.