AnyMath 文档

斯戈莱

Savitsky-Goley滤波器的计算。

库::`工程师`

语法

函数调用

* [参数:b],[参数:g]=sgolay(<参数:m>>,<参数:fl>>) —用多项式的阶数设计Savitsky-Goley平滑FIR滤波器 [参数:m] 和帧的长度 [参数:fl]. 返回FIR滤波器系数矩阵 [参数:b] 和差分滤波器的矩阵 [参数:g].

* [参数:b],[参数:g]=sgolay(<参数:m>>,<参数:fl>>,<参数:w>>) -也使用向量 [参数:w]，其中包含最小二乘最小化中使用的实正权重。

争论

输入参数

# *m*是 多项式的阶数

+ 一个非负整数

Details

多项式的阶数，作为非负整数给出。 意义 m 必须小于 [参数:fl]. 如果 M=<参数:fl>>-1，则所设计的滤波器不执行平滑。

# *fl* — 框架长度

+ 一个正奇数整数

Details

帧的长度，设为正奇数整数。 意义 弗尔 一定有超过 [参数:m].

# *w*是 权重的向量

+ 一(fl,1) （默认情况下）| 真正的正向量

Details

定义为具有实数正元素的向量的权重的向量。 权重向量的长度与 [参数:fl]，并用于最小二乘最小化。

输出参数

# *b* — 随时间变化的FIR滤波器系数

+ 矩阵

Details

Fir滤波器系数，时变，作为矩阵返回 [参数:fl]×[参数:fl]. 在平滑滤波器的实施中，最后 (<参数:fl>>-1)/2 行（每个都有一个滤波器）在初始瞬态期间应用于信号，并且第一 (<参数:fl>>-1)/2 行-在最后的过渡过程中。 中心线施加稳态信号。

# *g*是 差分滤波器的矩阵

+ 矩阵

Details

差分滤波器的矩阵。 每列 g 它是阶导数的差分滤波器 p-1，在哪里 p -列的索引。 如果设置了信号 x 长度 [参数:fl] 你可以找到导数的估计 p-顺序, xp碌鲁，其平均值从 xp((fl+1)/2)=阶乘(p)*(g[:,p+1]'*x).

例子：

使用Savitsky-Goley方法平滑噪声正弦

Details

我们将产生一个信号，代表一个频率为 0.2 Hz，以每秒五次的采样率将白高斯噪声添加到其中。 200 几秒钟。

import EngeeDSP.Functions: randn, sgolay,conv

dt = 1/5
t = (0:dt:200-dt)'
x = 5&ast;sin.(2&ast;pi&ast;0.2&ast;t) + randn(1, length(t))

我们使用 *sgolay* 来平滑信号。 我们使用框架从 21 计数和第四度的多项式。

order = 4
framelen = 21
b,g = sgolay(order,framelen)

通过与向量的中间行卷积来计算信号的稳态部分 b.

ycenter = conv(x,b[Int((framelen+1)/2), :],"valid")'

让我们计算瞬态。 我们使用矩阵的最后一行 b 开始和矩阵的第一行 b 完成。

ybegin = b[end:-1:Int((framelen+3)/2), :] &ast; x[framelen:-1:1]
yend = b[Int((framelen-1)/2):-1:1, :] &ast; x[end:-1:end-(framelen-1)]

让我们将信号的瞬态和稳态部分结合起来，得到一个完整的平滑信号。 让我们绘制原始信号和使用Savitsky—Goley方法平滑的信号。

y = [ybegin; ycenter; yend]

plot(x', label = "Noisy Sinusoid")
plot!(y, label = "S-G smoothed sinusoid")

sgolay 1

Savitsky-Goley过滤器与尺度

Details

让我们设计一个具有多项式的Savitsky—Goley平滑滤波器 5`以过滤频率为 `8192 赫兹。 我们用汉明窗。 帧长度为 101 倒计时。

import EngeeDSP.Functions: hamming, freqz, sgolay

Fs = 8192
m = 5
fl = 101
weights = hamming(fl)
b,g = sgolay(m,fl,weights)

我们来构造中心线的频率响应 b. 我们使用 1024 FFT点。

NFFT = 1024
b=b[Int((fl+1)/2), :]
freqz(b',1,NFFT,Fs,out=:plot)

sgolay 2

萨维茨基至戈利航线

Details

我们将产生一个信号，代表一个频率为 0.2 Hz，其中白高斯噪声被添加，以每秒四次的采样率在 20 几秒钟。

import EngeeDSP.Functions: randn, sgolay,conv

dt = 0.25
t = (0:dt:20-1)'
x = 5&ast;sin.(2&ast;pi&ast;0.2&ast;t) + 0.5&ast;randn(1, length(t))

让我们使用Savitsky—Goley方法评估正弦波的前三个导数。 我们使用 25-点窗口和五阶多项式。 让我们将列划分为度 dt的 以正确缩放导数。

b,g = sgolay(5,25)

dx = zeros(length(x), 4)
for p in 0:3
    kernel = factorial(p) / (-dt)^p &ast; g[:, p+1]
    dx[:, p+1] = conv(x, kernel, "same")
end

让我们绘制原始信号，平滑序列和导数。

plot(x', marker = :., label = "x")
plot!(dx, label = ["x (smoothed)" "x'" "x''" "x'''"])
title!("Savitzky-Golay Derivative Estimates")

sgolay 3

算法

Savitsky-Goley滤波器用于平滑具有较大频率范围的噪声信号。 Savitsky-Goley平滑滤波器通常比标准FIR平均滤波器滤除信号的高频内容更少。 然而，当噪声电平特别高时，它们在抑制噪声方面的效果较差。

一般而言，滤波包括用包含在以该点为中心的滑动窗口中的信号值的某种组合来替换每个信号点，基于附近点测量几乎相同的基值的假设。 例如，移动平均滤波器将每个数据点替换为周围数据点的局部平均值。 如果数据点有 左边的点和 右边的点，表示窗口的总长度 ，则移动平均滤波器进行替换


Savitsky-Goley滤波器通过使用最小二乘法近似阶多项式来推广这个想法。 到窗口中信号的值，并将计算得到的多项式曲线中心点作为新的平滑数据点。 对于给定的点

或者，就矩阵而言,


为了找到Savitsky-Goley估计值，我们使用pseudoinversion 来计算 ，然后在左边乘以 :

为了避免不良调理，功能 *sgolay* 使用QR分解计算经济分解 如何 ，就其而言 .

计算方法 这是必要的只有一次。 Savitsky-Golay对大部分信号点的估计是通过将信号与中心线卷积得到的。 . 结果是滤波信号的稳态部分。 第一个 排队 他们给出了最初的过渡过程，最后 排队 -最后的过渡过程。 降噪可以通过增加窗口长度来改善，但这会导致任何瞬变附近的吉布斯样振铃。

文学作品

  1. Orfanidis，Sophocles J.Introduction To Signal Processing. 恩格尔伍悬崖，新泽西州：Prentice Hall，1996。

  2. Press，William H.，Saul A.Teukolsky，William T.Vetterling和Brian P.Flannery。 C：科学计算的艺术。 繝ｼ繝ｫ縺ｧ縺呐

  3. Schafer，Ronald W."What_Is_A Savitzky-Golay Filter？ [讲义]。"_IEEE信号处理杂志_卷。 28，第4号，2011年7月，第111-117页。 https://doi.org/10.1109/MSP.2011.941097.