Документация Engee
Notebook

Фильтрация 1-D данных

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

В этом примере мы загрузим мат-файл с идеальным сигналом и шумной версией, а потом последовательно применим несколько популярных методов сглаживания и разложения: скользящее среднее, медианный фильтр, фильтр Савицкого–Голея, фильтр Хампеля, LOESS-сглаживание, TV-денойз, сингулярный спектральный анализ и эмпирическое разложение (EMD).

Подготовка к работе

Нужны: визуализация, чтение .mat, разные алгоритмы обработки и декомпозиции. Для этого импортируем необходимые пакеты

In [ ]: 
include("$(@__DIR__)/InstallPackages.jl")
In [ ]: 
using Plots
using MAT
using RollingFunctions
using SavitzkyGolay
using HampelOutliers
using Loess
using TVDenoise
using SingularSpectrumAnalysis
using EMD

Читаем файл, достаем матрицу result, из неё — эталон Y и шумный измеренный сигнал X. Сразу сравним их

In [ ]: 
# Чтение данных
path = "$(@__DIR__)/file.mat"
mat = matopen(path)     
data = read(mat)     
data = data["result"];      
Y = data["thetTgt"]
X = data["thetTgt_MFit"];
# Визуализация данных
plot(X', label="Шумные данные")
plot!(Y', label="Эталонные данные")
Out[0]:

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

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

In [ ]: 
nan_idxs = findall(isnan, X)
println("Количесвто пропусков в массиве Y = ", sum(findall(isnan, Y))[1],  " Количесвто пропусков в массиве X = ", sum(findall(isnan, X))[1])
ff = findfirst(x -> !isnan(x), X)
ll = findlast(x -> !isnan(x), X)
X = X[ff:ll]
Y = Y[ff:ll]

Количесвто пропусков в массиве Y = 0 Количесвто пропусков в массиве X = 9
Out[0]:
1×399 Matrix{Float64}:
 29.3893  29.4091  29.4282  29.4466  …  27.2312  27.2197  27.2081  27.1965

Алгоритм скользящего среднего

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

In [ ]: 
w = 10
x = vec(X)
raw = rollmean(x, w)       
pad = fill(NaN, w-1)          
meanX = vcat(pad, raw);   
plot([Y' vec(X) meanX], label=["Исходный сигнал" "Шумовой сигнал" "Отфильтрованный сигнал"], ylabel="", xlabel="t", legend=:topleft)
Out[0]:

Видно, что шум становится сглаживаются. Большие выбросы, конечно, никуда не делись. Этот метод самый простой.

Фильтр Савицкого–Голея

Фильтр Савицкого–Голея — это метод сглаживания, при котором для каждого положения окна фиксированной нечётной длины строится наименьшими квадратами полином заданной степени и вычисляются его коэффициенты. Затем значение в центральной точке заменяется предсказанием этого полинома.

In [ ]: 
y2_sg = savitzky_golay(vec(X), 101, 5);
plot([Y' vec(X) y2_sg.y], label=["Исходный сигнал" "Шумовой сигнал" "Отфильтрованный сигнал"], ylabel="", xlabel="t", legend=:topleft)
Out[0]:

Хампел-фильтр

Хампел-фильтр похож на медианный. Для каждого окна фиксированного размера он сначала находит медиану и оценивает разброс через медианное абсолютное отклонение (MAD), а потом проверяет каждую точку: если её расстояние до медианы больше заданного порога , то вместо неё ставят медиану окна, иначе оставляют без изменений. Так фильтр надёжно устраняет выбросы.

In [ ]: 
yHampel = Hampel.filter(vec(X), 30; boundary=:reflect, threshold=0);
plot([Y' vec(X) yHampel], label=["Исходный сигнал" "Шумовой сигнал" "Отфильтрованный сигналl"], ylabel="", xlabel="t", legend=:topleft)
Out[0]:

LOESS-сглаживание

LOESS-сглаживание на каждом шаге для окрестности точки строит локальную регрессию, которая предсказывает новое значение в центральной точке.

In [ ]: 
N = length(vec(X))
t = collect(1.0:N)
model = loess(reshape(t, N, 1), vec(X); span=0.2, degree=1);
X_loess = predict(model, t);
plot([Y' vec(X) X_loess], label=["Исходный сигнал" "Шумовой сигнал" "Отфильтрованный сигналl"], ylabel="", xlabel="t", legend=:topleft)
Out[0]:

Total Variation Denoising

Total Variation Denoising минимизирует разрывы и шум, решая задачу оптимизации - найти новый сигнал, близкий по L2-норме к исходному, но с минимальной общей вариацией.

In [ ]: 
kw = Dict(:isotropic=>false, :maxit=>100, :tol=>1e-4, :verbose=>false)
λ = 30.0
r = 3.0

denoised1d, hist1d = tvd(vec(X), λ, r;
                         kw...);
plot([Y' vec(X) denoised1d], label=["Исходный сигнал" "Шумовой сигнал" "Отфильтрованный сигналll"], ylabel="", xlabel="t", legend=:topleft)
Out[0]:

Сингулярный спектральный анализ (SSA)

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

Визуализируем извлеченный тренд

In [ ]: 
L = 40
trend, seasonal = analyze(vec(X), L)
plot(vec(X),         label="Шумный сигнал", alpha=0.5)
plot!(trend,     label="Тренд",       lw=2)
plot!(Y', label="Эталон", lw=2)
Out[0]:

Эмпирическое разложение

Эмпирическое модальное разложение разбивает сигнал на набор внутренних модальных функций — компоненты с осцилляциями разной частоты.

In [ ]: 
imf = EMD.emd(vec(X))


t = 1:length(vec(X))
plot(t, imf[7,:], label="IMF 7")
plot!(Y', label="Эталон", lw=2)
Out[0]:

Выводы

В данном примере были разобраны различные методы предобработки сигналов (фильтрации).