CSSED
Smooths one-dimensional data by error detection.
Required Arguments
XDATA — Array of length NDATA containing the abscissas of the data points. (Input)
FDATA — Array of length NDATA containing the ordinates (function values) of the data points. (Input)
DIS — Proportion of the distance the ordinate in error is moved to its interpolating curve. (Input)
It must be in the range 0.0 to 1.0. A suggested value for DIS is one.
SC — Stopping criterion. (Input)
SC should be greater than or equal to zero. A suggested value for SC is zero.
MAXIT — Maximum number of iterations allowed. (Input)
SDATA — Array of length NDATA containing the smoothed data. (Output)
Optional Arguments
NDATA — Number of data points. (Input)
Default: NDATA = size (XDATA,1).
FORTRAN 90 Interface
Generic: CALL CSSED (XDATA, FDATA, DIS, SC, MAXIT, SDATA [, …])
Specific: The specific interface names are S_CSSED and D_CSSED.
FORTRAN 77 Interface
Single: CALL CSSED (NDATA, XDATA, FDATA, DIS, SC, MAXIT, SDATA)
Double: The double precision name is DCSSED.
Description
The routine CSSED is designed to smooth a data set that is mildly contaminated with isolated errors. In general, the routine will not work well if more than 25% of the data points are in error. The routine CSSED is based on an algorithm of Guerra and Tapia (1974).
Setting NDATA = n, FDATA = f, SDATA = s and XDATA = x, the algorithm proceeds as follows. Although the user need not input an ordered XDATA sequence, we will assume that x is increasing for simplicity. The algorithm first sorts the XDATA values into an increasing sequence and then continues. A cubic spline interpolant is computed for each of the 6-point data sets (initially setting s = f)
(xj, sj) j = i − 3, …, i + 3 j ≠ i,
where i = 4, …, n − 3 using CSAKM. For each i the interpolant, which we will call Si, is compared with the current value of si, and a ‘point energy’ is computed as
pei = Si(xi) − si
Setting sc = SC, the algorithm terminates either if MAXIT iterations have taken place or if
∣pei∣ ≤ sc(xi+3 - xi-3)/6 i = 4,…, n - 3
If the above inequality is violated for any i, then we update the i-th element of s by setting si = si + d(pei), where d = DIS. Note that neither the first three nor the last three data points are changed. Thus, if these points are inaccurate, care must be taken to interpret the results.
The choice of the parameters d, sc and MAXIT are crucial to the successful usage of this subroutine. If the user has specific information about the extent of the contamination, then he should choose the parameters as follows: d = 1, sc = 0 and MAXIT to be the number of data points in error. On the other hand, if no such specific information is available, then choose d = .5, MAXIT ≤2n, and
In any case, we would encourage the user to experiment with these values.
Comments
1. Workspace may be explicitly provided, if desired, by use of C2SED/DC2SED. The reference is:
CALL C2SED (NDATA, XDATA, FDATA, DIS, SC, MAXIT, DATA, WK, IWK)
The additional arguments are as follows:
WK — Work array of length 4 * NDATA + 30.
IWK — Work array of length 2 * NDATA.
2. Informational error
|
Type |
Code |
Description |
|
3 |
1 |
The maximum number of iterations allowed has been reached. |
3. The arrays FDATA and SDATA may the same.
Example
We take 91 uniform samples from the function 5 + (5 + t2 sin t)/t on the interval [1, 10]. Then, we contaminate 10 of the samples and try to recover the original function values.
USE CSSED_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER NDATA
PARAMETER (NDATA=91)
!
INTEGER I, MAXIT, NOUT, ISB(10)
REAL DIS, F, FDATA(91), SC, SDATA(91), SIN, X, XDATA(91),&
RNOISE(10)
INTRINSIC SIN
!
DATA ISB/6, 17, 26, 34, 42, 49, 56, 62, 75, 83/
DATA RNOISE/2.5, -3.0, -2.0, 2.5, 3.0, -2.0, -2.5, 2.0, -2.0, 3.0/
!
F(X) = (X*X*SIN(X)+5.0)/X + 5.0
! EX. #1; No specific information
! available
DIS = 0.5
SC = 0.56
MAXIT = 182
! Set values for XDATA and FDATA
XDATA(1) = 1.0
FDATA(1) = F(XDATA(1))
DO 10 I=2, NDATA
XDATA(I) = XDATA(I-1) + .1
FDATA(I) = F(XDATA(I))
10 CONTINUE
! Contaminate the data
DO 20 I=1, 10
FDATA(ISB(I)) = FDATA(ISB(I)) + RNOISE(I)
20 CONTINUE
! Smooth data
CALL CSSED (XDATA, FDATA, DIS, SC, MAXIT, SDATA)
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99997)
! Write data
DO 30 I=1, 10
WRITE (NOUT,99999) F(XDATA(ISB(I))), FDATA(ISB(I)),&
SDATA(ISB(I))
30 CONTINUE
! EX. #2; Specific information
! available
DIS = 1.0
SC = 0.0
MAXIT = 10
! A warning message is produced
! because the maximum number of
! iterations is reached.
!
! Smooth data
CALL CSSED (XDATA, FDATA, DIS, SC, MAXIT, SDATA)
! Write heading
WRITE (NOUT,99998)
! Write data
DO 40 I=1, 10
WRITE (NOUT,99999) F(XDATA(ISB(I))), FDATA(ISB(I)),&
SDATA(ISB(I))
40 CONTINUE
!
99997 FORMAT (' Case A - No specific information available', /,&
' F(X) F(X)+NOISE SDATA(X)', /)
99998 FORMAT (' Case B - Specific information available', /,&
' F(X) F(X)+NOISE SDATA(X)', /)
99999 FORMAT (' ', F7.3, 8X, F7.3, 11X, F7.3)
END
Case A - No specific information available
F(X) F(X)+NOISE SDATA(X)
9.830 12.330 9.870
8.263 5.263 8.215
5.201 3.201 5.168
2.223 4.723 2.264
1.259 4.259 1.308
3.167 1.167 3.138
7.167 4.667 7.131
10.880 12.880 10.909
12.774 10.774 12.708
7.594 10.594 7.639
*** WARNING ERROR 1 from CSSED. Maximum number of iterations limit MAXIT
*** =10 exceeded. The best answer found is returned.
Case B - Specific information available
F(X) F(X)+NOISE SDATA(X)
9.830 12.330 9.831
8.263 5.263 8.262
5.201 3.201 5.199
2.223 4.723 2.225
1.259 4.259 1.261
3.167 1.167 3.170
7.167 4.667 7.170
10.880 12.880 10.878
12.774 10.774 12.770
7.594 10.594 7.592
