NAG Library Manual, Mark 28.3
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program c09fzfe

!     C09FZF Example Program Text
!     Mark 28.3 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: c09acf, c09fcf, c09fdf, c09fyf, c09fzf, dnrm2,    &
                             nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: mse, thresh
      Integer                          :: cindex, denoised, fr, i, ifail,      &
                                          ilev, j, k, lda, ldb, ldd, lenc, m,  &
                                          n, nf, nwcfr, nwcn, nwct, nwl, sda,  &
                                          sdb, sdd
      Character (10)                   :: mode, wavnam, wtrans
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:,:), an(:,:,:), b(:,:,:), c(:), &
                                          d(:,:,:), e(:,:,:)
      Integer, Allocatable             :: dwtlvfr(:), dwtlvm(:), dwtlvn(:)
      Integer                          :: icomm(260)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, log, real, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'C09FZF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
!     Read problem parameters
      Read (nin,*) m, n, fr
      Read (nin,*) wavnam, mode
      Write (nout,99999) wavnam, mode, m, n, fr

!     Allocate arrays to hold the original data, A, original data plus noise,
!     AN, reconstruction using denoised coefficients, B, and randomly
!     generated noise, X.
      lda = m
      ldb = m
      sda = n
      sdb = n
      Allocate (a(lda,lda,fr),an(lda,lda,fr),b(ldb,ldb,fr),e(m,n,fr))

!     Read in the original data
      Do k = 1, fr
        Do i = 1, m
          Read (nin,*) a(i,1:n,k)
        End Do
        If (k<fr) Then
          Read (nin,*)
        End If
      End Do

!     Output the original data
      Write (nout,99997)
      Do k = 1, fr
        Write (nout,99991) k
        Do i = 1, m
          Write (nout,99998) a(i,1:n,k)
        End Do
      End Do

!     Fill the array AN with the original data in A plus some noise
!     and return a VisuShrink denoising threshold, thresh.
      Call create_noise(a,an,lda,sda,m,n,fr,thresh)

!     Output the noisy data
      Write (nout,99996)
      Do k = 1, fr
        Write (nout,99991) k
        Do i = 1, m
          Write (nout,99998) an(i,1:n,k)
        End Do
      End Do

!     Query wavelet filter dimensions
!     For Multi-Resolution Analysis, decomposition, wtrans = 'M'
      wtrans = 'Multilevel'
      ifail = 0
      Call c09acf(wavnam,wtrans,mode,m,n,fr,nwl,nf,nwct,nwcn,nwcfr,icomm,      &
        ifail)

!     Allocate arrays to hold the coefficients, C, and the dimensions
!     of the coefficients at each level, DWTLVM, DWTLVN, DWTLVFR
      lenc = nwct
      Allocate (c(lenc),dwtlvm(nwl),dwtlvn(nwl),dwtlvfr(nwl))

!     Perform a forwards multi-level transform on the noisy data
      ifail = 0
      Call c09fcf(m,n,fr,an,lda,sda,lenc,c,nwl,dwtlvm,dwtlvn,dwtlvfr,icomm,    &
        ifail)

!     Reconstruct without thresholding of detail coefficients
      ifail = 0
      Call c09fdf(nwl,lenc,c,m,n,fr,b,ldb,sdb,icomm,ifail)

!     Calculate the Mean Square Error of the noisy reconstruction
      e(:,:,:) = a(:,:,:) - b(:,:,:)
      mse = dnrm2(m*n*fr,e,1)
      mse = mse**2
      mse = mse/real(m*n*fr,kind=nag_wp)
      Write (nout,99995) mse

!     Now perform the denoising by extracting each of the detail
!     coefficients at each level and applying hard thresholding

!     Allocate a 3D array to hold the detail coefficients
      ldd = dwtlvm(nwl)
      sdd = dwtlvn(nwl)
      Allocate (d(ldd,sdd,dwtlvn(nwl)))

      denoised = 0
!     For each level
      Do ilev = nwl, 1, -1

!       Select detail coefficients
        Do cindex = 1, 7

!         Extract coefficients into the 3D array D
          ifail = 0
          Call c09fyf(ilev,cindex,lenc,c,d,ldd,sdd,icomm,ifail)

