轉(zhuǎn)載http://blog.csdn.net/liyuanbhu/article/details/11119081 以五點三次平滑為例。取相鄰的5個數(shù)據(jù)點，可以擬合出一條3次曲線來，然后用3次曲線上相應(yīng)的位置的數(shù)據(jù)值作為濾波后結(jié)果。簡單的說就是 Savitzky-Golay 濾波器 。只不過Savitzky-Golay 濾波器并不特殊考慮邊界的幾個數(shù)據(jù)點，而這個算法還特意把邊上的幾個點的數(shù)據(jù)擬合結(jié)果給推導(dǎo)了出來。 首先是線性擬合平滑處理的代碼. 分別為三點線性平滑、五點線性平滑和七點線性平滑。 void linearSmooth3 ( double in[], double out[], int N ){ int i; if ( N < 3 ) { for ( i = 0; i <= N - 1; i++ ) { out[i] = in[i]; } } else { out[0] = ( 5.0 * in[0] + 2.0 * in[1] - in[2] ) / 6.0; for ( i = 1; i <= N - 2; i++ ) { out[i] = ( in[i - 1] + in[i] + in[i + 1] ) / 3.0; } out[N - 1] = ( 5.0 * in[N - 1] + 2.0 * in[N - 2] - in[N - 3] ) / 6.0; }}void linearSmooth5 ( double in[], double out[], int N ){ int i; if ( N < 5 ) { for ( i = 0; i <= N - 1; i++ ) { out[i] = in[i]; } } else { out[0] = ( 3.0 * in[0] + 2.0 * in[1] + in[2] - in[4] ) / 5.0; out[1] = ( 4.0 * in[0] + 3.0 * in[1] + 2 * in[2] + in[3] ) / 10.0; for ( i = 2; i <= N - 3; i++ ) { out[i] = ( in[i - 2] + in[i - 1] + in[i] + in[i + 1] + in[i + 2] ) / 5.0; } out[N - 2] = ( 4.0 * in[N - 1] + 3.0 * in[N - 2] + 2 * in[N - 3] + in[N - 4] ) / 10.0; out[N - 1] = ( 3.0 * in[N - 1] + 2.0 * in[N - 2] + in[N - 3] - in[N - 5] ) / 5.0; }}void linearSmooth7 ( double in[], double out[], int N ){ int i; if ( N < 7 ) { for ( i = 0; i <= N - 1; i++ ) { out[i] = in[i]; } } else { out[0] = ( 13.0 * in[0] + 10.0 * in[1] + 7.0 * in[2] + 4.0 * in[3] + in[4] - 2.0 * in[5] - 5.0 * in[6] ) / 28.0; out[1] = ( 5.0 * in[0] + 4.0 * in[1] + 3 * in[2] + 2 * in[3] + in[4] - in[6] ) / 14.0; out[2] = ( 7.0 * in[0] + 6.0 * in [1] + 5.0 * in[2] + 4.0 * in[3] + 3.0 * in[4] + 2.0 * in[5] + in[6] ) / 28.0; for ( i = 3; i <= N - 4; i++ ) { out[i] = ( in[i - 3] + in[i - 2] + in[i - 1] + in[i] + in[i + 1] + in[i + 2] + in[i + 3] ) / 7.0; } out[N - 3] = ( 7.0 * in[N - 1] + 6.0 * in [N - 2] + 5.0 * in[N - 3] + 4.0 * in[N - 4] + 3.0 * in[N - 5] + 2.0 * in[N - 6] + in[N - 7] ) / 28.0; out[N - 2] = ( 5.0 * in[N - 1] + 4.0 * in[N - 2] + 3.0 * in[N - 3] + 2.0 * in[N - 4] + in[N - 5] - in[N - 7] ) / 14.0; out[N - 1] = ( 13.0 * in[N - 1] + 10.0 * in[N - 2] + 7.0 * in[N - 3] + 4 * in[N - 4] + in[N - 5] - 2 * in[N - 6] - 5 * in[N - 7] ) / 28.0; }}

然后是利用二次函數(shù)擬合平滑。

