Мотивация

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

Иногда, когда вы исследуете входные данные, вы можете сглаживать данные в порядок, чтобы увидеть тренд в сигнале. В нашем примере у нас есть набор показаний температуры в Цельсии, взятых каждый час в аэропорту Логан в Бостоне за весь месяц января 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[1/2,1/2\right]}^{\mathit{n}}$. Этот тип фильтра аппроксимирует нормальную кривую для больших значений n. Это полезно для отфильтровывания высокой частоты шума для малых n. Чтобы найти коэффициенты для биномиального фильтра, сверните $\left[1/2,1/2\right]$ с самим собой и затем итерационно свернуть выход с $\left[1/2,1/2\right]$ заданное количество раз. В этом примере используйте пять общих итераций.

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, вы задаете сегмент данных нечетной длины и полином порядка строго меньше, чем длина сегмента. The 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')

Удаление выбросов через фильтр Hampel

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

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

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')

Из исходного сигнала удаляются только выбросы.

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

Для получения дополнительной информации о фильтрации и повторной дискретизации смотрите Signal Processing Toolbox.

Ссылка: Кендалл, Морис Г., Алан Стюарт и Дж. Кит Орден. Передовая теория статистики, том 3: Проект и анализ, и Timeseries. 4th Ed. London: Macmillan, 1983.