! !======================================================================= ! MODULE SMOOTHING ! ! This module contains smoothing routines for curves USE ACCURACY_REAL ! CONTAINS ! !======================================================================= ! SUBROUTINE SMOOFT(Y,N,PTS) ! ! This subroutine smoothes an array Y of length N, with a window ! whose full width is of order PTS neighboring points ! ! Based on J-P Moreau's implementation of a routine from ! "Numerical Recipes" by W.H. Press, B. P. Flannery, ! S.A. Teukolsky and W.T. Vetterling, Cambridge ! University Press, 1986 ! ! Note: the grid is assumed to be equally spaced ! ! ! Input parameters: ! ! * Y : function to be smoothed ! * N : size of array Y ! * PTS : number of neighbouring points taken into account ! ! ! Output parameters: ! ! * Y : smoothed function ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ZERO,ONE,FOURTH ! IMPLICIT NONE ! INTEGER, PARAMETER :: NMAX = 2048 ! double of maximum sample size ! INTEGER, INTENT(IN) :: N,PTS ! INTEGER :: LOGF INTEGER :: M,NMIN INTEGER :: J,MO2,K ! REAL (WP), INTENT(INOUT) :: Y(NMAX) ! REAL (WP) :: CONST,Y1,YN,RN1 REAL (WP) :: FAC ! REAL (WP) :: FLOAT,MAX ! LOGF = 6 ! log file unit ! M = 2 ! NMIN = N + 2 * PTS ! ! 1 IF(M < NMIN) THEN ! M = 2 * M ! GO TO 1 ! END IF ! ! IF(M > NMAX) THEN ! WRITE(LOGF,10) ! STOP ! END IF ! ! CONST = FLOAT((PTS / M)**2) ! Y1 = Y(1) ! YN = Y(N) ! RN1 = ONE / FLOAT(N-1) ! ! DO J = 1, N ! \ Y(J) = Y(J) - RN1 * ( Y1 * FLOAT(N-J) + YN* FLOAT(J-1) ) ! > remove linear trend END DO ! / ! IF(N+1 <= M) THEN ! DO J = N+1, M ! Y(J) = ZERO ! END DO ! END IF ! ! MO2 = M / 2 ! ! CALL REALFT(Y,MO2,1) ! Fourier transform ! Y(1) = Y(1) / FLOAT(MO2) ! FAC = ONE ! ! DO J = 1, MO2-1 ! K = 2 * J + 1 ! IF(FAC /= ZERO) THEN ! FAC = MAX(ZERO,(ONE - CONST * FLOAT(J)**2) / MO2) ! Y(K ) = FAC * Y(K) ! Y(K+1) = FAC * Y(K+1) ! ELSE ! Y(K) = ZERO ! Y(K+1) = ZERO ! END IF ! END DO ! ! FAC = MAX(ZERO,(ONE - FOURTH * FLOAT(PTS**2)) / MO2) ! last point Y(2) = FAC * Y(2) ! ! CALL REALFT(Y,MO2,-1) ! inverse Fourier transform ! DO J = 1, N ! restore linear trend Y(J) = RN1 * (Y1 * FLOAT(N-J) + YN * FLOAT(J-1)) + Y(J) ! END DO ! ! ! Formats ! 10 FORMAT(5X,'<<<<< SAMPLE TOO LARGE: INCREASE NMAX >>>>>',/, & 5X,'<<<<< IN SUBROUTINE SMOOFT (UTILITIES) >>>>>',//) ! CONTAINS ! !----------------------------------------------------------------------- ! SUBROUTINE FOUR1(DATA,NN,ISIGN) ! ! This subroutine replaces DATA(1:2*NN) by: ! ! * its discrete Fourier transform if ISIGN is 1 ! * NN times its inverse discrete Fourier transform if ISIGN is -1 ! ! DATA is either: ! ! * a complex array of length NN ! * a real array of length 2*NN ! ! Warning: NN MUST be an integer power of 2 (this is not checked for!) ! ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ZERO,ONE,TWO,HALF USE PI_ETC, ONLY : PI ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: NN,ISIGN ! INTEGER :: N,J,I,M INTEGER :: MMAX,ISTEP ! REAL (WP), INTENT(INOUT) :: DATA(2*NN) ! REAL (WP) :: TEMPR,TEMPI REAL (WP) :: THETA,WI,WPI,WPR,WR,WTEMP !for the trigonometric recurrences ! REAL (WP) :: FLOAT,SIN ! N = 2 * NN ! J = 1 ! ! DO I = 1, N, 2 ! this is the bit-reversal section of the routine IF(J > I)THEN ! TEMPR = DATA(J) ! exchange the two complex numbers TEMPI = DATA(J+1) ! DATA(J) = DATA(I) ! DATA(J+1) = DATA(I+1) ! DATA(I) = TEMPR ! DATA(I+1) = TEMPI ! END IF ! M = NN ! 1 IF((M > 2) .AND. (J > M)) THEN ! J = J - M ! M = M / 2 ! GO TO 1 ! END IF ! J = J + M ! END DO ! ! mmax = 2 ! here begins the Danielson-Lanczos ! ! section of the routine 2 IF (N > MMAX) THEN ! outer loop executed log2 NN times. ! ! ISTEP = 2 * MMAX ! THETA = TWO * PI / FLOAT(ISIGN * MMAX) ! initialize for the trigonometric recurrence. WPR = - TWO * SIN(HALF * THETA)**2 ! WPI = SIN(THETA) ! WR = ONE ! WI = ZERO ! ! ! DO M =1, MMAX, 2 ! here are the two nested inner loops. DO I=M,N,ISTEP ! J = I + MMAX ! this is the Danielson-Lanczos formula: TEMPR = WR * DATA(J) - WI * DATA(J+1) ! TEMPI = WR * DATA(J+1) + WI * DATA(J) ! DATA(J) = DATA(I) - TEMPR ! DATA(J+1) = DATA(I+1) - TEMPI ! DATA(I) = DATA(I) + TEMPR ! DATA(I+1) = DATA(I+1) + TEMPI ! END DO ! WTEMP = WR ! trigonometric recurrence. WR = WR * WPR - WI * WPI + WR ! WI = WI * WPR + WTEMP * WPI + WI ! END DO ! ! ! MMAX = ISTEP ! GO TO 2 ! not yet done. ! END IF ! all done. ! END SUBROUTINE FOUR1 ! !----------------------------------------------------------------------- ! SUBROUTINE REALFT(DATA,N,ISIGN) ! ! This subroutine calculates the Fourier transform of a set of ! N real-valued data points.Replaces this data (which is stored in ! array DATA(1:N)) by the positive frequency half of its complex Fourier ! transform. The real-valued first and last components of the ! complex transform are returned as elements DATA(1) and DATA(2), ! respectively. N must be a power of 2. This routine also calculates ! the inverse transform of a complex data array if it is the transform of real ! data. (Result in this case must be multiplied by 2/N.) ! ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ONE,TWO,HALF USE PI_ETC, ONLY : PI ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: N,ISIGN INTEGER :: I,I1,I2,I3,I4,N2P3 ! REAL (WP), INTENT(INOUT) :: DATA(N) ! REAL (WP) :: C1,C2 REAL (WP) :: H1I,H1R,H2I,H2R REAL (WP) :: WIS,WRS REAL (WP) :: THETA,WI,WR ! for the trigonometric REAL (WP) :: WPI,WPR,WTEMP ! recurrences ! THETA = PI / FLOAT(N / 2) ! initialize the recurrence. C1 = HALF ! ! IF(ISIGN == 1) THEN ! C2 = - HALF ! CALL FOUR1(DATA,N/2,+1) ! the forward transform is here. ELSE ! C2 = HALF ! otherwise set up for an inverse transform. THETA = - THETA ! END IF ! ! WPR = - TWO * SIN(HALF * THETA)**2 ! WPI = SIN(THETA) ! WR = ONE + WPR ! WI = WPI ! N2P3 = N + 3 ! ! DO I = 2, N / 4 ! case I=1 done separately below. I1 = 2 * I - 1 ! I2 = I1 + 1 ! I3 = N2P3 - I2 ! I4 = I3 + 1 ! WRS = WR ! WIS = WI ! H1R = C1 * (DATA(I1) + DATA(I3)) ! \ H1I = C1 * (DATA(I2) - DATA(I4)) ! > the two separate transforms H2R = - C2 * (DATA(I2) + DATA(I4)) ! > are separated out of DATA H2I = C2 * (DATA(I1) - DATA(I3)) ! / ! DATA(I1) = H1R + WRS * H2R - WIS * H2I ! \ DATA(I2) = H1I + WRS * H2I + WIS * H2R ! > here they are recombined DATA(I3) = H1R - WRS * H2R + WIS * H2I ! > to form the true transform DATA(I4) = - H1I + WRS * H2I + WIS * H2R ! / of the original real data ! WTEMP = WR ! the recurrence. WR = WR * WPR - WI * WPI + WR ! WI = WI * WPR + WTEMP * WPI + WI ! END DO ! ! IF(ISIGN == 1) THEN ! H1R = DATA(1) ! DATA(1) = H1R + DATA(2) ! squeeze the first and last data together DATA(2) = H1R - DATA(2) ! to get them all within the original array ! ELSE ! H1R = DATA(1) ! DATA(1) = C1 * (H1R + DATA(2)) ! DATA(2) = C1 * (H1R - DATA(2)) ! CALL FOUR1(DATA,N/2,-1) ! This is the inverse transform END IF ! for the case ISIGN = -1 ! END SUBROUTINE REALFT ! END SUBROUTINE SMOOFT ! !======================================================================= ! SUBROUTINE TSAVGOL(Y,N) ! ! This subroutine smoothes an array Y of length N, using ! a Savitzky-Golay filter ! ! Based on J-P Moreau's implementation of a routine from ! "Numerical Recipes" by W.H. Press, B. P. Flannery, ! S.A. Teukolsky and W.T. Vetterling, Cambridge ! University Press, 1986 ! ! Note: the grid is assumed to be equally spaced ! ! ! Input parameters: ! ! * Y : function to be smoothed ! * N : size of array Y ! ! ! Output parameters: ! ! * Y : smoothed function ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ZERO ! IMPLICIT NONE ! INTEGER, PARAMETER :: NMAX = 2048 ! double of maximum sample size INTEGER, PARAMETER :: NP = 1000 ! ! INTEGER, INTENT(IN) :: N ! INTEGER :: NL,NR,M INTEGER :: INDEX(NP) INTEGER :: I,J ! REAL (WP), INTENT(INOUT) :: Y(NMAX) ! REAL (WP) :: YSAVE(NMAX) REAL (WP) :: C(NP) ! YSAVE = Y ! save unsmoothed signal ! NL = 5 ! \ NR = 5 ! > see SAVGOL M = 4 ! / ! INDEX(1) = 0 ! seek shift index for given case NL, NR, M ! ! (see SAVGOL) J = 3 ! DO I = 2, NL+1 ! INDEX(I) = I - J ! J = J + 2 ! END DO ! ! ! (see SAVGOL) J = 2 ! DO I = NL+2, NL+NR+1 ! INDEX(I) = I - J ! J = J + 2 ! END DO ! ! ! Calculate Savitzky-Golay filter coefficients ! CALL SAVGOL(C,NL+NR+1,NL,NR,0,M) ! ! ! Apply filter to input data ! DO I = 1, N - NR ! Y(I) = ZERO ! DO J = 1, NL + NR + 1 ! IF(I + INDEX(J) > 0) THEN ! skip left points that do not exist Y(I) = Y(I) + C(J) * YSAVE(I+INDEX(J)) ! END IF ! END DO ! END DO ! ! CONTAINS ! !----------------------------------------------------------------------- ! SUBROUTINE SAVGOL(C,NP,NL,NR,LD,M) ! ! This subroutine returns in C(1:NP), in wrap-around order ! (see reference) consistent with the argument RESPNS ! in routine CONVLV, a set of Savitzky-Golay filter coefficients. ! ! ! Based on J-P Moreau's implementation of a routine from ! "Numerical Recipes" by W.H. Press, B. P. Flannery, ! S.A. Teukolsky and W.T. Vetterling, Cambridge ! University Press, 1986 ! ! ! Input parameters: ! ! * NL : number of data points to the left of each point > to include in the filter ! * NR : number of data points to the right of each point > ! the total number of data points used NL + NR + 1. ! * LD : order of the derivative desired ! LD = 0 for smoothed function ! LD = 1 for smoothed first derivative of the function ! * M : order of the smoothing polynomial ! equal to the highest conserved moment ! usual values are M = 2 or M = 4 ! Lower values for M will produce smoother results ! but may introduce bias ! Higher values for M will reduce the filter bias ! but may "over fit" the data and give a noisier result ! ! ! Output parameters: ! ! * C : Savitzky-Golay filter coefficients ! * NP : size of the C array ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ZERO,ONE ! IMPLICIT NONE ! INTEGER, PARAMETER :: MMAX = 6 ! INTEGER, INTENT(IN) :: NL,NR,LD,M INTEGER, INTENT(IN) :: NP ! INTEGER :: LOGF INTEGER :: D,ICODE,IMJ,IPJ INTEGER :: J,K,KK,MM INTEGER :: INDX(MMAX+1) ! REAL (WP), INTENT(OUT) :: C(NP) ! REAL (WP) :: FAC,SUM REAL (WP) :: A(MMAX+1,MMAX+1),B(MMAX+1) ! REAL (WP) :: FLOAT,MIN ! LOGF = 6 ! log file unit ! ! Testing the arguments ! IF(NP < NL+NR+1 .OR. NL < 0 .OR. NR < 0 .OR. LD > M .OR. & M > MMAX .OR. NL+NR < M) THEN ! WRITE(LOGF,10) STOP END IF ! ! ! Set up the normal equations of the desired least squares fit ! DO IPJ = 0, 2 * M ! SUM = ZERO ! IF(IPJ == 0) SUM = ONE ! DO K = 1, NR ! SUM = SUM + FLOAT(K)**IPJ ! END DO ! DO K = 1, NL ! SUM = SUM + FLOAT(-K)**IPJ ! END DO ! MM = MIN(IPJ,2*M-IPJ) ! DO IMJ = -MM, MM, 2 ! A(1+(IPJ+IMJ)/2,1+(IPJ-IMJ)/2) = SUM ! END DO ! END DO ! ! ! Solve them: LU decomposition. ! CALL LUDCMP(A,M+1,MMAX+1,INDX,D,ICODE) ! ! DO J = 1, M+1 ! B(J) = ZERO ! END DO ! ! ! Right-hand side vector is unit vector, ! depending on which derivative we want ! B(LD+1) = ONE ! ! ! Backsubstitute, giving one row of the inverse matrix ! CALL LUBKSB(A,M+1,MMAX+1,INDX,B) ! ! ! Zero the output array (it may be bigger ! than the number of coefficients) ! DO KK = 1 ,NP ! C(KK) = ZERO ! END DO ! ! ! Each Savitzky-Golay coefficient is the dot product ! of powers of an integer with the inverse matrix row ! DO K = -NL, NR ! SUM = B(1) ! FAC = ONE ! DO MM = 1, M ! FAC = FAC * K ! SUM = SUM + B(MM+1) * FAC ! END DO ! KK = MOD(NP-K,NP) + 1 ! Store in wrap-around order C(KK) = SUM ! END DO ! ! ! Formats: ! 10 FORMAT(5X,'<<<<< BAD ARGS IN SAVGOL >>>>>',//) ! ! END SUBROUTINE SAVGOL ! !----------------------------------------------------------------------- ! SUBROUTINE LUDCMP(A,N,NP,INDX,D,CODE) ! ! Given an N x N matrix A, this routine replaces it by the LU ! decomposition of a rowwise permutation of itself. ! ! ! Input parameters: ! ! * A : input matrix ! * N : dimensioning of matrix A ! * NP : physical dimension of matrix A ! ! ! Output parameters: ! ! * INDX : output vector recording the row permutation ! effected by the partial pivoting ! * D : = 1 number of row interchanges even ! = -1 number of row interchanges odd ! * CODE : return code (1 matrix is singular) ! ! This routine is used in combination with LUBKSB ! to solve linear equations or to invert a matrix. ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ZERO,ONE ! IMPLICIT NONE ! INTEGER, PARAMETER :: NMAX = 100 ! INTEGER, INTENT(IN) :: N,NP INTEGER, INTENT(OUT) :: INDX(N),D,CODE ! INTEGER :: AMAX,DUM,SUM INTEGER :: VV(NMAX) INTEGER :: I,J,K INTEGER :: IMAX ! REAL (WP), INTENT(OUT) :: A(NP,NP) ! REAL (WP) :: ABS ! REAL (WP), PARAMETER :: TINY = 1.0E-12_WP ! D = 1 ! CODE = 0 ! DO I = 1, N ! AMAX = ZERO ! DO J = 1, N ! IF(ABS(A(I,J)) > AMAX) AMAX = ABS(A(I,J)) ! END DO ! IF(AMAX < TINY) THEN ! CODE = 1 ! RETURN ! END IF ! VV(I) = ONE / AMAX ! END DO ! ! DO J = 1, N ! DO I = 1, J - 1 ! SUM = A(I,J) ! DO K = 1, I-1 ! SUM = SUM - A(I,K) * A(K,J) ! END DO ! A(I,J) = SUM ! END DO ! AMAX = ZERO ! DO I = J, N ! SUM = A(I,J) ! DO K = 1,J - 1 ! SUM = SUM - A(I,K) * A(K,J) ! END DO ! A(I,J) = SUM ! DUM = VV(I) * ABS(SUM) ! IF(DUM >= AMAX) THEN ! IMAX = I ! AMAX = DUM ! END IF ! END DO ! ! IF(J /= IMAX) THEN ! DO K = 1, N ! DUM = A(IMAX,K) ! A(IMAX,K) = A(J,K) ! A(J,K) = DUM ! END DO ! D = - D ! VV(IMAX) = VV(J) ! END IF ! ! INDX(J) = IMAX ! IF(ABS(A(J,J)) < TINY) A(J,J) = TINY ! ! IF(J /= N) THEN ! DUM = ONE / A(J,J) ! DO I = J + 1, N ! A(I,J) = A(I,J) * DUM ! END DO ! END IF ! END DO ! ! END SUBROUTINE LUDCMP ! !----------------------------------------------------------------------- ! SUBROUTINE LUBKSB(A,N,NP,INDX,B) ! ! This subroutine solves the set of N linear equations A * X = B ! ! ! Input parameters: ! ! * A : input LU decomposition from routine LUDCMP ! * N : same as in routine LUDCMP ! * NP : ! * INDX : same as in routine LUDCMP ! * B : right-handside vector ! ! ! Output parameters: ! ! * B : solution vector X ! ! ! This routine is also efficient for plain matrix inversion. ! ! ! Last modified (DS) : 16 Dec 2020 ! ! USE REAL_NUMBERS, ONLY : ZERO ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: N,NP INTEGER, INTENT(INOUT) :: INDX(N) ! INTEGER :: I,II,J,LL ! REAL (WP), INTENT(INOUT) :: A(NP,NP) REAL (WP), INTENT(INOUT) :: B(N) ! REAL (WP) :: SUM ! II = 0 ! ! DO I = 1, N ! LL = INDX(I) ! SUM = B(LL) ! B(LL) = B(I) ! IF(II /= 0) THEN ! DO J = II, I-1 ! SUM = SUM - A(I,J) * B(J) ! END DO ! ELSE IF(SUM /= ZERO) THEN ! II = I ! END IF ! B(I) = SUM ! END DO ! ! DO I = N, 1, -1 ! SUM = B(I) ! IF(I < N) THEN ! DO J = I+1, N ! SUM = SUM - A(I,J) * B(J) ! END DO ! END IF ! B(I) = SUM / A(I,I) ! END DO ! ! END SUBROUTINE LUBKSB ! END SUBROUTINE TSAVGOL ! END MODULE SMOOTHING