Program c09fzfe
! Mark 25 Release. NAG Copyright 2014.
! .. 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) Read (nin,*)
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 initialiser 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