什麼是滑動均值濾波
滑動平均濾波就是把連續取得的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;