RE: mean filter, another implementation

[andreas naessl]

hello, i experimented a bit with mean filter - implementations

i found my recmean_2 - function as fast as francescos solution, but the way is 
a bit different. see below. the conv-solution becomes a bit slow with long 
filter length... is there a performance difference between windows and linux 
versions ?

best regards, a. naessl


#> >>I'm looking for a *fast* implementation of a mean filter

0;

function m = meanfilter (x, ws)    #> Francesco Potorti`s solution
     l = length(x);
     a = cumsum(prepad(x, l+1));
     b = a(ws+1:l)-a(1:l-ws);
     m = b/ws;
endfunction

function r = rec_mean1(x, ws)   # completely filter based, slow, cause lots of 
b(i)=0
  a = [1.,-1.];                 # median filter by recursion IIR
  b = zeros(ws+1,1);            # solution 1 by a. naessl
  b(1) = 1./ws;
  b(ws+1) = -1./ws;
  r = filter(b,a,x);
endfunction

# my own favorite:

function r = rec_mean2(x, ws)    # equivalent, but much faster than rec_mean1, 
not dependend on filter-length
  l = length(x);                 # solution 2 by a. naessl
  x0 = (1./ws).*(postpad(x,l+ws).-prepad(x,l+ws));
  r = filter([1.],[1.,-1.],x0);
endfunction

function r = rec_mean3(x, ws)    # trial: try to get
  l = length(x);                 # initial state /comparable to francescos 
solution
  x1 = x(ws+1:l);                # where does the 1 sample time shift come from 
?
  x2 = x(1:l-ws);                # ... not really slow
  x3 = (1./ws).*(x1.-x2);        # solution 3 by a. naessl
  si = mean(x(1:ws));
  r = filter([1.],[1.,-1.],x3,si);
endfunction


ws = 500                # filter length
a = rand(1,1000000);
f = ones(1,ws)/ws;

tic; m1 = meanfilter(a, ws); toc
tic; m2 = conv(a, f); toc       # ### solution with conv. by Soren Hauberg
tic; m3 = rec_mean1(a, ws); toc
tic; m4 = rec_mean2(a, ws); toc
tic; m5 = rec_mean3(a, ws); toc

# gplot m1', m2', m3', m4', m5', a' w d

# ans = 0.28900    # oct. 2.1.50  winxp, 2nd run
# ans = 4.5480
# ans = 8.9120
# ans = 0.27800
# ans = 0.30800



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------