[整理]Matlab 之中心平滑濾波

滑動平均(moving average)：在地球物理異常圖上，選定某一尺寸的窗口，將窗口內的所有異常值做算術平均，將平均值作為窗口中心點的異常值。按點距或線距移動窗口，重復此平均方法，直到對整幅圖完成上述過程，這種過程稱為滑動平均。

滑動平均相當於低通濾波，在重力勘探和測井資料處理解釋中常用此方法。

如果滑動窗長為 n 的話，滑動平均就是讓數據通過一個 n 點的 FIR 濾波器，濾波器抽頭系數都是 1，這樣取滑動平均就是起到序列平滑的作用。

Matlab 中有多種滑動平均方法，比如 filter 和 tsmovavg 方法都可以實現。

普通滑動平均

基於filter的普通無權重滑動平均，有關於 filter 函數，可以 doc filter進行詳細信息的查看，這里我們由於只使用了簡單的滑動平均，在此記錄一種簡單的滑動平均方法。

seqOriginal = rand(1,100);
windowSize  = x;
seqFilter   = filter( ones(1, windowSize) / windowSize, 1, seqOriginal );

上述命令實際計算的是：

x表示seqOriginal, y表示seqFilter, a表示windowSize。
y(1) = (1 / a) * x(1);
y(2) = (1 / a) * x(2) + (1 / a) * x(1);
...
y(a) = (1 / a) * x(a) + (1 / a) * x(a - 1) + ... + (1 / a) * x(1);
...
y(i) = (1 / a) * x(i) + (1 / a) * x(i - 1) + ... + (1 / a) * x(i - a + 1);
...

注：該方法由於是計算該點和之前 windowSize 的點的平均值，故其輸出結果相對於原數據趨勢有一個滯后。如果數據量比較少，可能影響較大。

中心滑動平均

還有一種滑動平均的方法為中心滑動平均，數據點位於滑動窗的中心進行平均值的計算。則上述的計算變為（在此為了方便虛數，設a為奇數）:

y(1) = (1 / a) * x(1) + (1 / a) * x(2) + ... + (1 / a) * x((a+1) / 2);
y(2) = (1 / a) * x(1) + (1 / a) * x(2) + ... + (1 / a) * x((a+1) / 2 + 1);
...
y((a + 1) / 2) = (1 / a) * x(1) + (1 / a) * x(2) + ... + (1 / a) * x((a+1) / 2) + (1 / a) * x(a);
...
y(i) = (1 / a) * x(i - (a - 1) / 2) + (1 / a) * x(i - (a - 1) / 2 + 1) + ... + (1 / a) * x(i) + ... + (1 / a) * x(i + (a+1) / 2);
...

將上述式子轉換為 Matlab 語句為：

seqOriginal = rand(1,100);
windowSize  = x;
seq1        = filter( ones(1, windowSize / 2 + 1) / windowSize, 1, seqOriginal );
seq2        = filter( ones(1, windowSize / 2 + 1) / windowSize, 1, fliplr(seqOriginal) );
seqFilter   = seq1 + fliplr(seq2) - 1 / windowSize * seqOriginal;

為了方便使用，為其寫了一個 function 函數用來調用。

%fun_CenterAverageFilter
%--
%   seqFilter = fun_CenterAverageFilter(seqOriginal, windowSize)
%   中心滑動平均
%--
%   seqFilter   |Matrix|    濾波后輸出序列  
%   seqOriginal |Matrix|    原始序列
%   windowSize  |Double|    滑動窗口
%--
%Author:    Liu Tong
%History:
%----
%Rev:  1.0
%Date: 2016/12/22
%   create.
%--
function seqFilter = fun_CenterAverageFilter(seqOriginal, windowSize)
seq1        = filter( ones(1, windowSize / 2 + 1) / windowSize, 1, seqOriginal );
seq2        = filter( ones(1, windowSize / 2 + 1) / windowSize, 1, fliplr(seqOriginal) );
seqFilter   = seq1 + fliplr(seq2) - 1 / windowSize * seqOriginal;
end