hampel

Удаление выбросов с помощью идентификатора Хампеля

Описание

пример

y = hampel(x) применяет фильтр Hampel к входному вектору, x, обнаружить и удалить выбросы. Для каждой выборки x, функция вычисляет медиану окна, состоявшего из выборки и ее шести окружающих выборок, три на сторону. Это также оценивает стандартное отклонение каждой выборки о ее медиане окна использование среднего абсолютного отклонения. Если выборка отличается от медианы больше чем тремя стандартными отклонениями, это заменяется медианой. Если x матрица, затем hampel обработки каждый столбец x как независимый канал.

y = hampel(x,k) задает количество соседей, k, по обе стороны от каждой выборки x в окне измерения. k значения по умолчанию к 3.

пример

y = hampel(x,k,nsigma) задает много стандартных отклонений, nsigma, которым выборка x должен отличаться от локальной медианы для него, чтобы быть замененным медианой. nsigma значения по умолчанию к 3.

пример

[y,j] = hampel(___) также возвращает логическую матрицу, которая верна в местоположениях всех точек, идентифицированных как выбросы. Этот синтаксис принимает любой из входных параметров от предыдущих синтаксисов.

пример

[y,j,xmedian,xsigma] = hampel(___) также возвращает локальные медианы и предполагаемые стандартные отклонения для каждого элемента x.

hampel(___) без выходных аргументов строит отфильтрованный сигнал и аннотирует выбросы, которые были удалены.

Примеры

свернуть все

Сгенерируйте 100 выборок синусоидального сигнала. Замените шестые и двадцатые выборки на скачки.

x = sin(2*pi*(0:99)/100);
x(6) = 2;
x(20) = -2;

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

[y,i,xmedian,xsigma] = hampel(x);

Постройте отфильтрованный сигнал и аннотируйте выбросы.

n = 1:length(x);
plot(n,x)
hold on
plot(n,xmedian-3*xsigma,n,xmedian+3*xsigma)
plot(find(i),x(i),'sk')
hold off
legend('Original signal','Lower limit','Upper limit','Outliers')

Figure contains an axes. The axes contains 4 objects of type line. These objects represent Original signal, Lower limit, Upper limit, Outliers.

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

hampel(x,1)

Figure contains an axes. The axes contains 3 objects of type line. These objects represent original signal, filtered signal, outliers.

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

rng('default')

