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

Открыть Live Script

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

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

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

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();

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

  • A опорного сигнала

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

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

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

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

The 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

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

Постройте график полученной обесцененной синусоиды для каждого фильтра - фильтра Винера, адаптивного фильтра 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. The axes 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. The axes 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. The axes contains 4 objects of type line. These objects represent MMSE, EMSE, Predicted LMS learning curve, LMS learning curve.

Документация DSP System Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте