Медианный фильтр для удаления выбросов
В этой третьей части обработки сигналов с помощью Python я бы обсудил использование медианного фильтра для удаления больших пиковых сигналов. Вот ссылки на первую и вторую части. Во-первых, давайте еще раз проясним, в чем разница между средним и медианой ряда.
Чтобы вычислить медиану набора чисел, нам нужно
Медиана
- Сортировать числа по возрастанию
- a. Если есть нечетное количество значений в списке, медиана - это просто среднее число, b. Если есть четное количество значений в списке, медиана среднее из двух средних значений
Теперь рассмотрим список чисел. Как видите, большинство из этих чисел в этом конкретном списке меньше 50, тогда как очень немногие числа действительно большие. Среднее этих чисел равно
79.54545454545455
тогда как значение медианы
27.0
Таким образом, среднее значение набора сильно подтягивается несколькими большими значениями в списке, которые могут быть выбросами (нерепрезентативные точки данных). В отличие от этого, на медианное значение в значительной степени не влияет наличие выбросов.
Итак, когда мы применяем медианный фильтр к определенной точке данных, мы в основном заменяем эту точку данных на медианное значение ее соседних точек в пределах вручную установленного окна.
Прежде чем я подробнее расскажу о медианном фильтре, я хочу добавить, что медиана - это нелинейный фильтр. Таким образом, лучше всего применять медианный фильтр только к ограниченному количеству точек данных в вашем сигнале. Эти точки данных в большинстве случаев представляют собой необычные точки данных, выбросы или отсутствующие сигналы. Подробнее о линейных и нелинейных фильтрах я расскажу в одном из будущих постов.
Коды на Python
Давайте перейдем к кодированию здесь. Импортируйте необходимые библиотеки numpy и matplotlib. % Matplotlib inline позволяет визуализировать необходимые графики в записной книжке Jupyter.
import numpy as np import matplotlib.pyplot as plt import copy %matplotlib inline
Теперь сначала мы создаем набор случайных данных. Чтобы сгенерировать сигнал, мы будем использовать совокупную сумму списка случайных чисел, которая также известна как броуновское движение / случайное блуждание. В Numpy уже есть для этого встроенная функция cumsum (). Однако в более позднем посте я мог бы обсудить, как определить функцию для имитации случайного блуждания, просто для удовольствия.
#create signal #number of data points n = 1000 signal = np.cumsum(np.random.randn(n)) plt.plot(signal
Теперь я собираюсь установить часть точек данных на некоторое большое количество значений, которые будут отображаться как пиковый сигнал.
#find noise points #set proportion of data points for spiked noise prop_noise = 0.05 #randomize the indices noise_points = np.random.permutation(len(signal)) noise_points = noise_points[0:int(n*prop_noise)] #replace some proportion of time points with large spikes signal[noise_points] = 100 + np.random.randn(len(noise_points))*50 plt.plot(signal)
Вы должны иметь в виду, что фильтр сглаживания среднего или сглаживающий фильтр по Гауссу здесь не принесет особой пользы, потому что среднее значение будет в значительной степени взвешено по отношению к выбросам, которые имеют большие (положительные или отрицательные) значения. Но здесь наша цель - заменить выбросы некоторой комбинацией соседних точек. Именно для этого и предназначен медианный фильтр.
Затем мы должны решить, как определить выбросы для конкретного набора данных. Эта информация может быть заранее известна из вашего эксперимента, например, в каком диапазоне вы должны ожидать значения. Однако здесь мы собираемся определить порог выбросов путем визуального осмотра. Давайте посмотрим на гистограмму для всех значений данных, которые у нас есть.
plt.hist(signal, bins = 30, color='#0504aa',
alpha=0.7, rwidth=0.85)
Теперь, если мы увеличим масштаб оси x, мы обнаружим, что большинство значений в наборе данных лежат ниже 50. Таким образом, мы устанавливаем любое значение выше 50 как выброс.
#define threshold threshold = 50 # find data values above the threshold outliers = np.where( signal>threshold )[0]
Теперь я инициализирую отфильтрованный сигнал копией исходного сигнала.
# initialize filtered signal filtsig = copy.deepcopy(signal)
Теперь я просматриваю временные точки в серии, где значение превышает пороговое значение. Я заменяю эти значения в исходном сигнале средним значением от -k до + k соседних значений. Я установил здесь границу, потому что может случиться так, что одна из действительно ранних (в нашем случае ‹20-ая) или одна из последних точек данных (в нашем случае› (n-20) -я) имеет всплеск. В этом случае большое значение + 20 (или -20) будет вне диапазона нашего набора данных. Чтобы предотвратить эту ошибку в коде, мы устанавливаем нижнюю границу либо как 0, либо как выбросы [i] -k, в зависимости от того, что больше. Точно так же верхняя граница - выбросы [i] + k или n + 1, в зависимости от того, что меньше.
k=20
for i in range(len(outliers)):
# lower and upper bounds
lowbound = np.max((0,outliers[i]-k))
upbound = np.min((outliers[i]+k,n+1))
# compute median of surrounding points
filtsig[outliers[i]] = np.median(signal[lowbound:upbound])
Давайте теперь построим вместе исходный сигнал и обработанный сигнал.
Итак, мы в основном заменили пиковые точки данных на среднее значение их соседей, сохранив исходный сигнал без изменений. Именно для этого и полезны медианные фильтры.
