;+
;NAME:
;     WIENER
;PURPOSE:
;     Derives an optimal Wiener filter for a data set
;CATEGORY:
;CALLING SEQUENCE:
;     filter = wiener(data)
;INPUTS:
;     data = array (1D or 2D)  to be filtered
;OPTIONAL INPUT PARAMETERS:
;KEYWORD PARAMETERS
;     /quiet = work quietly
;OUTPUTS:
;     filter = the filter
;COMMON BLOCKS:
;SIDE EFFECTS:
;EXAMPLE:
;     IDL> filter = wiener(image)
;     IDL> f = fft(image,-1)
;     IDL> new_image = fft(f*filter,+1)
;RESTRICTIONS:
;     The size of the array must be even.  If the array is 2D, the size of 
;     both dimentions must be even.
;
;     For 2D arrays, the x and y directions are treated as the same.
;PROCEDURE:
;MODIFICATION HISTORY:
;     12-May-1994 T. Metcalf
;-

function wiener,a,fft=afft,quiet=quiet

   sa = size(a) & nsa = n_elements(sa)

   ;pa = abs(fft(a-mean(a),-1))
   if n_elements(afft) NE sa(nsa-1) then pa = abs(fft(a,-1)) else pa=abs(afft)

   if sa(0) EQ 2 then begin
      if (sa(1) MOD 2L) NE 0L OR (sa(2) MOD 2L) NE 0L then $
         message,'array must have even dimensions'
      nx = sa(1)/2L
      ny = sa(2)/2L
      f = (dist(sa(1),sa(2)))(0:nx-1L,0:ny-1L)
      pa = pa(0:nx-1L,0:ny-1L)
   endif else begin 
      if (sa(1) MOD 2L) NE 0 then message,'array must have even dimension'
      nx = sa(1)/2L
      f = [lindgen(nx-1L)]
      pa = pa(0:nx-1L)
   endelse

   n = n_elements(f)
   n2 = n/2L
   sf = sort(f)

   fs = f(sf)
   power = 1.0
   pa = pa(sf)
   ls = lindgen(n)
   unscramble = sort(ls(sf))

   if NOT keyword_set(quiet) then plot,fs,pa,/yst,/ylog

   ncoeff = polyfitw(fs(n2:*),pa(n2:*),sqrt(1.0+findgen(n2)),1)
   noise  = ncoeff(0) + fs*ncoeff(1)
   ;noise = noise > (max(noise)/2.0)
   scoeff = polyfitw(fs,alog((pa-noise)>(min(pa)/1000.)), $
                     sqrt(1.0/(findgen(n)+1.0)),3)
   signal = (exp(scoeff(0)+fs*(scoeff(1)+fs*(scoeff(2)+fs*scoeff(3))))) > 0.0

   if NOT keyword_set(quiet) then begin
      oplot,fs,noise,linestyle=1,color=!d.table_size/2
      oplot,fs,signal,linestyle=2,color=!d.table_size/2
      oplot,fs,signal+noise,linestyle=3,color=!d.table_size/2
   endif

   w = signal^2/(signal^2+noise^2)
   w = w(unscramble)
   if sa(0) EQ 2 then begin
      w = reform(w,nx,ny)
      w2 = make_array(size=sa,/float)
      w2(0:nx-1L,0:ny-1L) = w
      w2(nx:*,0:ny-1L) = rotate(w,5)
      w2(0:nx-1L,ny:*) = rotate(w,7)
      w2(nx:*,ny:*) = rotate(w,6)
      w = temporary(w2)
   endif else begin
      w = [w,rotate(w,2)]
   endelse

   return,w

end