Используя наименьший среднеквадратичный (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();
Для выполнения адаптации фильтру требуется два сигнала:
Опорный сигнал
Шумный сигнал, который содержит как требуемый сигнал, так и добавленную шумовую составляющую
Создание шумового сигнала с авторегрессионным шумом (определяется как 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);
Для сравнения вычислите оптимальный фильтр FIR Винера.
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
Наконец, сравните коэффициенты фильтра Винера с коэффициентами адаптивных фильтров. При адаптации адаптивные фильтры пытаются сходиться к коэффициентам Винера.
[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')
На этом графике показаны расчетные обучающие кривые для адаптивных фильтров 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')
