[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Roundoff error in SMOOTH
- Subject: Roundoff error in SMOOTH
- From: landsman(at)stars.gsfc.nasa.gov (Wayne Landsman)
- Date: 26 May 1999 15:49 EST
- Distribution: world
- News-Software: VAX/VMS VNEWS 1.50AXP
- Newsgroups: comp.lang.idl-pvwave
- Organization: NASA Goddard Space Flight Center -- Greenbelt, Maryland USA
- Xref: news.doit.wisc.edu comp.lang.idl-pvwave:14918
I was recently surprised to discover that applying the SMOOTH function to a
non-negative array could yield an array with negative numbers. I give an
example below. This problem is evidently due to some sort of roundoff error,
since it does not occur when using double precision. But it is not obvious to
me how averaging 9 non-negative numbers (for a 3x3 box smooth) could yield a
negative number, even allowing for roundoff error.
Although not obvious from my simple example, this has nothing to do with edge
effects -- I originally found the problem when 3x3 smoothing a 1024 x 1024
array.
--Wayne Landsman landsman@mpb.gsfc.nasa.gov
d = fltarr(11,11) ;Initialize array to zero
d[ [0,2,0,2],[2,3,4,5] ] = 0.48 ;Set a few values to a positive number
print,where(d lt 0) ;Make sure there are really no negative values
; -1
print,d[1:6,1:6] ;Print a block of numbers in array
; 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
; 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
; 0.00000 0.480000 0.00000 0.00000 0.00000 0.00000
; 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
; 0.480000 0.00000 0.00000 0.00000 0.00000 0.00000
; 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
dd = smooth(d,3) ;Smooth the array
print,dd[1:6,1:6] ;After smooothing there are negative values!
; 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
; 0.0533333 0.0533333 0.0533333 0.00000 0.00000 0.00000
; 0.106667 0.0533333 0.0533333 0.00000 0.00000 0.00000
; 0.160000 0.106667 0.0533333 -6.62274e-09 -6.62274e-09 -6.62274e-09
; 0.106667 0.0533333 0.00000 0.00000 0.00000 0.00000
; 0.0533333 0.0533333 0.00000 0.00000 0.00000 0.00000
dd = smooth(double(d),3) ;Negative values do not occur in double precision
print,float(dd[1:6,1:6])
; 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
; 0.0533333 0.0533333 0.0533333 0.00000 0.00000 0.00000
; 0.106667 0.0533333 0.0533333 0.00000 0.00000 0.00000
; 0.160000 0.106667 0.0533333 0.00000 0.00000 0.00000
; 0.106667 0.0533333 0.00000 0.00000 0.00000 0.00000
; 0.0533333 0.0533333 0.00000 0.00000 0.00000 0.00000