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 4. 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.

 

#$=head1 NAME
#$
#$pefest  - find prediction error filter, avoid bursty noise
#$
#$=head1 SYNOPSIS
#$
#$C<call pefest1(dd,aa,niter)> 
#$
#$=head1 INPUT PARAMETERS
#$
#$=over 4
#$
#$=item dd  -  C<real(:)>   
#$
#$      Input data
#$
#$=item niter - integer
#$
#$      Number of itterations
#$
#$=back
#$
#$=head1 OUTPUT PARAMETERS
#$
#$=over 4
#$
#$=item  aa - type(filter)
#$
#$       output filter
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Find prediction-error filter (helix magic) avoid bursts
#$
#$=head1 SEE ALSO
#$
#$L<npef>,L<hconest>,L<solver>,L<cgstep>,L<pef>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$

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) {
    integer, intent( in)         :: niter
    real, dimension( :), pointer :: yy
    type( filter)                :: aa
    real, dimension( size( yy))  :: rr
    real                         :: rbar
    integer                      :: stat
    call helicon_init( aa)                          # starting guess
    stat = helicon_lop( .false., .false., yy, rr)
    rbar = quantile( size( yy)/3, abs( rr))         # rbar=(r safe below rbar)
    where( aa%mis)  yy = 0.
    call find_mask(( yy /= 0. .and. abs( rr) < 5 * rbar), aa)
    call find_pef( yy, aa, 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 helicon
  use quantile_mod
contains
  subroutine fixbad1 (niter, aa, yy) {
    integer,       intent (in)	    :: niter
    type( filter), intent (in)      :: aa
    real,    dimension (:)	    :: yy
    real,    dimension (size (yy))  :: rr
    logical, dimension (size (yy))  :: known
    integer                         :: stat
    call helicon_init( aa)
    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, known, .true.)
 }
}