2015年11月9日 星期一

[Algorithm] 移動平均濾波 Moving Average

平均法是濾波上最容易實現的一種方法之一,常常被拿來應用,像是每 10 筆資料做一次平均值之類的,但這種方式有下面幾個缺點:

1、假設每秒取樣 n 筆資料時,每 m 筆資料做平均值,更新速度會等於 n/m Hz,低於取樣速度。
2、資料曲線不夠平滑,每筆平均數間容易產生落差。


移動平均顧名思義就是移動 + 平均,此方法可以解決上述問題,在不影響更新速度的同時增加平滑度,運算量也不大,下面將介紹簡單的移動平均與加權的移動平均兩種方法。


簡單移動平均 Simple Moving Average, SMA

移動平均顧名思義就是移動 + 平均,在固定的"視窗大小"下取平均值,而這個視窗會移動,其實就是摺積(Convolution)

gif 圖片截自 wikipedia - Convolution https://en.wikipedia.org/wiki/Convolution


下面是一個實際的範例,紅色框框表示視窗,程式上可以理解成 FIFO,大小為 4 筆資料,初值為 0,總共有 8 筆資料,有問號的 8 筆資料是移動平均算出來的資料,一開始因為 FIFO 的資料沒有被填滿,所以誤差比較大。



大致的概念就是現一個長度 N 的 FIFO,最新的資料進去,最舊的資料出來,FIFO 裡面始終保持 N 筆最新的資料,對裡面 N 筆資料做平均值就是移動平均了,N 的大小會影響到訊號的響應速度,N 越大,越容易抗雜訊,但響應速度會變慢,N 越小,越不易抗雜訊,但響應速度越快,當 N = 1 時,其實就是實際的資料。

下面是自己用 MATLAB 撰寫的移動平均,之前也有寫過一個,不過是直接讓 FIFO 的資料移動,也就是頭跟尾始終保持最新和最舊的資料,但在這種方法下,更新 FIFO 的部分運算量會隨著 FIFO 的長度增加而增加,所以就另外寫了一個直接將最舊的資料替代成最新的資料,更新 FIFO 的運算量在不同長度下將不會變,效率應該也高於前者。


簡單移動平均 moveAve_s.m
function [aveData, p_front] = moveAve_s( newData, moveAveFIFO, fifoLens, windowLens )

    p_front = moveAveFIFO(1);
    p_rear  = mod(p_front - windowLens + fifoLens, fifoLens) + 1;

    sumData = 0;
    if p_front > p_rear
        for i = p_front : -1 : p_rear
            sumData = sumData + moveAveFIFO(i + 1);
        end
    else
        for i = p_front : -1 : 1
            sumData = sumData + moveAveFIFO(i + 1);
        end
        for i = fifoLens : -1 : p_rear
            sumData = sumData + moveAveFIFO(i + 1);
        end
    end

    aveData = sumData / windowLens;
    p_front = mod(p_front, fifoLens) + 1;

end


加權移動平均 Weighted Moving Average, WMA

加權移動平均與簡單的移動平均只差在每次做累加的時候會乘上權重,透過權重來調整資料的重要性或影響性。


加權移動平均 moveAve_w.m
function [aveData, p_front] = moveAve_w( newData, moveAveFIFO, weighted, fifoLens, windowLens )

    p_front = moveAveFIFO(1);
    p_rear  = mod(p_front - windowLens + fifoLens, fifoLens) + 1;

    sumData = 0;
    count = 1;
    if p_front > p_rear
        for i = p_front : -1 : p_rear
            sumData = sumData + moveAveFIFO(i + 1) * weighted(count + 1);
            count = count + 1;
        end
    else
        for i = p_front : -1 : 1
            sumData = sumData + moveAveFIFO(i + 1) * weighted(count + 1);
            count = count + 1;
        end
        for i = fifoLens : -1 : p_rear
            sumData = sumData + moveAveFIFO(i + 1) * weighted(count + 1);
            count = count + 1;
        end
    end

    aveData = sumData / weighted(1);
    p_front = mod(p_front, fifoLens) + 1;

end


驗證兩種移動平均方法 testMoveAve.m
clear all;

timeScal = 0.01;
runTimes = 10;

fifoLens = 256;
windowLens = 32;
fifo = zeros(1, fifoLens+1);
fifo(1) = windowLens;

weighted = zeros(1, windowLens + 1);
for i = 1 : windowLens
    weighted(i + 1) = 2^((windowLens - i) / 4);
end
weighted(1) = sum(weighted(2 : end));

saveLens = 0;
for t = 0 : timeScal : runTimes
    newData = sin(t) + randn();

    p_front = fifo(1);
    fifo(p_front + 1) = newData;
    [sma_data, p_front] = moveAve_s(newData, fifo, fifoLens, windowLens);
    [wma_data, p_front] = moveAve_w(newData, fifo, weighted, fifoLens, windowLens);
    fifo(1) = p_front;

    saveLens = saveLens + 1;
    save_org(saveLens)  = sin(t);
    save_orgn(saveLens) = newData;
    save_sma(saveLens)  = sma_data;
    save_wma(saveLens)  = wma_data;
end

hold on;
grid on;
plot([0 : timeScal : runTimes], save_org,  'g');
plot([0 : timeScal : runTimes], save_orgn, 'r');
plot([0 : timeScal : runTimes], save_sma,  'b');
plot([0 : timeScal : runTimes], save_wma,  'c');

xlabel('t(s)');
ylabel('data');
legend('org', 'org + noise', 'sma', 'wma');


Github→https://github.com/Hom-Wang/MATLAB/tree/master/moveAve