Фильтр Гаусса к временному ряду

В наши дни Интернет наводнен ресурсами, которые помогают ориентироваться в науке о данных и/или машинном обучении. Часто вы могли бы найти младшего ученого, такого как я, погруженного в кучу данных и пытающегося разобраться в них (что вы также можете назвать анализом временных рядов, используя причудливые термины). Только когда я рыскал в Интернете в поисках руководств по анализу научных данных, я заметил нехватку ресурсов. Таким образом, я подумал, что было бы здорово, если бы я смог поделиться некоторыми методами, которые я изучил/изучаю на протяжении всего курса процедуры анализа данных/сигналов с использованием Python в качестве языка программирования.

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

  1. Анализ научных данных (спектральный анализ сигналов для спектроскопии/биологических наук)
  2. Обработка аудиосигнала, цифровая обработка сигнала

3. Обработка речевого сигнала,

4. Обработка изображений

5. Обработка финансовых данных

6. Сейсмология и т. д.

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

В этой конкретной серии (следите за обновлениями для будущих частей, часть 1, часть 3) я хочу затронуть темы, касающиеся нескольких инструментов, необходимых для обработки сигналов. Они перечислены ниже.

  • Временная фильтрация
  • Свертка
  • Вейвлет-преобразование
  • Спектральный анализ (быстрое преобразование Фурье)
  • Анализ частотно-временной области
  • Очистка/шумоподавление данных
  • Повторная выборка (вверх и вниз)
  • Интерполяция и экстраполяция
  • Обнаружение функций
  • Отношение сигнал/шум

Ядро Гаусса (теория):

В части I серии я говорил о фильтре среднего ряда, где мы рассматривали каждую точку данных как среднее число k (порядок) соседних точек данных. Здесь мы собираемся использовать функцию Гаусса в качестве ядра для сглаживания, где для каждой точки данных мы будем взвешивать предыдущую и последующую точки данных по кривой Гаусса. Здесь нет коэффициента нормализации, как это было для среднего фильтра. Это связано с тем, что мы уже предполагаем, что весовая функция gнормализована. Другими словами, площадь под функцией Гаусса равна одной единице.

Давайте посмотрим на стандартную функцию Гаусса. Вот формула для создания функции Гаусса.

Наиболее важным параметром гауссовой функции является атрибут, называемый полной шириной и половиной максимума (FWHM), w. FWHM гауссовой диаграммы — это расстояние между двумя точками, наиболее близкими к 50-процентному усилению. Полуширина связана со стандартным отклонением гауссовой диаграммы следующим образом.

Кодирование на Python

Теперь позвольте мне показать, как выглядит код в Python. Здесь я предполагаю базовый уровень знакомства читателей с python.

Импорт библиотек

Во-первых, давайте импортируем библиотеки, необходимые для запуска кодов.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Сгенерируйте искусственный сигнал с шумом

sigRate = 2000 #Hz
time = np.arange(0,3, 1/sigRate)
n = len(time)
p = 15 #poles for random interpolation
ampl = np.interp(np.linspace(0,p,n),np.arange(0,p),np.random.rand(p)*30)
noiseamp = 5
noise = noiseamp*np.random.randn(n)
signal = ampl + noise
plt.plot(time, ampl, 'ro-', markersize = 5)
plt.plot(time, signal, alpha=0.5)

Теперь мы генерируем ядро ​​Гаусса. gauss_time — это вектор времени, который мы используем для фильтра. Значение k должно быть таким, чтобы гауссиана опускалась до нуля на обоих своих крыльях.

fwhm = 25 # in ms
k = 50
gauss_time = 2000*np.arange(-k,k)/sigRate
# create Gaussian window
gauswin = np.exp( -(4*np.log(2)*gauss_time**2) / fwhm**2 )

Например, если k слишком мало (по сравнению с fwhm), крылья на самом деле не распрямляются.

Тогда как для k = 25 он выглядит гладким. Подробнее об этом краевом эффекте мы поговорим в одной из наших будущих публикаций.

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

pstPeakHalf = k + np.argmin( (gauss_win[k:]-.5)**2 )
prePeakHalf = np.argmin( (gauss_win-.5)**2 )
empFWHM = gauss_time[pstPeakHalf] - gauss_time[prePeakHalf]

Остальное очень похоже на приложение среднего фильтра, где мы создаем переменную filtSig и инициализируем ее нулем. Затем примените к нему фильтр Гаусса.

# initialize filtered signal vector
filtSig_Gauss = np.zeros(n)
# # implement the running mean filter
for i in range(k+1,n-k-1):
    # each point is the weighted average of k surrounding points
    filtSig_Gauss[i] = np.sum( signal[i-k:i+k]*gauss_win )
    
plt.plot(time,signal,'r',label='Original')
plt.plot(time,filtSig_Gauss,'k',label='Gaussian-filtered')
plt.xlabel('Time (s)')
plt.ylabel('amp. (a.u.)')
plt.legend()
plt.title('Gaussian smoothing filter')

Подробнее о выборе фильтра я расскажу в одном из следующих постов.