sgolay

R2026b

萨维茨基-戈雷滤波器设计

全页折叠

语法

b = sgolay(m,fl)

b = sgolay(m,fl,w)

[b,g] = sgolay(___)

描述

b = sgolay(m,fl) 设计具有多项式阶数 m 和帧长度 fl 的萨维茨基-戈雷 FIR 平滑滤波器。

示例

b = sgolay(m,fl,w) 指定加权向量 w，其中包含在最小二乘最小化期间要使用的正实数值权重。

示例

[b,g] = sgolay(___) 返回微分滤波器的矩阵 g。您可以将这些输出参量与任何先前的输入语法结合使用。

示例

全部折叠

含噪正弦波的萨维茨基-戈雷平滑

打开实时脚本

生成信号，该信号由嵌入在高斯白噪声中的 0.2 Hz 正弦波组成，并以每秒五次的采样率进行 200 秒采样。

dt = 1/5;
t = (0:dt:200-dt)';

x = 5*sin(2*pi*0.2*t) + randn(size(t));

使用 sgolay 对信号进行平滑。使用包含 21 个采样的帧和四阶多项式。

order = 4;
framelen = 21;

b = sgolay(order,framelen);

通过将信号与 b 的中心行进行卷积来计算信号的稳态部分。

ycenter = conv(x,b((framelen+1)/2,:),'valid');

计算瞬变。使用 b 的最后几行用于启动阶段，使用 b 的前几行用于终止阶段。

ybegin = b(end:-1:(framelen+3)/2,:) * x(framelen:-1:1);
yend = b((framelen-1)/2:-1:1,:) * x(end:-1:end-(framelen-1));

将瞬变部分与稳态部分串联起来以生成完整平滑信号。绘制原始信号和萨维茨基-戈雷估计。

y = [ybegin; ycenter; yend];
plot([x y])
legend('Noisy Sinusoid','S-G smoothed sinusoid')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Noisy Sinusoid, S-G smoothed sinusoid.

带权重的萨维茨基-戈雷滤波器

打开实时脚本

设计具有五阶多项式的萨维茨基-戈雷平滑滤波器，对以 $F_{s} = 8192$ Hz 采样的含噪信号进行滤波。使用汉明窗。帧长度为 101 个采样。

Fs = 8192;
m = 5;
fl = 101;

weights = hamming(fl);

b = sgolay(m,fl,weights);

绘制 b 的中心行的频率响应。使用 1024 个 FFT 点。

NFFT = 1024;
freqz(b((fl+1)/2,:),1,NFFT,Fs)

Figure contains 2 axes objects. Axes object 1 with title Phase, xlabel Frequency (kHz), ylabel Phase (degrees) contains an object of type line. Axes object 2 with title Magnitude, xlabel Frequency (kHz), ylabel Magnitude (dB) contains an object of type line.

有关加权向量及其对信号滤波性能影响的详细信息，请参阅Design and Analyze Savitzky-Golay Filters。

萨维茨基-戈雷微分

打开实时脚本

生成信号，该信号由嵌入在高斯白噪声中的 0.2 Hz 正弦波组成，并以每秒四次的采样率进行 20 秒采样。

dt = 0.25;
t = (0:dt:20-1)';

x = 5*sin(2*pi*0.2*t)+0.5*randn(size(t));

使用萨维茨基-戈雷方法估计正弦波的前三个导数。使用包含 25 个采样的帧和五阶多项式。将各列除以 dt 的幂以对导数进行正确缩放。

[b,g] = sgolay(5,25);

dx = zeros(length(x),4);
for p = 0:3
  dx(:,p+1) = conv(x, factorial(p)/(-dt)^p * g(:,p+1), 'same');
end

绘制原始信号、平滑后的序列和导数估计值。

plot(x,'.-')
hold on
plot(dx)
hold off

