?1.Savitzky-Golay 濾波器

Savitzky-Golay濾波器（通常簡稱為S-G濾波器）最初由Savitzky和Golay于1964年提出，發表于Analytical Chemistry 雜志。之后被廣泛地運用于數據流平滑除噪，是一種在時域內基于局域多項式最小二乘法擬合的濾波方法。這種濾波器最大的特點在于在濾除噪聲的同時可以確保信號的形狀、寬度不變。
使用平滑濾波器對信號濾波時，實際上是擬合了信號中的低頻成分，而將高頻成分平滑出去了。如果噪聲在高頻端，那么濾波的結果就是去除了噪聲，反之，若噪聲在低頻段，那么濾波的結果就是留下了噪聲。 最成功的應用在于心電圖的基線漂移現象：
基線漂移現象是有很低的頻率信號引起的（體位的移動、呼吸作用），這樣的低頻信號有時我們也把它稱為信號中的趨勢項。可以利用平滑濾波器擬合該低頻信號，再將其從原始的心電信號中減去，就可以抑制趨勢項。 總之，平滑濾波是光譜分析中常用的預處理方法之一。用Savitzky．Golay方法進行平滑濾波，可以提高光譜的平滑性，并降低噪音的干擾。S-G平滑濾波的效果，隨著選取窗寬不同而不同，可以滿足多種不同場合的需求。

2.信號的最小二乘平滑

信號的最小二乘平滑的基本思想可以通過下圖來說明：

一列數據x[n] 在圖中用實心的圓點表示。考慮一組以n=0為中心的2M+1個數據，可以用如下的多項式來擬合：

所以只需要獲得擬合多項式的常數項。

Savitzky和Golay?發現計算a0相當于對原始數據進行一次FIR?濾波。也就是說可以利用卷積運算來實現。

也可以說是對輸入數據進行了加權平均。圖1中x?表示的就是加權系數。下面就來說說如何獲得a0。

由基本的微積分知識可知，若要ε最小，ε對各個參數的偏導數都應為0，也就是：

因此，Savitzky-Golay卷積平滑算法是移動平滑算法的改進：

3.matlab仿真測試

noise = wgn(1,101,2); for n=1:1:101s(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:98y(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濾波器的處理結果： 依據算法思想進行卷積濾波處理結果：

4.算法的C語言實現

要用C 語言從頭寫起這個代碼會非常的長，光一個 ployfit 函數實現起來就很麻煩。程序中直接調用了 GSL 的相關函數。先說說 polyfit 函數。多項式函數擬合可以轉換成超定線性代數方程組的最小二乘解的問題。標準解法是 SVD 分解。不過GSL庫中提供了gsl_multifit_linear函數進行線性模型的擬合，可以直接使用。下面給出代碼。很簡單，就不多做解釋了。這個例子也可以作為 gsl_multifit_linear 函數用法的一個小例子。 /** * 多項式擬合函數 * 根據 x y 擬合出一個N次多項式。返回多項式的系數。 */ 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 函數就可以計算 SG 濾波器的系數了。Matlab 中的 polyval 函數 可以用 gsl_poly_eval 來實現。
/** * 計算 Savitzky-Golay 濾波器系數 * SG 濾波器的階數為 2M+1，多項式的最高次數為 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.平滑濾波的尷尬

當然，這是針對所有FIR濾波器而言，都會有以下幾個尷尬，主要在模板值大小上： 1.N足夠大，就可以獲得足夠小的NRR（Noice Reduce Rate） 2.N過大會使濾波器具有過大的延遲 3.N過大會使其主瓣的單邊的帶寬大大降低，這就有可能在濾波時使有用的信號(或某一個頻率組份)也受到損失；因此，在平滑中，N不宜取得過大。 致謝：本文的大部分工作都是在CSDN?liyuanbhu博主研究基礎之上進行的，感謝他的辛苦勞動。

總結

以上是生活随笔為你收集整理的Savitzky-Golay 滤波器详解及C/matlab语言程序设计的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。