Используя наименьший средний квадрат (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
Наконец, сравните коэффициенты фильтра Винера с коэффициентами адаптивных фильтров. Во время адаптации адаптивные фильтры пытаются сходиться к коэффициентам Винера.
[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')
