Program c09ezfe
! C09EZF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
! .. Use Statements ..
Use nag_library, Only: c09abf, c09ecf, c09edf, c09eyf, c09ezf, 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, i, ifail, ilev, j, &
lda, ldb, ldd, lenc, m, n, nf, nwcn, &
nwct, nwl
Character (10) :: mode, wavnam, wtrans
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:,:), an(:,:), b(:,:), c(:), &
d(:,:), e(:,:)
Integer, Allocatable :: dwtlvm(:), dwtlvn(:)
Integer :: icomm(180)
! .. Intrinsic Procedures ..
Intrinsic :: abs, log, real, sqrt
! .. Executable Statements ..
Write (nout,*) 'C09EZF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
! Read problem parameters
Read (nin,*) m, n
Read (nin,*) wavnam, mode
Write (nout,99999) wavnam, mode, m, n
! 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
Allocate (a(lda,n),an(lda,n),b(ldb,n),e(m,n))
! Read in the original data
Do i = 1, m
Read (nin,*) a(i,1:n)
End Do
! Output the original data
Write (nout,99997)
Do i = 1, m
Write (nout,99998) a(i,1:n)
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,m,n,thresh)
! Output the noisy data
Write (nout,99996)
Do i = 1, m
Write (nout,99998) an(i,1:n)
End Do
! Query wavelet filter dimensions
! For Multi-Resolution Analysis, decomposition, wtrans = 'M'
wtrans = 'Multilevel'
ifail = 0
Call c09abf(wavnam,wtrans,mode,m,n,nwl,nf,nwct,nwcn,icomm,ifail)
! Allocate arrays to hold the coefficients, C, and the dimensions
! of the coefficients at each level, DWTLVM, DWTLVN
lenc = nwct
Allocate (c(lenc),dwtlvm(nwl),dwtlvn(nwl))
! Perform a forwards multi-level transform on the noisy data
ifail = 0
Call c09ecf(m,n,an,lda,lenc,c,nwl,dwtlvm,dwtlvn,icomm,ifail)
! Reconstruct without thresholding of detail coefficients
ifail = 0
Call c09edf(nwl,lenc,c,m,n,b,ldb,icomm,ifail)
! Calculate the Mean Square Error of the noisy reconstruction
e(:,:) = a(:,:) - b(:,:)
mse = dnrm2(m*n,e,1)
mse = mse**2
mse = mse/real(m*n,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 2D array to hold the detail coefficients
ldd = dwtlvm(nwl)
Allocate (d(ldd,dwtlvn(nwl)))
denoised = 0
! For each level
Do ilev = nwl, 1, -1
! Select detail coefficients
Do cindex = 1, 3
! Extract coefficients into the 2D array D
ifail = 0
Call c09eyf(ilev,cindex,lenc,c,d,ldd,icomm,ifail)
! Perform the hard thresholding operation
Do j = 1, dwtlvn(nwl-ilev+1)
Do i = 1, dwtlvm(nwl-ilev+1)
If (abs(d(i,j))<thresh) Then
d(i,j) = 0.0_nag_wp
denoised = denoised + 1
End If
End Do
End Do
! Insert the denoised coefficients back into C
ifail = 0
Call c09ezf(ilev,cindex,lenc,c,d,ldd,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)
! Reconstruct original data following thresholding of detail coefficients
ifail = 0
Call c09edf(nwl,lenc,c,m,n,b,ldb,icomm,ifail)
! Calculate the Mean Square Error of the denoised reconstruction
e(:,:) = a(:,:) - b(:,:)
mse = dnrm2(m*n,e,1)
mse = mse**2
mse = mse/real(m*n,kind=nag_wp)
Write (nout,99993) mse
! Output the denoised reconstruction
Write (nout,99992)
Do i = 1, m
Write (nout,99998) b(i,1:n)
End Do
99999 Format (1X,' MLDWT :: Wavelet : ',A,/,1X,' End mode : ',A,/, &
1X,' M : ',I4,/,1X,' N : ',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 : ')
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,m,n,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) :: lda, m, n
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: a(lda,n)
Real (Kind=nag_wp), Intent (Out) :: an(lda,n)
! .. Local Scalars ..
Real (Kind=nag_wp) :: var, xmu
Integer :: genid, i, ifail, 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))
! Generate the noise variates
ifail = 0
Do i = 1, n
Call g05skf(m,xmu,var,state,x(1,i),ifail)
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,kind=nag_wp)))
End Subroutine create_noise
End Program c09ezfe