n = 59;
x = sin(pi./[15 10]'*(1:n)+pi/3)';

spk = randi(2*n,9,1);
x(spk) = x(spk)*2;
x(randi(2*n,6,1)) = NaN;

plot(x)

Figure contains an axes. The axes contains 2 objects of type line.

Отфильтруйте сигнал с помощью hampel с настройками по умолчанию.

y = hampel(x);
plot(y)

Figure contains an axes. The axes contains 2 objects of type line.

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

y = hampel(x,4,2);
plot(y)

Figure contains an axes. The axes contains 2 objects of type line.

Выведите рабочую медиану для каждого канала. Наложите медианы на графике сигнала.

[y,j,xmd,xsd] = hampel(x,4,2);
plot(x)
hold on
plot(xmd,'--')

Figure contains an axes. The axes contains 4 objects of type line.

Сгенерируйте многоканальный сигнал, который состоит из двух синусоид различных частот, встроенных в белый Гауссов шум модульного отклонения.

rng('default')

t = 0:60;
x = sin(pi./[10;2]*t)'+randn(numel(t),2);

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

k = 4;
nsig = 2;

[y,h] = hampel(x,k,nsig);

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

for k = 1:2
    hk = h(:,k);
    ax = subplot(2,1,k);
    plot(t,x(:,k))
    hold on
    plot(t,y(:,k))
    plot(t(hk),x(hk,k),'*')
    hold off
    ax.XTick = t(hk);
end

Figure contains 2 axes. Axes 1 contains 3 objects of type line. Axes 2 contains 3 objects of type line.

Сгенерируйте 100 выборок синусоидального сигнала. Замените шестые и двадцатые выборки на скачки.

n = 1:100;
x = sin(2*pi*n/100);
x(6) = 2;
x(20) = -2;

Используйте hampel вычислить локальное среднее и предполагаемое стандартное отклонение для каждой выборки. Используйте значения по умолчанию входных параметров:

  • Размер окна 2×3+1=7.

  • Вопросы, которые отличаются от их медианы окна больше чем тремя стандартными отклонениями, рассматриваются выбросы.

Постройте результат.

[y,i,xmedian,xsigma] = hampel(x);

plot(n,x)
hold on
plot(n,[1;1]*xmedian+3*[-1;1]*xsigma)
plot(find(i),x(i),'sk')
hold off
legend('Signal','Lower','Upper','Outliers')

Figure contains an axes. The axes contains 4 objects of type line. These objects represent Signal, Lower, Upper, Outliers.

Повторите вычисление с помощью размера окна 2×10+1=21 и два стандартных отклонения как критерии идентификации выбросов.

sds = 2;
adj = 10;
[y,i,xmedian,xsigma] = hampel(x,adj,sds);

plot(n,x)
hold on
plot(n,[1;1]*xmedian+sds*[-1;1]*xsigma)
plot(find(i),x(i),'sk')
hold off
legend('Signal','Lower','Upper','Outliers')

Figure contains an axes. The axes contains 4 objects of type line. These objects represent Signal, Lower, Upper, Outliers.

Входные параметры

свернуть все

Входной сигнал в виде вектора или матрицы. Если x матрица, затем hampel обработки каждый столбец x как независимый канал.

Пример: cos(pi/4*(0:159))+randn(1,160) одноканальный сигнал вектора-строки.

Пример: cos(pi./[4;2]*(0:159))'+randn(160,2) двухканальный сигнал.

Типы данных: single | double

Количество соседей по обе стороны от демонстрационного xs в виде целочисленного скаляра. Выборки близко к ребрам сигнала, которые имеют меньше, чем k выборки на одной стороне сравниваются с медианой меньшего окна.

Типы данных: single | double

Количество стандартных отклонений, который выборка x должен отличаться от его локальной медианы, которая будет рассмотрена выбросом. Задайте nsigma как действительный скаляр. Функция оценивает стандартное отклонение путем масштабирования локального среднего абсолютного отклонения (MAD) на коэффициент κ=12erf1(1/2)1.4826.

Типы данных: single | double

Выходные аргументы

свернуть все

Фильтрованный сигнал, возвращенный как вектор или матрица одного размера с x.

Типы данных: single | double

Индекс выброса, возвращенный как вектор или матрица одного размера с x.

Типы данных: логический

Локальные медианы, возвращенные как вектор или матрица одного размера с x.

Типы данных: single | double

Предполагаемые стандартные отклонения, возвращенные как вектор или матрица одного размера с x.

Типы данных: single | double

Больше о

свернуть все

Идентификатор Хампеля

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

Учитывая последовательность x 1, x 2, x 3, …, xn и раздвижное окно длины k, задает оценочное использование среднего и стандартного отклонения "точка-точка":

  • Локальная медиана — mi=median(xik,xik+1,xik+2,,xi,,xi+k2,xi+k1,xi+k)

  • Стандартное отклонение σi=κmedian(|xikmi|,,|xi+kmi|), где κ=12erf1(1/2)1.4826

Количество σi/κ известно как среднее абсолютное отклонение (MAD).

Если демонстрационный xi таков что

|ximi|>nσσi

для заданного порога затем идентификатор Хампеля объявляет xi выброс и заменяет его на mi.

Около конечных точек последовательности функция обрезает окно, использованное для расчета mi и σi:

  • i <k + 1

    mi=median(x1,x2,x3,,xi,,xi+k2,xi+k1,xi+k)

    σi=κmedian(|x1m1|,,|xi+kmi|)

  • i> nk

    mi=median(xik,xik+1,xik+2,,xi,,xn2,xn1,xn)

    σi=κmedian(|xikmi|,,|xnmn|)

Ссылки

[1] Лю, Hancong, Сириш Шах и Вэй Цзян. “Онлайновое определение выбросов и очистка данных”. Компьютеры и Химическое машиностроение. Издание 28, март 2004, стр 1635–1647.

Введенный в R2015b
Для просмотра документации необходимо авторизоваться на сайте