このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

sgolay

Savitzky-Golay フィルター設計

構文

b = sgolay(order,framelen)

b = sgolay(order,framelen,weights)

[b,g] = sgolay(___)

説明

b = sgolay(order,framelen) では、多項式次数 order およびフレーム長 framelen を持つ Savitzky-Golay FIR 平滑化フィルターが設計されます。

b = sgolay(order,framelen,weights) では、最小二乗の最小化で使用される正の実数値の重みを含む、重みベクトル weights が指定されます。

[b,g] = sgolay(___) では、微分フィルターの行列 g が返されます。以前の任意の入力構文と共にこれらの出力引数を使用できます。

例

ノイズの多い正弦波の Savitzky-Golay 平滑化

ホワイト ガウス ノイズに含まれる 0.2 Hz の正弦波で構成される信号を生成し、200 秒間、毎秒 5 回サンプリングします。

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

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

sgolay を使用して信号を平滑化します。21 サンプルのフレームと 4 次多項式を使用します。

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));

過渡状態と定常状態部分を連結して、完全に平滑化された信号を生成します。元の信号および Savitzky-Golay 推定をプロットします。

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

Savitzky-Golay 微分

ホワイト ガウス ノイズに含まれる 0.2 Hz の正弦波で構成される信号を生成し、20 秒間、毎秒 4 回サンプリングします。

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

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

Savitzky-Golay 法を使用して、正弦波の最初の 3 つの導関数を推定します。25 サンプルのフレームと 5 次多項式を使用します。列を 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')

入力引数

order — 多項式の次数

正の整数

多項式の次数。正の整数として指定します。order の値は framelen より小さくなければなりません。order = framelen – 1 の場合、設計されるフィルターでは平滑化は行われません。

framelen — フレーム長

正の奇数の整数

フレーム長。正の奇数の整数で指定します。framelen の値は order より大きくなければなりません。

weights — 重みベクトル

正の実数ベクトル

重みベクトル。正の実数ベクトルとして指定します。重みベクトルの長さは framelen と同じで、最小二乗の最小化を実行するために使用されます。

出力引数

b — 時変 FIR フィルター係数

行列

時変 FIR フィルター係数。framelen 行 framelen 列の行列として指定します。平滑化フィルターを実装 (たとえば、sgolayfilt を実装) する際には、最後の (framelen-1)/2 行 (個々の FIR フィルター) は、開始時の過渡状態中に信号に適用されます。また、最初の (framelen-1)/2 行は、終了時の過渡状態中に信号に適用されます。中心の行は、安定状態で信号に適用されます。

g — 微分フィルターの行列

行列

微分フィルターの行列。行列として指定します。g の各列は、次数 p-1 の導関数に対応する差分フィルターです。ここで、p は列インデックスです。長さが framelen の信号 x が与えられた場合、xp((framelen+1)/2) = (factorial(p)) * g(:,p+1)' * x より、中央値である p 次の導関数 xp の推定を求めることができます。

アルゴリズム

Savitzky-Golay フィルターは、周波数範囲が広いノイズ信号の平滑化に使用されます。Savitzky-Golay 平滑化フィルターは、標準の平均化 FIR フィルターに比べ、信号の高周波数成分を維持しやすいという傾向があります。ただし、ノイズ レベルが著しく高い場合のノイズの除去には優れていません。

一般に、フィルター処理では、信号の各点を、その点を中心とする移動ウィンドウに含まれる信号値のいくつかの組み合わせで置き換えます。この処理は、近傍点はほぼ同じ基礎値を測定するという仮定に基づいています。たとえば、移動平均フィルターは、各データ点を周囲のデータ点の局所的な平均で置き換えます。与えられたデータ点が左に k 点、右に k 点をもち、ウィンドウの長さの合計が L = 2k + 1 である場合、移動平均フィルターは次の置き換えを行います。

$$x_s\to {\widehat{x}}_s=\frac{1}{L}{\displaystyle \sum _{r=-k}^k{x}_{s+r}.}$$

Savitzky-Golay フィルターは、ウィンドウ内の信号値で n 次多項式を最小二乗近似し、近似された多項式曲線について計算された中心点を新たに平滑化されたデータ点として取得することで、この理念を一般化します。与えられた点 x_s について、次のようになります。

$$\left[\begin{array}{c}{x}_{s-k}\\ ⋮\\ x_{s-1}\\ x_s\\ x_{s+1}\\ ⋮\\ x_{s+k}\end{array}\right]=\left[\begin{array}{c}{b}_0+b_1\left(t_s-k\Delta t\right)+b_2{\left(t_s-k\Delta t\right)}^2+\cdots +b_n{\left(t_s-k\Delta t\right)}^n\\ ⋮\\ b_0+b_1\left(t_s-1\Delta t\right)+b_2{\left(t_s-1\Delta t\right)}^2+\cdots +b_n{\left(t_s-1\Delta t\right)}^n\\ b_0+b_1\left(t_s-0\Delta t\right)+b_2{\left(t_s-0\Delta t\right)}^2+\cdots +b_n{\left(t_s-0\Delta t\right)}^n\\ b_0+b_1\left(t_s+1\Delta t\right)+b_2{\left(t_s+1\Delta t\right)}^2+\cdots +b_n{\left(t_s+1\Delta t\right)}^n\\ ⋮\\ b_0+b_1\left(t_s+k\Delta t\right)+b_2{\left(t_s+k\Delta t\right)}^2+\cdots +b_n{\left(t_s+k\Delta t\right)}^n\end{array}\right]=\left[\begin{array}{c}{a}_0+a_1\left(-k\right)+a_2{\left(-k\right)}^2+\cdots +a_n{\left(-k\right)}^n\\ ⋮\\ a_0+a_1\left(-1\right)+a_2{\left(-1\right)}^2+\cdots +a_n{\left(-1\right)}^n\\ a_0+a_1\left(0\right)+a_2{\left(0\right)}^2+\cdots +a_n{\left(0\right)}^n\\ a_0+a_1\left(1\right)+a_2{\left(1\right)}^2+\cdots +a_n{\left(1\right)}^n\\ ⋮\\ a_0+a_1\left(k\right)+a_2{\left(k\right)}^2+\cdots +a_n{\left(k\right)}^n\end{array}\right]$$

または、行列という観点では次のようになります。

$$x=\left[\begin{array}{ccccc}1& -k& {\left(-k\right)}^2& \cdots & {\left(-k\right)}^n\\ ⋮& ⋮& ⋮& ⋰& ⋮\\ 1& -2& {\left(-2\right)}^2& \cdots & {\left(-2\right)}^n\\ 1& -1& {\left(-1\right)}^2& \cdots & {\left(-1\right)}^n\\ 1& 0& 0& \cdots & 0\\ 1& 1& 1^2& \cdots & 1^n\\ 1& 2& 2^2& \cdots & 2^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& k& k^2& \cdots & k^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ ⋮\\ a_n\end{array}\right]\equiv Ha.$$

Savitzky-Golay 推定値を求めるには、H の疑似逆数を使用して a を計算してから、H でプリマルチプライ処理を行います。

$$\widehat{x}=H{\left(H^{T}H\right)}^{-1}{H}^{T}x=Bx.$$

悪条件を回避するために、sgolay では、関数 qr を使用して、B = QQ^T という観点で H の無駄のないサイズの分解を H = QR として計算します。

B は 1 回のみ計算する必要があります。大半の信号点についての Savitzky-Golay 推定値は、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++ コードを生成します。

バージョン履歴

R2006a より前に導入

参考

fir1 | firls | filter | sgolayfilt