void quadraticSmooth5(double in[], double out[], int N){ int i; if ( N < 5 ) { for ( i = 0; i <= N - 1; i++ ) { out[i] = in[i]; } } else { out[0] = ( 31.0 * in[0] + 9.0 * in[1] - 3.0 * in[2] - 5.0 * in[3] + 3.0 * in[4] ) / 35.0; out[1] = ( 9.0 * in[0] + 13.0 * in[1] + 12 * in[2] + 6.0 * in[3] - 5.0 *in[4]) / 35.0; for ( i = 2; i <= N - 3; i++ ) { out[i] = ( - 3.0 * (in[i - 2] + in[i + 2]) + 12.0 * (in[i - 1] + in[i + 1]) + 17 * in[i] ) / 35.0; } out[N - 2] = ( 9.0 * in[N - 1] + 13.0 * in[N - 2] + 12.0 * in[N - 3] + 6.0 * in[N - 4] - 5.0 * in[N - 5] ) / 35.0; out[N - 1] = ( 31.0 * in[N - 1] + 9.0 * in[N - 2] - 3.0 * in[N - 3] - 5.0 * in[N - 4] + 3.0 * in[N - 5]) / 35.0; }}void quadraticSmooth7(double in[], double out[], int N){ int i; if ( N < 7 ) { for ( i = 0; i <= N - 1; i++ ) { out[i] = in[i]; } } else { out[0] = ( 32.0 * in[0] + 15.0 * in[1] + 3.0 * in[2] - 4.0 * in[3] - 6.0 * in[4] - 3.0 * in[5] + 5.0 * in[6] ) / 42.0; out[1] = ( 5.0 * in[0] + 4.0 * in[1] + 3.0 * in[2] + 2.0 * in[3] + in[4] - in[6] ) / 14.0; out[2] = ( 1.0 * in[0] + 3.0 * in [1] + 4.0 * in[2] + 4.0 * in[3] + 3.0 * in[4] + 1.0 * in[5] - 2.0 * in[6] ) / 14.0; for ( i = 3; i <= N - 4; i++ ) { out[i] = ( -2.0 * (in[i - 3] + in[i + 3]) + 3.0 * (in[i - 2] + in[i + 2]) + 6.0 * (in[i - 1] + in[i + 1]) + 7.0 * in[i] ) / 21.0; } out[N - 3] = ( 1.0 * in[N - 1] + 3.0 * in [N - 2] + 4.0 * in[N - 3] + 4.0 * in[N - 4] + 3.0 * in[N - 5] + 1.0 * in[N - 6] - 2.0 * in[N - 7] ) / 14.0; out[N - 2] = ( 5.0 * in[N - 1] + 4.0 * in[N - 2] + 3.0 * in[N - 3] + 2.0 * in[N - 4] + in[N - 5] - in[N - 7] ) / 14.0; out[N - 1] = ( 32.0 * in[N - 1] + 15.0 * in[N - 2] + 3.0 * in[N - 3] - 4.0 * in[N - 4] - 6.0 * in[N - 5] - 3.0 * in[N - 6] + 5.0 * in[N - 7] ) / 42.0; }} 最后是三次函數(shù)擬合平滑。

