DEBURST

We can use the same technique to throw out fitting equations from defective data that we use for missing data. Recall the theory and discussion leading up to Figure . There we identified defective data by its lack of continuity. We used the fitting equations $0\approx w_i (y_{i+1} -2y_i + y_{i-1})$where the weights w_i were chosen to be approximately the inverse to the residual (y_i+1 -2y_i + y_i-1) itself.

Here we will first use the second derivative (Laplacian in 1-D) to throw out bad points, while we determine the PEF. Having the PEF, we use it to fill in the missing data.

 

module pefest {      # Estimate a PEF avoiding zeros and bursty noise on input.
  use quantile_mod
  use helicon
  use misinput
  use pef
contains
  subroutine pefest1 (niter, yy, aa, lag, mx) {
    integer, intent (in)             :: niter
    real,    dimension (:), pointer  :: yy, aa
    integer, dimension (:), pointer  :: lag
    logical, dimension (:), pointer  :: mx
    real,    dimension (size (yy))   :: rr, yfre
    real                             :: rbar
    integer                          :: stat
    call helicon_init (aa, lag)                     # starting guess
    stat = helicon_lop (.false., .false., yy, rr)
    rbar = quantile( size( yy)/3, abs( rr))         # rbar=(r safe below rbar)
    where( yy == 0. .or. abs( rr) > 5 * rbar)
       yfre = 0.                                    # risky
    elsewhere
       yfre = 1.                                    # safe enuf
    where (mx) yy = 0.
    call find_mask( yfre, aa, lag, mx)
    call find_pef( yy, aa, lag, mx, niter)  
  }
}

The result of this ``PEF-deburst'' processing is shown in Figure 6.

 

pefdeburst90
Figure 6 Top is synthetic data with noise spikes and bursts. (Some bursts are fifty times larger than shown.) Next is after running medians. Next is Laplacian filter Cauchy deburst processing. Last is PEF-deburst processing.

Given the PEF that comes out of pefest1(), subroutine fixbad1() below convolves it with the data and looks for anomalous large outputs. For each that is found, the input data is declared defective and set to zero. Then subroutine mis1() is invoked to replace the zeroed values by reasonable ones.  

module fixbad { # Given a PEF, find bad data and restore it.
  use mis2
  use quantile_mod
  use helicon
contains
  subroutine fixbad1 (niter, aa, lag, yy) {
    integer, intent (in)	    :: niter
    real,    dimension (:), pointer :: aa
    integer, dimension (:), pointer :: lag
    real,    dimension (:)	    :: yy
    real,    dimension (size (yy))  :: rr
    logical, dimension (size (yy))  :: known
    integer                         :: stat
    call helicon_init (aa, lag)
    stat = helicon_lop (.false., .false., yy, rr); rr = abs (rr)
    known = ( yy > 0.) .and. ( rr < 4. * quantile( size(rr)/2, rr))
    call mis1 (niter, yy, aa, lag, known, .true., .true.)
 }
}