Мотивация

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

Иногда, когда вы исследуете входные данные, можно хотеть сглаживать данные для того, чтобы видеть тренд в сигнале. В нашем примере у нас есть набор температурных показаний в Цельсия, бравшемся каждый час в Аэропорту Логана в Бостоне в течение целого месяца января 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)')

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

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

Фильтры Savitzky-Golay

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

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

Как удобство, можно использовать функциональный sgolayfilt реализовать Savitzky-Golay сглаживание фильтра. Использовать 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')

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

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

Существуют примерно 1 000 / 60 = 16,667 выборок в полном цикле 60 Гц, когда произведено на уровне 1 000 Гц. Давайте попытаемся "окружить" и использовать фильтр с 17 точками. Это даст нам максимальную фильтрацию на основной частоте 1 000 Гц / 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 Гц = 1 020 Гц, мы можем использовать наши 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')

Средний фильтр

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

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

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

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

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

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

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

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

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

Только выбросы удалены из исходного сигнала.

Дополнительные материалы для чтения

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

Ссылка: Кендалл, Морис Г., Алан Стюарт и порядок Дж. Кита. Усовершенствованная теория статистики, издания 3: проект и анализ и timeseries. 4-й Эд. Лондон: Макмиллан, 1983.