comp.lang.idl-pvwave archive: archive » Re: algorithm question. Can I get rid of the for loop?

Home » Public Forums » archive » Re: algorithm question. Can I get rid of the for loop?

Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend

Re: algorithm question. Can I get rid of the for loop? [message #83749 is a reply to message #83748]

Fri, 22 March 2013 10:01

cgguido
Messages: 195
Registered: August 2005

Senior Member

What about using the fact that Median neglects NaNs?

Stick a NaN in 'Y' for every 'X' you are missing... then median it straight up?

G

On Thursday, March 21, 2013 4:32:15 PM UTC-5, Søren Frimann wrote:
> Hi All.
>
>
>
> I have an implementation of a hampel filter (see e.g. http://exploringdatablog.blogspot.dk/2012/01/moving-window-f ilters-and-pracma.html) in IDL.
>
>
>
> My implementation looks like this:
>
>
>
> ##############################################
>
> FUNCTION hampel, x, y, dx, THRESHOLD=threshold
>
>
>
> Compile_Opt idl2
>
>
>
> IF N_Elements(threshold) EQ 0 THEN $
>
> threshold = 3
>
>
>
> ;initialize arrays
>
> s0 = FltArr(N_Elements(y))
>
> y0 = FltArr(N_Elements(y))
>
> yy = y
>
>
>
> FOR i=0,N_Elements(y)-1 DO BEGIN
>
> index = Where((x GE x[i] - dx) AND (x LE x[i] + dx))
>
> y0[i] = Median(y[index]) ; Median filtering
>
> s0[i] = 1.4826*Median(Abs(y[index] - y0[i])) ;estimating uncertainty
>
> ENDFOR
>
>
>
> ol = Where(Abs(y - y0) GE threshold*s0) ;Index of outliers
>
> yy[ol] = y0[ol]
>
>
>
> result = Create_Struct('y',yy, $
>
> 'sigma',s0)
>
>
>
> RETURN, result
>
>
>
> END
>
> ##############################################
>
>
>
> the filter runs a moving window of width 2*dx measured in the same units as x.
>
> x is generally not uniformly spaced (so there's not a constant number of points inside the window as it moves).
>
> x and y can be quite long vectors so the filter takes a long time to run.
>
> Can anyone see any method for speeding the code up?
>
> Any help would be much appreciated!
>
>
>
> Cheers,
>
> Søren

Report message to a moderator

[Message index]

		Re: algorithm question. Can I get rid of the for loop? By: Jeremy Bailin on Wed, 27 March 2013 11:55
		Re: algorithm question. Can I get rid of the for loop? By: Jeremy Bailin on Wed, 27 March 2013 10:00
		Re: algorithm question. Can I get rid of the for loop? By: Jeremy Bailin on Wed, 27 March 2013 09:53
		Re: algorithm question. Can I get rid of the for loop? By: Søren Frimann on Wed, 27 March 2013 06:31
		Re: algorithm question. Can I get rid of the for loop? By: bobgstockwell on Fri, 22 March 2013 10:29
		Re: algorithm question. Can I get rid of the for loop? By: cgguido on Fri, 22 March 2013 10:01
		Re: algorithm question. Can I get rid of the for loop? By: Heinz Stege on Fri, 22 March 2013 06:25
		Re: algorithm question. Can I get rid of the for loop? By: Jeremy Bailin on Thu, 21 March 2013 20:39

Previous Topic:	Rotating Tickmark labels
Next Topic:	Re: Inscribed circle

-=] Back to Top [=-

[ Syndicate this forum (XML) ] [

] [

]

Current Time: Fri Oct 10 10:36:12 PDT 2025

Total time taken to generate the page: 0.72405 seconds