空间平滑滤波matlab代码,[整理]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