Мотивация

Сглаживание - это то, как мы обнаруживаем важные закономерности в наших данных, оставляя неважные вещи (т.е. шум). Для сглаживания используется фильтрация. Цель сглаживания состоит в медленном изменении стоимости, чтобы было легче видеть тенденции в наших данных.

Иногда при изучении входных данных может потребоваться сглаживание данных, чтобы увидеть тренд в сигнале. В нашем примере мы имеем набор показаний температуры по Цельсию, снятых каждый час в аэропорту Логан в Бостоне в течение всего января 2011 года.

load bostemp
days = (1:31*24)/24;
plot(days, tempC)
axis tight
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

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

Фильтр скользящего среднего

В простейшей форме фильтр скользящего среднего длины N принимает среднее из каждых N последовательных выборок формы сигнала.

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

hoursPerDay = 24;
coeff24hMA = ones(1, hoursPerDay)/hoursPerDay;

avg24hTempC = filter(coeff24hMA, 1, tempC);
plot(days,[tempC avg24hTempC])
legend('Hourly Temp','24 Hour Average (delayed)','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Задержка фильтра

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

Любой симметричный фильтр длиной N будет иметь задержку (N-1 )/2 выборок. Мы можем учесть эту задержку вручную.

fDelay = (length(coeff24hMA)-1)/2;
plot(days,tempC, ...
     days-fDelay/24,avg24hTempC)
axis tight
legend('Hourly Temp','24 Hour Average','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Извлечение средних разниц

Альтернативно, мы также можем использовать фильтр скользящего среднего, чтобы получить лучшую оценку того, как время суток влияет на общую температуру. Для этого сначала вычитают сглаженные данные из почасовых измерений температуры. Затем разделите разностные данные на дни и возьмите среднее значение за все 31 день месяца.

figure
deltaTempC = tempC - avg24hTempC;
deltaTempC = reshape(deltaTempC, 24, 31).';

plot(1:24, mean(deltaTempC))
axis tight
title('Mean temperature differential from 24 hour average')
xlabel('Hour of day (since midnight)')
ylabel('Temperature difference (\circC)')

Извлечение пиковой оболочки

Иногда нам также хотелось бы иметь плавно изменяющуюся оценку того, как ежедневно изменяются максимумы и минимумы нашего температурного сигнала. Для этого мы можем использовать envelope функция для подключения экстремальных максимумов и минимумов, обнаруженных за подмножество 24-часового периода. В этом примере мы гарантируем наличие по крайней мере 16 часов между каждым экстремальным максимумом и экстремальным минимумом. Мы также можем понять, как происходят максимумы и минимумы, взяв среднее значение между двумя крайностями.

[envHigh, envLow] = envelope(tempC,16,'peak');
envMean = (envHigh+envLow)/2;

plot(days,tempC, ...
     days,envHigh, ...
     days,envMean, ...
     days,envLow)
   
axis tight
legend('Hourly Temp','High','Mean','Low','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Фильтры взвешенного скользящего среднего

Другие типы фильтров скользящего среднего не взвешивают каждый образец одинаково.

Другой общий фильтр следует биномиальному расширению ${\left[\mathrm{}\mathrm{}\mathrm{}\mathrm{1/2,1/2}\right]}^{\mathit{}}$n. Этот тип фильтра аппроксимирует нормальную кривую для больших значений n. Он полезен для фильтрации высокочастотного шума для малых n. Чтобы найти коэффициенты для биномиального фильтра, сверните $\mathrm{[}\mathrm{}\mathrm{}\mathrm{}1/2,1/2$] с собой и затем итеративно сверните выходной сигнал с $\mathrm{[}\mathrm{}\mathrm{}\mathrm{}1/2,1/2$] заданное количество раз. В этом примере используйте пять итераций.

h = [1/2 1/2];
binomialCoeff = conv(h,h);
for n = 1:4
    binomialCoeff = conv(binomialCoeff,h);
end

figure
fDelay = (length(binomialCoeff)-1)/2;
binomialMA = filter(binomialCoeff, 1, tempC);
plot(days,tempC, ...
     days-fDelay/24,binomialMA)
axis tight
legend('Hourly Temp','Binomial Weighted Average','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

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

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

alpha = 0.45;
exponentialMA = filter(alpha, [1 alpha-1], tempC);
plot(days,tempC, ...
     days-fDelay/24,binomialMA, ...
     days-1/24,exponentialMA)

axis tight
legend('Hourly Temp', ...
       'Binomial Weighted Average', ...
       'Exponential Weighted Average','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Увеличьте изображение показаний за один день.

axis([3 4 -5 2])

Фильтры Савицкого-Голая

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

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

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

cubicMA   = sgolayfilt(tempC, 3, 7);
quarticMA = sgolayfilt(tempC, 4, 7);
quinticMA = sgolayfilt(tempC, 5, 9);
plot(days,[tempC cubicMA quarticMA quinticMA])
legend('Hourly Temp','Cubic-Weighted MA', 'Quartic-Weighted MA', ...
       'Quintic-Weighted MA','location','southeast')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')
axis([3 5 -5 2])

Передискретизация

Иногда полезно выполнить повторную выборку сигнала, чтобы правильно применить скользящее среднее.

В нашем следующем примере мы произвели выборку напряжения разомкнутого контура на входе аналогового прибора при наличии помех от шума линии переменного тока 60 Гц. Мы произвели выборку напряжения с частотой дискретизации 1 кГц.

load openloop60hertz
fs = 1000;
t = (0:numel(openLoopVoltage)-1) / fs;
plot(t,openLoopVoltage)
ylabel('Voltage (V)')
xlabel('Time (s)')
title('Open-loop Voltage Measurement')

Попытаемся устранить эффект линейного шума с помощью фильтра скользящего среднего.

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

Существует приблизительно 1000/60 = 16,667 выборок в полном цикле 60 Гц при выборке при 1000 Гц. Попробуем «округлить» и использовать 17-точечный фильтр. Это даст нам максимальную фильтрацию при основной частоте 1000 Гц/17 = 58,82 Гц.

plot(t,sgolayfilt(openLoopVoltage,1,17))
ylabel('Voltage (V)')
xlabel('Time (s)')
title('Open-loop Voltage Measurement')
legend('Moving average filter operating at 58.82 Hz', ...
       'Location','southeast')

Заметим, что хотя напряжение значительно сглажено, оно все еще содержит небольшую пульсацию 60 Гц.

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

Если мы проведем повторную выборку сигнала при 17 * 60 Гц = 1020 Гц, мы можем использовать наш 17-точечный фильтр скользящего среднего, чтобы удалить шум линии 60 Гц.

fsResamp = 1020;
vResamp = resample(openLoopVoltage, fsResamp, fs);
tResamp = (0:numel(vResamp)-1) / fsResamp;
vAvgResamp = sgolayfilt(vResamp,1,17);
plot(tResamp,vAvgResamp)
ylabel('Voltage (V)')
xlabel('Time (s)')
title('Open-loop Voltage Measurement')
legend('Moving average filter operating at 60 Hz', ...
    'Location','southeast')

Медианный фильтр

Скользящее среднее, взвешенное скользящее среднее и фильтры Савицки-Голая сглаживают все данные, которые они фильтруют. Однако это не всегда может быть желаемым. Например, что если наши данные взяты из тактового сигнала и имеют острые края, которые мы не хотим сглаживать? Фильтры, обсуждавшиеся до сих пор, работают не так хорошо:

load clockex

yMovingAverage = conv(x,ones(5,1)/5,'same');
ySavitzkyGolay = sgolayfilt(x,3,5);
plot(t,x, ...
     t,yMovingAverage, ...
     t,ySavitzkyGolay)

legend('original signal','moving average','Savitzky-Golay')

Скользящее среднее и Savitzky-Golay фильтруют соответственно недостаточно корректно и чрезмерно корректно вблизи краев тактового сигнала.

Простой способ сохранить кромки, но все же сгладить уровни - использовать медианный фильтр:

yMedFilt = medfilt1(x,5,'truncate');
plot(t,x, ...
     t,yMedFilt)
legend('original signal','median filter')

Удаление отклонений через фильтр Хампеля

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

Чтобы увидеть это, загрузите аудиозапись свистка поезда и добавьте несколько искусственных скачков шума:

load train
y(1:400:end) = 2.1;
plot(y)

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

hold on
plot(medfilt1(y,3))
hold off
legend('original signal','filtered signal')

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

hampel(y,13)
legend('location','best')

Из исходного сигнала удаляются только отклонения.

Дальнейшее чтение

Для получения дополнительной информации о фильтрации и повторной выборке см. Панель инструментов обработки сигналов.

Справка: Кендалл, Морис Г., Алан Стюарт и Дж. Кит Орд. Продвинутая теория статистики, том 3: Дизайн и анализ, и Временной ряд. 4-й эд. Лондон: Макмиллан, 1983.