Engee 文档
Notebook

过滤1-D数据

数据过滤是训练神经网络之前的重要步骤。 如果没有干净的信号，模型可以从噪声和异常值中学习，这将大大降低预测的质量。 经过过滤后，数据将变得更加干净，神经网络将能够专注于主要的事情-那些真正重要的模式。

在这个例子中，我们将上传一个具有理想信号和噪声版本的mat文件，然后依次应用几种流行的平滑和分解方法：移动平均、中值滤波器、Savitsky-Golay滤波器、Hampel滤波器、黄土平滑、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]:

可以看出，噪声信号在参考信号周围波动非常强烈。 平滑RF噪声波动会很好。

为了使过滤算法正常工作，数据中没有间隙是必要的。 在测量数据中经常发现信号边缘的间隙。 如果它们存在于所研究的信号中，只需取一部分数据即可。

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]:

可以看出，噪声正在变得平滑。 当然，大量排放并没有消失。 这种方法是最简单的。

Savitsky-Goley过滤器

Savitsky-Goley滤波器是一种平滑方法，其中为固定奇数长度的窗口的每个位置构造给定程度的最小二乘多项式并计算其系数。 然后将中心点处的值替换为此多项式的预测。

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

汉普尔过滤器

Hampel滤波器类似于中值滤波器。 对于每个固定大小的窗口，它首先找到中位数并通过中位数绝对偏差（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]:

黄土-平滑

在每个步骤中，黄土平滑都会为点的邻域构建局部回归，从而预测中心点处的新值。

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]:

总变异去噪

总变异去噪通过解决寻找接近原始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]:

结论

在这个例子中，分析了各种信号预处理（滤波）方法。