/** * 五點三次平滑 * */void cubicSmooth5 ( double in[], double out[], int N ){ int i; if ( N < 5 ) { for ( i = 0; i <= N - 1; i++ ) out[i] = in[i]; } else { out[0] = (69.0 * in[0] + 4.0 * in[1] - 6.0 * in[2] + 4.0 * in[3] - in[4]) / 70.0; out[1] = (2.0 * in[0] + 27.0 * in[1] + 12.0 * in[2] - 8.0 * in[3] + 2.0 * in[4]) / 35.0; for ( i = 2; i <= N - 3; i++ ) { out[i] = (-3.0 * (in[i - 2] + in[i + 2])+ 12.0 * (in[i - 1] + in[i + 1]) + 17.0 * in[i] ) / 35.0; } out[N - 2] = (2.0 * in[N - 5] - 8.0 * in[N - 4] + 12.0 * in[N - 3] + 27.0 * in[N - 2] + 2.0 * in[N - 1]) / 35.0; out[N - 1] = (- in[N - 5] + 4.0 * in[N - 4] - 6.0 * in[N - 3] + 4.0 * in[N - 2] + 69.0 * in[N - 1]) / 70.0; } return;}void cubicSmooth7(double in[], double out[], int N){ int i; if ( N < 7 ) { for ( i = 0; i <= N - 1; i++ ) { out[i] = in[i]; } } else { out[0] = ( 39.0 * in[0] + 8.0 * in[1] - 4.0 * in[2] - 4.0 * in[3] + 1.0 * in[4] + 4.0 * in[5] - 2.0 * in[6] ) / 42.0; out[1] = ( 8.0 * in[0] + 19.0 * in[1] + 16.0 * in[2] + 6.0 * in[3] - 4.0 * in[4] - 7.0* in[5] + 4.0 * in[6] ) / 42.0; out[2] = ( -4.0 * in[0] + 16.0 * in [1] + 19.0 * in[2] + 12.0 * in[3] + 2.0 * in[4] - 4.0 * in[5] + 1.0 * in[6] ) / 42.0; for ( i = 3; i <= N - 4; i++ ) { out[i] = ( -2.0 * (in[i - 3] + in[i + 3]) + 3.0 * (in[i - 2] + in[i + 2]) + 6.0 * (in[i - 1] + in[i + 1]) + 7.0 * in[i] ) / 21.0; } out[N - 3] = ( -4.0 * in[N - 1] + 16.0 * in [N - 2] + 19.0 * in[N - 3] + 12.0 * in[N - 4] + 2.0 * in[N - 5] - 4.0 * in[N - 6] + 1.0 * in[N - 7] ) / 42.0; out[N - 2] = ( 8.0 * in[N - 1] + 19.0 * in[N - 2] + 16.0 * in[N - 3] + 6.0 * in[N - 4] - 4.0 * in[N - 5] - 7.0 * in[N - 6] + 4.0 * in[N - 7] ) / 42.0; out[N - 1] = ( 39.0 * in[N - 1] + 8.0 * in[N - 2] - 4.0 * in[N - 3] - 4.0 * in[N - 4] + 1.0 * in[N - 5] + 4.0 * in[N - 6] - 2.0 * in[N - 7] ) / 42.0; }}

Matlab代碼

加高斯白噪

function [Y,NOISE] = noisegen(X,SNR)
% noisegen add white Gaussian noise to a signal.
% [Y, NOISE] = NOISEGEN(X,SNR) adds white Gaussian NOISE to X. The SNR is in dB.
%X是純信號，SNR是要求的信噪比，Y是帶噪信號，NOISE是疊加在信號上的噪聲。
NOISE=randn(size(X));
NOISE=NOISE-mean(NOISE);
signal_power = 1/length(X)*sum(X.*X);
noise_variance = signal_power / ( 10^(SNR/10) );
NOISE=sqrt(noise_variance)/std(NOISE)*NOISE;
Y=X+NOISE;

五點三次平滑濾波

function y = smooth5_3(x, m)
% x為需要的數(shù)據(jù)
% m 為平滑次數(shù)
n=length(x);
a=x;
for k=1: m
b(1) = (69*a(1) +4*(a(2) +a(4)) -6*a(3) -a(5)) /70;
b(2) = (2* (a(1) +a(5)) +27*a(2) +12*a(3) -8*a(4)) /35;
for j=3:n-2
b (j) = (-3*(a(j-2) +a(j+2)) +12*(a(j-1) +a(j+1)) +17*a(j)) /35;
end
b (n-1) = (2*(a(n) +a(n-4)) +27*a(n-1) +12*a(n-2) -8*a(n-3)) /35;
b (n) = (69*a(n) +4* (a(n-1) +a(n-3)) -6*a(n-2) -a(n-4)) /70;
a=b;
end
y =a;

總結(jié)

以上是生活随笔為你收集整理的五点三次平滑法滤波 C 和 matlab代码的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。