Open in the MATLAB editor: c09fz_example
function c09fz_example fprintf('c09fz example results\n\n'); % 3D Data m = int64(4); n = int64(4); f = int64(4); a = zeros(m,n,f); a(1:2:m,:,1:2:f) = 0.01; a(2:2:m,:,2:2:f) = 0.01; a(2:2:m,:,1:2:f) = 1; a(1:2:m,:,2:2:f) = 1; % Add Noise genid = int64(3); subid = int64(0); seed(1) = int64(642521); [state, ifail] = g05kf(genid, subid, seed); [state, x, ifail] = g05sk(m*n*f, 0, 1.0e-4, state); an = a + reshape(x,[m,n,f]); fprintf('\nInput data a:\n'); disp(a); fprintf('\nNoisy data an:\n'); disp(an); % Wavelet setup wavnam = 'Haar'; mode = 'Period'; wtrans = 'Multilevel'; [nwl, nf, nwct, nwcn, nwcf, icomm, ifail] = ... c09ac(... wavnam, wtrans, mode, m, n, f); % Multi-level wavelet transform on noisy data [c, dwtlvm, dwtlvn, dwtlvfr, icomm, ifail] = ... c09fc(... n, f, an, nwct, nwl, icomm); % Reconstruct without thresholding of detail coefficients [b, ifail] = c09fd(nwl, c, m, n, f, icomm); % Mean square error of noisy reconstruction mse = (norm(reshape(a-b,[m*n*f,1]))^2)/(double(m*n*f)); fprintf('Without denoising Mean Square Error is %9.6f\n',mse); % De-noise by applying hard threshold to detail coefficients thresh = 0.01*sqrt(2*log(double(m*n*f))); nt = 0; nnt = 0; for ilev = 1:nwl level = int64(nwl - ilev + 1); for detail = int64(1:7) % Extract the selected set of coefficients. [d, icomm, ifail] = c09fy(... level, detail, c, icomm); % Threshold d1 = dwtlvm(ilev); d2 = dwtlvn(ilev); d3 = dwtlvfr(ilev); for i = 1:d1 for j = 1:d2 for k = 1:d3 if abs(d(i,j,k))<thresh d(i,j,k) = 0; nt = nt + 1; end nnt = nnt + 1; end end end % Insert de-noised coefficients back into c [c, icomm, ifail] = c09fz(... level, detail, c, d, icomm); end end fprintf('\nNumber of coefficients denoised is %3d out of %3d\n',nt,nnt); % Reconstruct data after threholding [b, ifail] = c09fd(nwl, c, m, n, f, icomm); % Mean square error of de-noised reconstruction mse = (norm(reshape(a-b,[m*n*f,1]))^2)/(double(m*n*f)); fprintf('With denoising Mean Square Error is %9.6f\n\n',mse); disp('Reconstruction of denoised input: '); disp(b);
c09fz example results Input data a: (:,:,1) = 0.0100 0.0100 0.0100 0.0100 1.0000 1.0000 1.0000 1.0000 0.0100 0.0100 0.0100 0.0100 1.0000 1.0000 1.0000 1.0000 (:,:,2) = 1.0000 1.0000 1.0000 1.0000 0.0100 0.0100 0.0100 0.0100 1.0000 1.0000 1.0000 1.0000 0.0100 0.0100 0.0100 0.0100 (:,:,3) = 0.0100 0.0100 0.0100 0.0100 1.0000 1.0000 1.0000 1.0000 0.0100 0.0100 0.0100 0.0100 1.0000 1.0000 1.0000 1.0000 (:,:,4) = 1.0000 1.0000 1.0000 1.0000 0.0100 0.0100 0.0100 0.0100 1.0000 1.0000 1.0000 1.0000 0.0100 0.0100 0.0100 0.0100 Noisy data an: (:,:,1) = 0.0135 -0.0093 -0.0004 0.0378 1.0015 0.9842 1.0007 0.9889 -0.0017 0.0139 0.0138 -0.0049 0.9899 1.0070 1.0049 0.9983 (:,:,2) = 1.0094 1.0080 0.9921 0.9902 0.0105 -0.0009 0.0160 0.0197 0.9994 1.0044 0.9956 1.0014 0.0091 -0.0084 0.0187 0.0023 (:,:,3) = 0.0058 -0.0053 0.0011 0.0159 1.0113 0.9894 1.0018 0.9992 0.0106 0.0082 0.0093 0.0153 1.0023 1.0157 1.0084 0.9834 (:,:,4) = 0.9969 1.0010 0.9904 0.9968 0.0227 0.0022 0.0062 0.0214 0.9948 0.9981 0.9951 0.9968 0.0121 0.0103 0.0114 0.0206 Without denoising Mean Square Error is 0.000081 Number of coefficients denoised is 55 out of 63 With denoising Mean Square Error is 0.000015 Reconstruction of denoised input: (:,:,1) = 0.0053 0.0053 0.0166 0.0166 1.0026 1.0026 0.9913 0.9913 0.0055 0.0055 0.0077 0.0077 1.0025 1.0025 1.0003 1.0003 (:,:,2) = 1.0026 1.0026 0.9913 0.9913 0.0053 0.0053 0.0166 0.0166 1.0025 1.0025 1.0003 1.0003 0.0055 0.0055 0.0077 0.0077 (:,:,3) = 0.0073 0.0073 0.0110 0.0110 1.0006 1.0006 0.9969 0.9969 0.0078 0.0078 0.0131 0.0131 1.0002 1.0002 0.9949 0.9949 (:,:,4) = 1.0006 1.0006 0.9969 0.9969 0.0073 0.0073 0.0110 0.0110 1.0002 1.0002 0.9949 0.9949 0.0078 0.0078 0.0131 0.0131