legend('x','x (smoothed)','x''','x''''', 'x''''''')
title('Savitzky-Golay Derivative Estimates')

Figure contains an axes object. The axes object with title Savitzky-Golay Derivative Estimates contains 5 objects of type line. These objects represent x, x (smoothed), x', x'', x'''.

输入参数

全部折叠

`m` — 多项式阶数
非负整数

多项式阶数，指定为非负整数。m 的值必须小于 fl。如果 m = fl – 1，则设计的滤波器不产生平滑。

`fl` — 帧长度
正奇数

帧长度，指定为正奇数。fl 的值必须大于 m。

`w` — 加权向量
`ones(fl,1)` (默认) | 正实数向量

加权向量，指定为正实数向量。加权向量与 fl 具有相同长度，并用于执行最小二乘最小化。

输出参量

全部折叠

`b` — 时变 FIR 滤波器系数
矩阵

时变 FIR 滤波器系数，指定为 fl×fl 矩阵。在平滑滤波器实现中（例如，sgolayfilt），最后 (fl-1)/2 行（每行是一个 FIR 滤波器）在启动瞬变期间应用于信号，前 (fl-1)/2 行在终止瞬变期间应用于信号。中心行在稳态下应用于信号。

`g` — 微分滤波器的矩阵
矩阵

微分滤波器的矩阵，指定为矩阵。g 的每列是阶数为 p-1 的导数的微分滤波器，其中 p 是列索引。给定长度为 fl 的信号 x，您可以根据 xp((fl+1)/2) = (factorial(p)) * g(:,p+1)' * x 求得其中间值的 p阶导数 xp 的估计值。

算法

萨维茨基-戈雷滤波器用于平滑频率跨度较大的含噪信号。与标准平均值 FIR 滤波器相比，萨维茨基-戈雷平滑滤波器往往会滤除更少的信号高频成分。然而，当噪声电平特别高时，它们在抑制噪声方面不太成功。

通常，滤波包括用移动窗中包含的信号值的某种组合替换信号的每个点，该窗以该点为中心，假设附近点测量的是几乎相同的基础值。例如，移动平均滤波器用周围数据点的局部平均值替换每个数据点。如果给定数据点左侧有 k 个点，右侧有 k 个点，总窗长度为 L = 2k + 1，则移动平均滤波器进行替换

$x_{s} \to {\hat{x}}_{s} = \frac{1}{L} \sum_{r = - k}^{k} x_{s + r} .$

萨维茨基-戈雷滤波器通过最小二乘拟合穿过窗中信号值的 n 阶多项式并将拟合多项式曲线的计算中心点作为新平滑数据点，泛化此思路。对于给定点 x_s，

$[\begin{matrix} x_{s - k} \\ ⋮ \\ x_{s - 1} \\ x_{s} \\ x_{s + 1} \\ ⋮ \\ x_{s + k} \end{matrix}] = [\begin{matrix} b_{0} + b_{1} (t_{s} - k Δ t) + b_{2} {(t_{s} - k Δ t)}^{2} + \dots + b_{n} {(t_{s} - k Δ t)}^{n} \\ ⋮ \\ b_{0} + b_{1} (t_{s} - 1 Δ t) + b_{2} {(t_{s} - 1 Δ t)}^{2} + \dots + b_{n} {(t_{s} - 1 Δ t)}^{n} \\ b_{0} + b_{1} (t_{s} - 0 Δ t) + b_{2} {(t_{s} - 0 Δ t)}^{2} + \dots + b_{n} {(t_{s} - 0 Δ t)}^{n} \\ b_{0} + b_{1} (t_{s} + 1 Δ t) + b_{2} {(t_{s} + 1 Δ t)}^{2} + \dots + b_{n} {(t_{s} + 1 Δ t)}^{n} \\ ⋮ \\ b_{0} + b_{1} (t_{s} + k Δ t) + b_{2} {(t_{s} + k Δ t)}^{2} + \dots + b_{n} {(t_{s} + k Δ t)}^{n} \end{matrix}] = [\begin{matrix} a_{0} + a_{1} (- k) + a_{2} {(- k)}^{2} + \dots + a_{n} {(- k)}^{n} \\ ⋮ \\ a_{0} + a_{1} (- 1) + a_{2} {(- 1)}^{2} + \dots + a_{n} {(- 1)}^{n} \\ a_{0} + a_{1} (0) + a_{2} {(0)}^{2} + \dots + a_{n} {(0)}^{n} \\ a_{0} + a_{1} (1) + a_{2} {(1)}^{2} + \dots + a_{n} {(1)}^{n} \\ ⋮ \\ a_{0} + a_{1} (k) + a_{2} {(k)}^{2} + \dots + a_{n} {(k)}^{n} \end{matrix}]$

或者，用矩阵表示，

$x = [\begin{matrix} 1 & - k & {(- k)}^{2} & \dots & {(- k)}^{n} \\ ⋮ & ⋮ & ⋮ & ⋰ & ⋮ \\ 1 & - 2 & {(- 2)}^{2} & \dots & {(- 2)}^{n} \\ 1 & - 1 & {(- 1)}^{2} & \dots & {(- 1)}^{n} \\ 1 & 0 & 0 & \dots & 0 \\ 1 & 1 & 1^{2} & \dots & 1^{n} \\ 1 & 2 & 2^{2} & \dots & 2^{n} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & k & k^{2} & \dots & k^{n} \end{matrix}] [\begin{matrix} a_{0} \\ ⋮ \\ a_{n} \end{matrix}] \equiv H a .$

要求得萨维茨基-戈雷估计值，请使用 H 的伪逆计算 a，然后前乘 H：

$\hat{x} = H {(H^{T} H)}^{- 1} H^{T} x = B x .$

为避免病态，sgolay 使用 qr 函数将 H 的精简分解计算为 H = QR，据此 B = QQ^T。

只需计算一次 B。大多数信号点的萨维茨基-戈雷估计值来自将信号与 B 的中心行进行卷积。这将得到滤波信号的稳态部分。B 的前 k 行产生初始瞬变，B 的最后 k 行产生最终瞬变。有关示例，请参阅 sgolayfilt。可以通过增大窗长度来改进噪声抑制，但这会在任何瞬变附近引入类似于吉布斯现象的振铃效应。

参考

[1] Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1996.

[2] Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. New York: Cambridge University Press, 1992.

[3] Schafer, Ronald W. “What Is a Savitzky-Golay Filter? [Lecture Notes].” IEEE Signal Processing Magazine Vol. 28, Number 4, July 2011, pp. 111–117. https://doi.org/10.1109/MSP.2011.941097.

扩展功能

全部展开

C/C++ 代码生成
使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成
使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

sgolay 函数完全支持基于线程的环境。有关详细信息，请参阅在基于线程的环境中运行 MATLAB 函数。

版本历史记录

在 R2006a 之前推出

全部展开

R2026b: 使用基于线程的环境

您可以在基于线程的环境中使用 sgolay。有关详细信息，请参阅在基于线程的环境中运行 MATLAB 函数。

sgolay

语法

描述

示例

含噪正弦波的萨维茨基-戈雷平滑

带权重的萨维茨基-戈雷滤波器

萨维茨基-戈雷微分

输入参数

`m` — 多项式阶数
非负整数

`fl` — 帧长度
正奇数

`w` — 加权向量
`ones(fl,1)` (默认) | 正实数向量

输出参量

`b` — 时变 FIR 滤波器系数
矩阵

`g` — 微分滤波器的矩阵
矩阵

算法

参考

扩展功能

C/C++ 代码生成
使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成
使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

版本历史记录

R2026b: 使用基于线程的环境

另请参阅

函数

主题

sgolay

语法

描述

示例

含噪正弦波的萨维茨基-戈雷平滑

带权重的萨维茨基-戈雷滤波器

萨维茨基-戈雷微分

输入参数

m — 多项式阶数 非负整数

fl — 帧长度 正奇数

w — 加权向量 ones(fl,1) (默认) | 正实数向量

输出参量

b — 时变 FIR 滤波器系数 矩阵

g — 微分滤波器的矩阵 矩阵

算法

参考

扩展功能

C/C++ 代码生成 使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成 使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

基于线程的环境 使用 MATLAB® backgroundPool 在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool 加快代码运行速度。

版本历史记录

R2026b: 使用基于线程的环境

另请参阅

函数

主题

`m` — 多项式阶数
非负整数

`fl` — 帧长度
正奇数

`w` — 加权向量
`ones(fl,1)` (默认) | 正实数向量

`b` — 时变 FIR 滤波器系数
矩阵

`g` — 微分滤波器的矩阵
矩阵

C/C++ 代码生成
使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成
使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。