如何实现滑动平均滤波?

什么是滑动均值滤波

滑动平均滤波就是把连续取得的N个采样值看成一个队列,队列的长度固定为N,每次采样得到一个新数据放到队尾,并丢掉原来队首的一次数据,把队列中的N个数据进行平均运算,就可以获得新的滤波结果。


具体的matlab代码

clear

clc

load boxinfo.mat  %载入音频数据

T = data;

figure(1)

plot(T,'-*')

title('原始数据')

hold on;

%% 

%滑动平滑滤波

L = length(T);

N=10;  % 窗口大下

k = 0;

m =0 ;

for i = 1:L

m = m+1;

if i+N-1 > L

break

else

for j = i:N+i-1

k = k+1;

W(k) = T(j) ;

end

T1(m) = mean(W);

k = 0;

end

end

plot(T1,'r-o')

grid

legend('原始数据','滤波之后')


滤波前后对比图

滑动平均滤波


简单分析一下

经过滑动滤波之后,波形整体变得平滑,这里我们重点关注一下x轴附近的点,可以发现,在波形与x轴交叉的地方,波形都平稳过度,这极大方便的我们后期进行统计。


窗口大小选择

从代码中我们可以发现窗口大小我们选择的是10,如何选择窗口大小,这里我们需要进行一些简单的分析和测试。如果x轴附近的噪点数量(一上一下)比较多,那么窗口大小就应该大一些,反之,小一些。但是过大又会出现过拟合的现象,所以可以多取几个值,然后对比一下,选择一个最好的即可。


不同的窗口大小对比图

如何实现滑动平均滤波?


简单分析一下

从图中我们可以很明显的看出,当N=4的时候,滤波效果还不是很好,在x轴附近依然有噪点(一上一下),当N=7的时候,已经基本满足我们的要求,图形已经可以很平稳的过度了,但是从右边的标记处可以看出还是不是很平稳,所以可以继续提高N值,当N=10的时候,波形就完全能够达到我们的要求,所以取10即可。


滑动平均(moving average):在地球物理异常图上,选定某一尺寸的窗口,将窗口内的所有异常值做算术平均,将平均值作为窗口中心点的异常值。按点距或线距移动窗口,重复此平均方法,直到对整幅图完成上述过程,这种过程称为滑动平均。


滑动平均相当于低通滤波,在重力勘探和测井资料处理解释中常用此方法。 如果滑动窗长为n的话,滑动平均就是让数据通过一个n点的FIR滤波器,滤波器抽头系数都是1,这样取滑动平均就是起到序列平滑的作用。


利用filter函数求滑动平均

Matlab有多种计算滑动平均的方法,现介绍基于filter函数的计算方法。设原始数据为x,平均窗口设为a(a为正整数),那么无权重滑动平均后的数据y为:

windowSize =a;

y=filter(ones(1,windowSize)/windowSize,1,x);

上述命令实际上计算的是:

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);

... ....


可以看出,计算某一位置处的平均值时,窗口的前端位于该处。有时为了将窗口中部放在所计算的位置处,这样上述计算方式则变为(为叙述方便起见,设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的计算语句为:

windowSize =a;

y1=filter(ones(1,a/2+1)/windowSize,1,x);

y2=filter(ones(1,a/2+1)/windowSize,1,fliplr(x));

y=y1+fliplr(y2)-(1/a)*x;


如利用1-2-1 滤波器计算有权重的中心滑动平均,其Matlab语句为:

y1=filter([0.50.25],1,x);

y2=filter([0.5 0.25],1,fliplr(x));

y=y1+fliplr(y2)-0.5*x;


图片加载中...

在线留言

◎欢迎您的留言,您也可以通过以下方式联系我们:

◎客户服务热线:021-51095123

◎邮箱:xin021@126.com

展开