!         Perform the hard thresholding operation
          Do k = 1, dwtlvfr(nwl-ilev+1)
            Do j = 1, dwtlvn(nwl-ilev+1)
              Do i = 1, dwtlvm(nwl-ilev+1)
                If (abs(d(i,j,k))<thresh) Then
                  d(i,j,k) = 0.0_nag_wp
                  denoised = denoised + 1
                End If
              End Do
            End Do
          End Do

!         Insert the denoised coefficients back into C
          ifail = 0
          Call c09fzf(ilev,cindex,lenc,c,d,ldd,sdd,icomm,ifail)

        End Do

      End Do

!     Output the number of coefficients that were set to zero
      Write (nout,99994) denoised, nwct - dwtlvm(1)*dwtlvn(1)*dwtlvfr(1)

!     Reconstruct original data following thresholding of detail coefficients
      ifail = 0
      Call c09fdf(nwl,lenc,c,m,n,fr,b,ldb,sdb,icomm,ifail)

!     Calculate the Mean Square Error of the denoised reconstruction
      e(:,:,:) = a(:,:,:) - b(:,:,:)
      mse = dnrm2(m*n*fr,e,1)
      mse = mse**2
      mse = mse/real(m*n*fr,kind=nag_wp)
      Write (nout,99993) mse

!     Output the denoised reconstruction
      Write (nout,99992)
      Do k = 1, fr
        Write (nout,99991) k
        Do i = 1, m
          Write (nout,99998) b(i,1:n,k)
        End Do
      End Do

99999 Format (1X,' MLDWT :: Wavelet  : ',A,/,1X,'          End mode : ',A,/,   &
        1X,'          M        : ',I4,/,1X,'          N        : ',I4,/,1X,    &
        '          FR       : ',I4)
99998 Format (8(F8.4,1X),:)
99997 Format (/,1X,' Original data            A  : ')
99996 Format (/,1X,' Original data plus noise AN : ')
99995 Format (/,1X,' Without denoising Mean Square Error is ',F9.6)
99994 Format (/,1X,' Number of coefficients denoised is ',I3,' out of ',I3)
99993 Format (/,1X,' With denoising Mean Square Error is ',F9.6)
99992 Format (/,1X,' Reconstruction of denoised input D : ')
99991 Format (1X,' Frame ',I2,' : ')

    Contains

!     Subroutine fills the output array AN with the data in A
!     plus some noise taken from a normal distribution, and
!     returns the VisuShrink denoising threshold, thresh.

      Subroutine create_noise(a,an,lda,sda,m,n,fr,thresh)

!       .. Use Statements ..
        Use nag_library, Only: g05kff, g05skf
!       .. Parameters ..
        Integer, Parameter             :: lseed = 1
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: thresh
        Integer, Intent (In)           :: fr, lda, m, n, sda
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: a(lda,sda,fr)
        Real (Kind=nag_wp), Intent (Out) :: an(lda,sda,fr)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: var, xmu
        Integer                        :: genid, i, ifail, j, lstate, subid
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: x(:,:,:)
        Integer                        :: seed(lseed)
        Integer, Allocatable           :: state(:)
!       .. Executable Statements ..

!       Set up call to g05skf in order to create some random noise from
!       a normal distribution to add to the original data.
!       Initial call to RNG initializer to get size of STATE array
        seed(1) = 642521
        genid = 3
        subid = 0
        lstate = 0
        Allocate (state(lstate))
        ifail = 0
        Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!       Reallocate STATE
        Deallocate (state)
        Allocate (state(lstate))

!       Initialize the generator to a repeatable sequence
        ifail = 0
        Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!       Set the distribution parameters for the random noise.
        xmu = 0.0_nag_wp
        var = 0.1E-3_nag_wp

        Allocate (x(m,n,fr))

!       Generate the noise variates
        ifail = 0
        Do j = 1, fr
          Do i = 1, n
            Call g05skf(m,xmu,var,state,x(1,i,j),ifail)
          End Do
        End Do

!       Add the noise to the original input and save in AN
        an(:,:,:) = a(:,:,:) + x(:,:,:)

!       Calculate the threshold based on VisuShrink denoising
        thresh = sqrt(var)*sqrt(2._nag_wp*log(real(m*n*fr,kind=nag_wp)))

      End Subroutine create_noise

    End Program c09fzfe