pro vdecr, image_2, im, w
;              image_2 - signal debruite
;              im      - image originale 512 x 512
;              w       - transformee en ondelettes
;                        de la coupe 255 en x de l'image im

   print,'estimation de la vitesse de decroissance des coefficients ...'

   n1 =640
   n2 =512
   pwr=512

   im_init   = bytarr(n1,n2)
   im        = bytarr(pwr,pwr)


   pwr1  = pwr-1
   nbpix = LONG(pwr)*pwr

   print,'number of files to process ?'
   read, nwindow

   nb = nwindow - 1


   name1 = 'spi_c242'

   for i = 0, nb do begin

   number = string( format='(2I)', i)
   no     = strcompress( number,/remove_all)
   ss =  no
   if (i lt 10) then ss = '0' + no

   file = name1 + ss + '.bin'

   openr, 1, file
   readu, 1, im_init
   close,1

   print, 'Read unformatted initial image (640,512) from file ', file

   im = im_init(0:511,0:511)

;---------------------------------------------------------------------
  coupe = 255
  ; coupe en x
  image_1 = im( *, coupe )

  coeffs = 4
  ;2 moments nuls

  w     = wtn(image_1,coeffs)

  ; separer le signal du bruit

  aw = abs( w )
  print,'abs(coefficients d ondelettes):',min(aw),max(aw)

  eps = 0.001
; -----------
  ind = where( aw lt eps )
  aw( ind ) = eps

  z = aw(0:511)
 ;-------------

  ;ordre decroissant
  ;=================
  y = z( reverse( sort( z )))

  p = fltarr(500)
  k = -1
  for i = 10, 509 do begin
  k = k+1
  p(k) = ( alog(i)/alog(2) )/abs( alog( y(i) )/alog(2) )
  endfor

  print,'p:', p(0:20)

  set_plot,'x'
  window,0,xsize=pwr,ysize= pwr 
  plot,/ylog,yrange=[0.01,2.], p(0:240), title='Vitesse de decroissance des ABS(Coefficients d ondelettes)'  

  set_plot,'ps'
  device,filename='dessin1.ps'
  device,xsize=15.,ysize=15.
  plot,/ylog,yrange=[0.01,2.], p(0:240), title='Vitesse de decroissance des ABS(Coefficients d ondelettes)'  
  device,/close

  indices = indgen(512)+1
  xx = alog(indices)/alog(2)
  yy = alog(y)/alog(2)

  set_plot,'x'
  window,1,xsize=pwr,ysize= pwr 
  plot, xx, yy, title='Daub4 wavelet basis', xtitle='Amplitude of the sorted coefficients'

  set_plot,'ps'
  device,filename='dessin2.ps'
  device,xsize=15.,ysize=15.
  plot, xx, yy, title='Daub4 wavelet basis', xtitle='Amplitude of the sorted coefficients'
  device,/close

  ; histogramme
  ; ===========
  coeff = 100./512.
  y = coeff*histogram( z, min=1.,max=30. )
  print, 'histogramme:',min(y),max(y)

  set_plot,'x'
  window,2,xsize=pwr,ysize= pwr 
  plot, y, title='Histogramme des coefficients d ondelettes(1.,30.)',ytitle='%'

  set_plot,'ps'
  device,filename='dessin3.ps'
  device,xsize=15.,ysize=15.
  plot, y, title='Histogramme des coefficients d ondelettes(1.,30.)',ytitle='%'
  device,/close

;----------------------------------------------------------------------------
  rep='c'
  print,'pour continuer,tapez un caractere'
  read,rep

  print,'thres=?'
  read,thres

  ind = where( abs(w) lt thres )
  w( ind ) = 0.0

  image_2 = wtn(w, coeffs, /inverse)

  count = pwr - n_elements( ind )

  print,'The signal is reconstructed from:'
  print,100.*count/pwr,'% of original image data.'

  set_plot,'ps'
  device,filename='dessin4.ps'
  device,xsize=15.,ysize=15.
  plot,image_1,title='Original signal'
  device,/close

  set_plot,'x'
  window,0,xsize=pwr,ysize=pwr
  plot,image_1,title='Original signal'

  window,1,xsize=pwr,ysize=pwr
  plot,image_2,title='Denoised signal'

  set_plot,'ps'
  device,filename='dessin5.ps'
  device,xsize=15.,ysize=15.
  plot,image_2,title='Denoised signal'
  device,/close

  set_plot,'x'
  endfor
  end