FUNCTION MEDSMOOTH,ARRAY,WINDOW ;+ ; NAME: ; MEDSMOOTH ; ; PURPOSE: ; Median smoothing of a vector, including points near its ends. ; ; CALLING SEQUENCE: ; SMOOTHED = MEDSMOOTH( VECTOR, WINDOW_WIDTH ) ; ; INPUTS: ; VECTOR = The (1-d numeric) vector to be smoothed ; WINDOW = Odd integer giving the full width of the window over which ; the median is determined for each point. (If WINDOW is ; specified as an even number, then the effect is the same as ; using WINDOW+1) ; ; OUTPUT: ; Function returns the smoothed vector ; ; PROCEDURE: ; Each point is replaced by the median of the nearest WINDOW of points. ; The width of the window shrinks towards the ends of the vector, so that ; only the first and last points are not filtered. These points are ; replaced by forecasting from smoothed interior points. ; ; EXAMPLE: ; Create a vector with isolated high points near its ends ; IDL> a = randomn(seed,40) & a[1] = 10 & a[38] = 10 ; Now do median smoothing with a 7 point window ; IDL> b = medsmooth(a,7) ; Note that, unlike MEDIAN(), that MEDSMOOTH will remove the isolated ; high points near the ends. ; REVISION HISTORY: ; Written, H. Freudenreich, STX, 12/89 ; H.Freudenreich, 8/90: took care of end-points by shrinking window. ; Speed up using vector median when possible W. Landsman February 2002 ;- LEND = N_ELEMENTS(ARRAY)-1 IF (LEND+1) LT WINDOW THEN BEGIN message,/CON, \$ 'ERROR - Size of smoothing window must be smaller than array size' RETURN,ARRAY ENDIF OFFSET = FIX(WINDOW/2) smoothed = median(array, window ) ; Fix the ends: NUMLOOP = (WINDOW-1)/2 - 1 IF NUMLOOP GT 0 THEN BEGIN FOR J=1,NUMLOOP DO BEGIN LEN = 2*J+1 SMOOTHED[J] = MEDIAN(ARRAY[0:LEN-1]) SMOOTHED[LEND-J] = MEDIAN(ARRAY[LEND-LEN+1:LEND]) ENDFOR ENDIF ; Now replace the very last and first points: Y0 = 3.*ARRAY[0]-2.*ARRAY[1] ; Predicted value of point -1 SMOOTHED[0] = MEDIAN([Y0,ARRAY[0],ARRAY[1]]) Y0 = 3.*ARRAY[LEND]-2.*ARRAY[LEND-1] ; Predicted value of point LEND+1 SMOOTHED[LEND] = MEDIAN([Y0,ARRAY[LEND],ARRAY[LEND-1]]) RETURN,SMOOTHED END