Design a Savitzky-Golay polynomial smoothing filter.
Compute a Savitzky-Golay polynomial smoothing filter. The resulting filter is placed in matrix H. Only the center row of H is used for filtering, the upper and lower part of the H matrix are applied to the transition region when the signal starts and stops. Diff contains differentiation filters. Weights contain weights for the least square minization. Order must be less then FrameSize and FrameSize must be Odd.
The resulting H vector contains FIR type impulse response, which can be passed to the FirInit routine.
Example 1
And excerpt from the Savitzky-Golay demo for a single block filter:
uses MtxExpr, Math387, MtxVec, SignalUtils, MtxVecTee, MtxVecEdit;
procedure TForm1.Button1Click(Sender: TObject);
var H,D: Matrix;
Data, Data1: Vector;
Frs,Ord: integer;
begin
Data.Size(100);
Data.RandGauss;
Data1.Copy(Data);
Frs := 10;
if not Odd(Frs) then Inc(Frs);
Ord := 5;
SavGolayImpulse(H,D,Frs,Ord,nil); //single block processing
SavGolayFilter(Data,H);
DrawIt([Data,Data1],['Filtered data','Original data']);
end;
#include "MtxVecCPP.h" //MtxVecCPP.cpp must be included in the project
#include "MtxVecEdit.hpp"
#include "MtxVecTee.hpp"
#include "SignalUtils.hpp"
void __fastcall TForm1::BitBtn1Click(TObject *Sender)
{
Matrix H,D;
Vector Data, Data1;
Data->Size(100);
Data->RandGauss();
Data1->Copy(Data);
int Frs = 10;
if ((Frs%2) == 0) Frs++;
int Ord = 5;
SavGolayImpulse(H,D,Frs,Ord,NULL); //single block processing
SavGolayFilter(Data,H);
DrawIt(OPENARRAY(TVec*,(Data,Data1)),OPENARRAY(AnsiString,("Filtered data","Original data")));
}
using Dew.Math;
using Dew.Math.Editors;
using Dew.Math.Units;
using Dew.Signal;
using Dew.Signal.Units;
using Dew.Math.Tee;
using Dew.Signal.Tee;
private void button1_Click(object sender, EventArgs e)
{
Matrix h = new Matrix(0,0);
Matrix d = new Matrix(0,0);
Vector data = new Vector(0);
Vector data1 = new Vector(0);
int frs = 10; //try changing frs here (must be odd)
int ord = 5;
data.Size(100);
data.RandGauss();
data1.Copy(data);
if (frs % 2 == 0) frs++;
ord = 5;
SignalUtils.SavGolayImpulse(h,d,frs,ord,null); //single block processing
SignalUtils.SavGolayFilter(data,h);
TeeChart.DrawIt(new TVec[2] { data, data1 }, new string[2] { "Filtered data", "Original data" },"Savitzky Golay",false);
}
Streamed filtering. A vector with a sine signal is broken down in to smaller pieces and they are filtered one by one.
procedure TForm1.Button1Click(Sender: TObject);
var h,b,c: Vector;
n,i: integer;
State: TFirState;
begin
SavGolayImpulse(h,15,7);
FirInit(h,State);
try
b := Ramp(300,0,0.1);
c.Size(b);
n := 10;
for i := 0 to (b.Length div n) - 1 do
begin
b.SetSubRange(i*n,n);
c.SetSubRange(i*n,n);
FirFilter(b,c,State);
end;
b.SetFullRange;
c.SetFullRange;
DrawIt([b,c],['Filtered data','Original data']);
finally
FirFree(State);
end;
end;
void __fastcall TForm1::BitBtn1Click(TObject *Sender)
{
Vector h,b,c;
int n,i;
TFirState State;
SavGolayImpulse(h,15,7);
FirInit(h,State);
try
{
b = Ramp(300,0,0.1);
c->Size(b);
n = 10;
for (i = 0; i< (b->Length/n); i++)
{
FirFilter(b(i*n, i*n+n-1),c(i*n, i*n+n-1),State);
}
DrawIt(OPENARRAY(TVec*,(b,c)),OPENARRAY(AnsiString,("Filtered data","Original data")));
}
__finally
{
FirFree(State);
}
}
using Dew.Math;
using Dew.Math.Editors;
using Dew.Math.Units;
using Dew.Signal;
using Dew.Signal.Units;
using Dew.Math.Tee;
using Dew.Signal.Tee;
private void button1_Click(object sender, EventArgs e)
{
Vector h = new Vector(0);
Vector b = MtxExpr.Ramp(300, 0, 0.1);
Vector c = new Vector(b.Length);
TFirState state = new TFirState();
int n = 10;
int i;
SignalUtils.SavGolayImpulse(h,15,7,null);
SignalUtils.FirInit(h,ref state,1,0,1,0);
try
{
c.Size(b);
n = 10;
int bLength = b.Length; //must be outside of the "for" to prevent reevaluation
for (i = 0; i< (bLength/n); i++)
{
b.SetSubRange(i * n, n); //select only a small subvector of the vector
c.SetSubRange(i * n, n);
SignalUtils.FirFilter(b,c,ref state);
}
TeeChart.DrawIt(new TVec[2] { b, c }, new string[2] { "Filtered data", "Original data" }, "Savitzky Golay", false);
}
finally
{
SignalUtils.FirFree(ref state);
}
}