轉(zhuǎn)自 http://blog.csdn.net/shenziheng1/article/details/53391422
1.Savitzky-Golay 濾波器
Savitzky-Golay濾波器（通常簡稱為S-G濾波器）最初由Savitzky和Golay于1964年提出征唬，發(fā)表于Analytical Chemistry 雜志。之后被廣泛地運用于數(shù)據(jù)流平滑除噪茁彭，是一種在時域內(nèi)基于局域多項式最小二乘法擬合的濾波方法总寒。這種濾波器最大的特點在于在濾除噪聲的同時可以確保信號的形狀、寬度不變理肺。
使用平滑濾波器對信號濾波時摄闸，實際上是擬合了信號中的低頻成分善镰，而將高頻成分平滑出去了。如果噪聲在高頻端年枕，那么濾波的結(jié)果就是去除了噪聲炫欺，反之，若噪聲在低頻段熏兄，那么濾波的結(jié)果就是留下了噪聲品洛。
最成功的應(yīng)用在于心電圖的基線漂移現(xiàn)象：
image基線漂移現(xiàn)象是有很低的頻率信號引起的（體位的移動、呼吸作用）摩桶，這樣的低頻信號有時我們也把它稱為信號中的趨勢項桥状。可以利用平滑濾波器擬合該低頻信號硝清，再將其從原始的心電信號中減去辅斟，就可以抑制趨勢項。
總之芦拿，平滑濾波是光譜分析中常用的預(yù)處理方法之一士飒。用Savitzky．Golay方法進(jìn)行平滑濾波，可以提高光譜的平滑性防嗡，并降低噪音的干擾变汪。S-G平滑濾波的效果，隨著選取窗寬不同而不同蚁趁，可以滿足多種不同場合的需求。
2.信號的最小二乘平滑
信號的最小二乘平滑的基本思想可以通過下圖來說明：
一列數(shù)據(jù)x[n] 在圖中用實心的圓點表示实胸∷眨考慮一組以n=0為中心的2M+1個數(shù)據(jù)，可以用如下的多項式來擬合：imageimage所以只需要獲得擬合多項式的常數(shù)項庐完。
Savitzky和Golay 發(fā)現(xiàn)計算a0相當(dāng)于對原始數(shù)據(jù)進(jìn)行一次FIR 濾波钢属。也就是說可以利用卷積運算來實現(xiàn)。
image也可以說是對輸入數(shù)據(jù)進(jìn)行了加權(quán)平均门躯。圖1中x 表示的就是加權(quán)系數(shù)淆党。下面就來說說如何獲得a0。
由基本的微積分知識可知讶凉，若要ε最小染乌，ε對各個參數(shù)的偏導(dǎo)數(shù)都應(yīng)為0****，也就是：
image因此懂讯，Savitzky-Golay卷積平滑算法是移動平滑算法的改進(jìn)：
image
3.matlab仿真測試
noise = wgn(1,101,2);
for n=1:1:101
s(1,n) = 3*sin(0.2*n)+3*sin(0.4*n);
end
x = s + noise;
%濾波器
% y = sgolayfilt(x,3,7);
h = [-2, 3, 6, 7, 6, 3, -2 ]/21;
y = zeros(1,101);
y(1,1:3) = s(1,1:3);
y(1,99:101) = s(1,99:101);
for k=4:1:98
y(1,k) = x(1,k-3)*h(1)+x(1,k-2)*h(2)+x(1,k-1)*h(3)+x(1,k)*h(4)+...
x(1,k+1)*h(5)+x(1,k+2)*h(6)+x(1,k+3)*h(7);
end
plot(x,'b');
hold on
plot(y,'r');
hold on
plot(s,'g');
hold off
MATLAB自帶的sgolayfilt濾波器的處理結(jié)果：
圖片.png
依據(jù)算法思想進(jìn)行卷積濾波處理結(jié)果：
4.算法的C語言實現(xiàn)
要用C 語言從頭寫起這個代碼會非常的長荷憋，光一個 ployfit 函數(shù)實現(xiàn)起來就很麻煩。程序中直接調(diào)用了 GSL 的相關(guān)函數(shù)褐望。先說說 polyfit 函數(shù)勒庄。多項式函數(shù)擬合可以轉(zhuǎn)換成超定線性代數(shù)方程組的最小二乘解的問題串前。標(biāo)準(zhǔn)解法是 SVD 分解。不過GSL庫中提供了gsl_multifit_linear函數(shù)進(jìn)行線性模型的擬合实蔽，可以直接使用荡碾。下面給出代碼。很簡單局装，就不多做解釋了坛吁。這個例子也可以作為 gsl_multifit_linear 函數(shù)用法的一個小例子。
/**
* 多項式擬合函數(shù)
* 根據(jù) x y 擬合出一個N次多項式贼邓。返回多項式的系數(shù)阶冈。
*/
gsl_vector* PolyFit(gsl_vector *x, gsl_vector *y, int N)
{
int i, j;
int M = GSL_MIN(x->size, y->size);
double chisq, xi;
gsl_matrix *XX = gsl_matrix_alloc(M, N + 1);
gsl_vector *c = gsl_vector_alloc(N + 1);
gsl_matrix *cov = gsl_matrix_alloc(N + 1, N + 1);
for(i = 0; i < M; i++)
{
gsl_matrix_set(XX, i, 0, 1.0);
xi = gsl_vector_get(x, i);
for(j = 1; j <= N; j++)
{
gsl_matrix_set(XX, i, j, pow(xi, j));
}
}
gsl_multifit_linear_workspace *workspace = gsl_multifit_linear_alloc(M, N + 1);
gsl_multifit_linear(XX, y, c, cov, &chisq, workspace);
gsl_matrix_free(XX);
gsl_matrix_free(cov);
gsl_multifit_linear_free(workspace);
return c;
}
有了 PolyFit 函數(shù)就可以計算 SG 濾波器的系數(shù)了。Matlab 中的 polyval 函數(shù) 可以用 gsl_poly_eval 來實現(xiàn)塑径。
/**
* 計算 Savitzky-Golay 濾波器系數(shù)
* SG 濾波器的階數(shù)為 2M+1女坑，多項式的最高次數(shù)為 N
*/
gsl_vector* SG_FilterCreate(int M, int N /* Poly Order */)
{
int i;
gsl_vector *x = gsl_vector_alloc(2 * M + 1);
gsl_vector *y = gsl_vector_alloc(2 * M + 1);
gsl_vector_set_zero(y);
gsl_vector_set(y, M, 1);
for(i= -M; i <= M; i++)
{
gsl_vector_set(x, i + M, i);
}
gsl_vector *c = PolyFit(x, y, N);
gsl_vector_free(x);
gsl_vector_free(y);
gsl_vector *fir = gsl_vector_alloc(2 * M + 1);
for(i = -M; i <= M; i++)
{
gsl_vector_set(fir, i + M, gsl_poly_eval (c->data, N + 1, i));
}
gsl_vector_free(c);
return fir;
}
5.平滑濾波的尷尬
當(dāng)然，這是針對所有FIR濾波器而言统舀，都會有以下幾個尷尬匆骗，主要在模板值大小上：
1.N足夠大，就可以獲得足夠小的NRR（Noice Reduce Rate）
2.N過大會使濾波器具有過大的延遲
3.N過大會使其主瓣的單邊的帶寬大大降低誉简，這就有可能在濾波時使有用的信號(或某一個頻率組份)也受到損失碉就；因此，在平滑中闷串，N不宜取得過大瓮钥。
致謝：本文的大部分工作都是在CSDN liyuanbhu博主研究基礎(chǔ)之上進(jìn)行的，感謝他的辛苦勞動烹吵。
