hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_wav_3d_init (c09ac)

Purpose

nag_wav_3d_init (c09ac) returns the details of the chosen three-dimensional discrete wavelet filter. For a chosen mother wavelet, discrete wavelet transform type (single-level or multi-level DWT) and end extension method, this function returns the maximum number of levels of resolution (appropriate to a multi-level transform), the filter length, the total number of coefficients and the number of wavelet coefficients in the second and third dimensions for the single-level case. This function must be called before any of the three-dimensional transform functions in this chapter.

Syntax

[nwl, nf, nwct, nwcn, nwcfr, icomm, ifail] = c09ac(wavnam, wtrans, mode, m, n, fr)
[nwl, nf, nwct, nwcn, nwcfr, icomm, ifail] = nag_wav_3d_init(wavnam, wtrans, mode, m, n, fr)

Description

Three-dimensional discrete wavelet transforms (DWT) are characterised by the mother wavelet, the end extension method and whether multiresolution analysis is to be performed. For the selected combination of choices for these three characteristics, and for given dimensions (m × n × frm×n×fr) of data array AA, nag_wav_3d_init (c09ac) returns the dimension details for the transform determined by this combination. The dimension details are: lmaxlmax, the maximum number of levels of resolution that would be computed were a multi-level DWT applied; nfnf, the filter length; nctnct the total number of wavelet coefficients (over all levels in the multi-level DWT case); ncnncn, the number of coefficients in the second dimension for a single-level DWT; and ncfrncfr, the number of coefficients in the third dimension for a single-level DWT. These values are also stored in the communication array icomm, as are the input choices, so that they may be conveniently communicated to the three-dimensional transform functions in this chapter.

References

None.

Parameters

Compulsory Input Parameters

1:     wavnam – string
The name of the mother wavelet. See the C09 Chapter Introduction for details.
wavnam = 'HAAR'wavnam='HAAR'
Haar wavelet.
wavnam = 'DBn'wavnam='DBn', where n = 2,3,,10n=2,3,,10
Daubechies wavelet with nn vanishing moments (2n2n coefficients). For example, wavnam = 'DB4'wavnam='DB4' is the name for the Daubechies wavelet with 44 vanishing moments (88 coefficients).
wavnam = 'BIORxwavnam='BIORx.y'y', where xx.yy can be one of 1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5 or 3.7
Biorthogonal wavelet of order xx.yy. For example wavnam = 'BIOR3.1'wavnam='BIOR3.1' is the name for the biorthogonal wavelet of order 3.13.1.
Constraint: wavnam = 'HAAR'wavnam='HAAR', 'DB2''DB2', 'DB3''DB3', 'DB4''DB4', 'DB5''DB5', 'DB6''DB6', 'DB7''DB7', 'DB8''DB8', 'DB9''DB9', 'DB10''DB10', 'BIOR1.1''BIOR1.1', 'BIOR1.3''BIOR1.3', 'BIOR1.5''BIOR1.5', 'BIOR2.2''BIOR2.2', 'BIOR2.4''BIOR2.4', 'BIOR2.6''BIOR2.6', 'BIOR2.8''BIOR2.8', 'BIOR3.1''BIOR3.1', 'BIOR3.3''BIOR3.3', 'BIOR3.5''BIOR3.5' or 'BIOR3.7''BIOR3.7'.
2:     wtrans – string (length ≥ 1)
The type of discrete wavelet transform that is to be applied.
wtrans = 'S'wtrans='S'
Single-level decomposition or reconstruction by discrete wavelet transform.
wtrans = 'M'wtrans='M'
Multiresolution, by a multi-level DWT or its inverse.
Constraint: wtrans = 'S'wtrans='S' or 'M''M'.
3:     mode – string (length ≥ 1)
The end extension method.
mode = 'P'mode='P'
Periodic end extension.
mode = 'H'mode='H'
Half-point symmetric end extension.
mode = 'W'mode='W'
Whole-point symmetric end extension.
mode = 'Z'mode='Z'
Zero end extension.
Constraint: mode = 'P'mode='P', 'H''H', 'W''W' or 'Z''Z'.
4:     m – int64int32nag_int scalar
The number of elements, mm, in the first dimension (number of rows of each two-dimensional frame) of the input data, AA.
Constraint: m2m2.
5:     n – int64int32nag_int scalar
The number of elements, nn, in the second dimension (number of columns of each two-dimensional frame) of the input data, AA.
Constraint: n2n2.
6:     fr – int64int32nag_int scalar
The number of elements, frfr, in the third dimension (number of frames) of the input data, AA.
Constraint: fr2fr2.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     nwl – int64int32nag_int scalar
The maximum number of levels of resolution, lmaxlmax, that can be computed if a multi-level discrete wavelet transform is applied (wtrans = 'M'wtrans='M'). It is such that 2lmax min (m,n,fr) < 2lmax + 1 2lmax min(m,n,fr) <2lmax+1, for lmaxlmax an integer.
If wtrans = 'S'wtrans='S', nwl is not set.
2:     nf – int64int32nag_int scalar
The filter length, nfnf, for the supplied mother wavelet. This is used to determine the number of coefficients to be generated by the chosen transform.
3:     nwct – int64int32nag_int scalar
The total number of wavelet coefficients, nctnct, that will be generated. When wtrans = 'S'wtrans='S' the number of rows required (i.e., the first dimension of each two-dimensional frame) in each of the output coefficient arrays can be calculated as ncm = nct / (8 × ncn × ncfr)ncm=nct/(8×ncn×ncfr). When wtrans = 'M'wtrans='M' the length of the array used to store all of the coefficient matrices must be at least nctnct.
4:     nwcn – int64int32nag_int scalar
For a single-level transform (wtrans = 'S'wtrans='S'), the number of coefficients that would be generated in the second dimension, ncnncn, for each coefficient type. For a multi-level transform (wtrans = 'M'wtrans='M') this is set to 11.
5:     nwcfr – int64int32nag_int scalar
For a single-level transform (wtrans = 'S'wtrans='S'), the number of coefficients that would be generated in the third dimension, ncfrncfr, for each coefficient type. For a multi-level transform (wtrans = 'M'wtrans='M') this is set to 11.
6:     icomm(260260) – int64int32nag_int array
Contains details of the wavelet transform and the problem dimension which is to be communicated to the two-dimensional discrete transform functions in this chapter.
7:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry, wavnam had an illegal value.
  ifail = 2ifail=2
On entry, wtrans had an illegal value.
  ifail = 3ifail=3
On entry, mode had an illegal value.
  ifail = 4ifail=4
Constraint: fr2fr2.
Constraint: m2m2.
Constraint: n2n2.

Accuracy

Not applicable.

Further Comments

None.

Example

function nag_wav_3d_init_example
m  = int64(8);
n  = int64(8);
fr = int64(8);
wavnam = 'DB4';
mode = 'zero';
wtrans = 'Multilevel';
a = zeros(m, n, fr);
a(:, :, 1) = [10, 31, 04, 10, 13, 15, 04, 06;
              26, 24, 03, 18, 17, 22, 20, 05;
              06, 05, 06, 11, 22, 23, 23, 01;
              09, 15, 18, 01, 30, 24, 08, 01;
              18, 04, 26, 20, 31, 21, 04, 06;
              25, 23, 25, 14, 13, 03, 03, 29;
              22, 29, 07, 29, 13, 31, 03, 12;
              22, 03, 30, 05, 10, 04, 01, 19];
a(:, :, 2) = [01, 02, 14, 31, 19, 28, 06, 15;
              26, 25, 25, 04, 05, 15, 24, 05;
              01, 29, 08, 18, 22, 18, 31, 23;
              08, 04, 16, 21, 14, 02, 02, 21;
              10, 03, 14, 03, 25, 10, 24, 15;
              03, 16, 26, 21, 16, 19, 25, 27;
              28, 29, 01, 20, 03, 24, 31, 28;
              31, 28, 14, 30, 13, 29, 20, 04];
a(:, :, 3) = [31, 26, 23, 05, 22, 01, 16, 08;
              21, 01, 29, 10, 23, 14, 09, 03;
              20, 10, 11, 22, 26, 31, 03, 21;
              09, 24, 19, 03, 04, 01, 13, 29;
              18, 16, 05, 06, 09, 16, 08, 16;
              32, 19, 32, 01, 06, 04, 01, 17;
              29, 29, 02, 29, 27, 25, 31, 06;
              28, 15, 15, 22, 18, 01, 18, 14];
a(:, :, 4) = [15, 09, 04, 14, 26, 10, 03, 28;
              21, 24, 32, 27, 01, 27, 08, 16;
              10, 27, 29, 15, 13, 01, 05, 16;
              04, 01, 08, 31, 14, 06, 05, 27;
              01, 19, 11, 31, 12, 31, 17, 26;
              27, 01, 16, 06, 18, 02, 17, 17;
              30, 09, 15, 32, 32, 29, 16, 02;
              03, 11, 26, 02, 23, 08, 10, 31];
a(:, :, 5) = [12, 07, 06, 12, 01, 13, 30, 26;
              27, 27, 20, 16, 30, 28, 13, 30;
              29, 15, 15, 05, 01, 13, 31, 02;
              31, 21, 27, 30, 08, 07, 11, 03;
              17, 04, 06, 01, 09, 25, 03, 15;
              12, 18, 16, 05, 09, 16, 06, 13;
              03, 05, 26, 30, 19, 11, 32, 24;
              06, 16, 07, 15, 31, 10, 20, 14];
a(:, :, 6) = [20, 07, 17, 11, 04, 21, 25, 17;
              18, 22, 22, 06, 01, 05, 15, 17;
              25, 24, 16, 13, 19, 16, 23, 10;
              01, 31, 05, 13, 11, 12, 01, 18;
              01, 27, 09, 05, 29, 26, 23, 13;
              02, 17, 17, 14, 31, 21, 16, 05;
              26, 21, 10, 21, 09, 11, 01, 15;
              08, 15, 18, 04, 16, 09, 03, 29];
a(:, :, 7) = [26, 02, 30, 26, 07, 04, 09, 01;
              15, 02, 10, 22, 16, 15, 04, 03;
              04, 07, 32, 27, 07, 05, 17, 04;
              22, 30, 06, 18, 32, 02, 01, 31;
              15, 19, 20, 12, 10, 28, 27, 03;
              26, 31, 21, 02, 27, 10, 22, 13;
              32, 03, 27, 23, 01, 11, 04, 26;
              03, 01, 31, 21, 27, 21, 14, 09];
a(:, :, 8) = [02, 16, 16, 23, 23, 09, 27, 12;
              15, 17, 20, 27, 05, 04, 18, 16;
              29, 32, 20, 08, 14, 32, 11, 04;
              28, 01, 15, 19, 14, 09, 30, 18;
              20, 02, 08, 11, 20, 24, 14, 03;
              18, 15, 16, 03, 23, 01, 19, 31;
              32, 27, 28, 09, 15, 23, 09, 13;
              01, 24, 30, 04, 18, 11, 01, 22];


% Query wavelet filter dimensions
[lmax, nf, nwct, nwcn, nwcfr, icomm, ifail] = ...
      nag_wav_3d_init(wavnam, wtrans, mode, m, n, fr);

% Transform one less than the max possible number of levels.
nwl = lmax - 1;

% Perform Discrete Wavelet transform
[c, dwtlvm, dwtlvn, dwtlvfr, icomm, ifail] =  ...
      nag_wav_3d_multi_fwd(n, fr, a, nwct, nwl, icomm);

% nag_wav_3d_init returns nwct based on max levels, so recalculate.
nwct = sum(7*dwtlvm(1:nwl).*dwtlvn(1:nwl).*dwtlvfr(1:nwl)) + ...
       dwtlvm(1)*dwtlvn(1)*dwtlvfr(1);

fprintf(' Number of Levels :                     %10d\n\n', nwl);
fprintf(' Length of wavelet filter :             %10d\n', nf);
fprintf(' Total number of wavelet coefficients : %10d\n\n', nwct);
fprintf(' Number of coefficients in 1st dimension for each level:\n');
fprintf(' %8d', dwtlvm(1:nwl));
fprintf('\n');
fprintf(' Number of coefficients in 2nd dimension for each level:\n');
fprintf(' %8d', dwtlvn(1:nwl));
fprintf('\n');
fprintf(' Number of coefficients in 3rd dimension for each level:\n');
fprintf(' %8d', dwtlvfr(1:nwl));
fprintf('\n');

% Select the deepest level.
want_level = nwl;

% Select the approximation coefficients.
want_coeffs = 0;

% Identify each set of coefficients in c
for ilevel = nwl:-1:1

  if ilevel ~= want_level
    continue
  end

  nwcm = dwtlvm(nwl-ilevel+1);
  nwcn = dwtlvn(nwl-ilevel+1);
  nwcfr = dwtlvfr(nwl-ilevel+1);

  fprintf('\n--------------------------------\n');
  fprintf(' Level %d output is %d by %d by %d.\n', ilevel, nwcm, nwcn, nwcfr);
  fprintf('--------------------------------\n\n');

  for itype_coeffs = 0:7

    if itype_coeffs ~= want_coeffs
      continue
    end

    % Unless we're looking at the deepest level of nesting, which contains
    % approximation coefficients, advance the pointer on past the preceding
    % levels
    if ilevel == nwl
      locc = 0;
    else
      locc = 8*dwtlvm(1)*dwtlvn(1)*dwtlvfr(1);
      for i = ilevel + 1 : nwl - 1
        locc = locc + 7*dwtlvm(nwl-i+1)*dwtlvn(nwl-i+1)*dwtlvfr(nwl-i+1);
      end
    end

    % Now decide which coefficient type we are considering
    switch (itype_coeffs)
      case {0}
        if (ilevel==nwl)
          fprintf('Approximation coefficients (LLL)\n');
          locc = locc + 1;
        end
      case {1}
        fprintf('Detail coefficients (LLH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients
          locc = locc + nwcm*nwcn*nwcfr + 1;
        else
          locc = locc + 1;
        end
      case {2}
        fprintf('Detail coefficients (LHL)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 1 set of
          % detail coefficients
          locc = locc + 2*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 1 set of detail coefficients
          locc = locc + nwcm*nwcn*nwcfr + 1;
        end
      case {3}
        fprintf('Detail coefficients (LHH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 2 sets of
          % detail coefficients
          locc = locc + 3*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 2 sets of detail coefficients
          locc = locc + 2*nwcm*nwcn*nwcfr + 1;
        end
      case {4}
        fprintf('Detail coefficients (HLL)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 3 sets of
          % detail coefficients
          locc = locc + 4*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 3 sets of detail coefficients
          locc = locc + 3*nwcm*nwcn*nwcfr + 1;
        end
      case {5}
        fprintf('Detail coefficients (HLH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 4 sets of
          % detail coefficients
          locc = locc + 5*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 4 sets of detail coefficients
          locc = locc + 4*nwcm*nwcn*nwcfr + 1;
        end
      case {6}
        fprintf('Detail coefficients (HHL)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 5 sets of
          % detail coefficients
          locc = locc + 6*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 4 sets of detail coefficients
          locc = locc + 5*nwcm*nwcn*nwcfr + 1;
        end
      case {7}
        fprintf('Detail coefficients (HHH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 6 sets of
          % detail coefficients
          locc = locc + 7*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 5 sets of detail coefficients
          locc = locc + 6*nwcm*nwcn*nwcfr + 1;
        end
      end

  if itype_coeffs > 0 || ilevel == nwl

    if (itype_coeffs==0)
      % For a multi level transform approx coeffs stored as
      % nwcm x nwcn x nwcfr
      i1 = locc;
      for k = 1:nwcfr
        for j = 1:nwcn
          for i = 1:nwcm
            d(i,j,k) = c(i1);
            i1 = i1 + 1;
          end
        end
      end
    else
      % ... but detail coefficients are stored as ncwfr x nwcm x nwcn
      for k = 1:nwcfr
        for j = 1:nwcn
          for i = 1:nwcm
            i1 = locc - 1 + (j-1)*nwcfr*nwcm + (i-1)*nwcfr + k;
            d(i,j,k) = c(i1);
          end
        end
      end
    end

    % Print out the selected set of coefficients
    fprintf('Level %d, Coefficients %d:\n', ilevel, itype_coeffs);
    format compact;
    for k = 1:nwcfr
      fprintf('Frame %d:\n', k);
      for i = 1:nwcm
        disp(d(i,1:nwcn,k));
      end
    end

  end

  end
end

% Reconstruct original data
[b, ifail] = nag_wav_3d_mxolap_multi_inv(nwl, c, m, n, fr, icomm);

fprintf('\n Reconstruction       b : \n');
% Convert to int16 to get more compact output
disp(int16(b));
 
 Number of Levels :                              2

 Length of wavelet filter :                      8
 Total number of wavelet coefficients :       5145

 Number of coefficients in 1st dimension for each level:
        7        7
 Number of coefficients in 2nd dimension for each level:
        7        7
 Number of coefficients in 3rd dimension for each level:
        7        7

--------------------------------
 Level 2 output is 7 by 7 by 7.
--------------------------------

Approximation coefficients (LLL)
Level 2, Coefficients 0:
Frame 1:
   1.0e-03 *
   -0.0001   -0.0010    0.0024    0.0005    0.1251    0.0186    0.0018
   1.0e-03 *
   -0.0010   -0.0051    0.0376   -0.0781    0.0339   -0.7497   -0.0167
    0.0000    0.0000   -0.0001   -0.0002   -0.0020    0.0036   -0.0002
   -0.0000   -0.0000   -0.0002    0.0021    0.0025   -0.0124    0.0010
    0.0001   -0.0000   -0.0017    0.0009    0.0928    0.1155    0.0004
    0.0002   -0.0007    0.0013   -0.0063    0.1584    0.0931    0.0096
    0.0000   -0.0001    0.0003   -0.0006    0.0123    0.0061    0.0014
Frame 2:
   -0.0000    0.0000    0.0000   -0.0000   -0.0010   -0.0005   -0.0000
    0.0000   -0.0000    0.0001   -0.0006    0.0026    0.0035    0.0004
    0.0001   -0.0000   -0.0008    0.0027    0.0133   -0.0064   -0.0032
   -0.0002    0.0000    0.0032   -0.0067   -0.0708    0.0073    0.0148
   -0.0003    0.0035   -0.0155    0.0406   -0.3676   -0.3434   -0.0682
   -0.0011    0.0004    0.0241   -0.0866   -0.4993   -0.5807   -0.0674
   -0.0002   -0.0003    0.0048   -0.0128   -0.0800   -0.0731   -0.0045
Frame 3:
    0.0000    0.0000   -0.0002    0.0005    0.0006    0.0027    0.0005
   -0.0000    0.0002   -0.0012    0.0037   -0.0224    0.0005   -0.0006
   -0.0002   -0.0011    0.0067   -0.0126    0.0447   -0.0734    0.0068
    0.0008    0.0025   -0.0141   -0.0008    0.0872    0.3261   -0.0494
    0.0012   -0.0173    0.0687   -0.0681    0.5915   -0.1717    0.3943
    0.0016    0.0123   -0.1221    0.4190   -0.5269    1.2295    0.1617
    0.0003    0.0028   -0.0182    0.0396    0.1154    0.2823    0.0102
Frame 4:
   -0.0000   -0.0002    0.0011   -0.0030    0.0059   -0.0102   -0.0026
    0.0000   -0.0010    0.0042   -0.0106    0.0948   -0.0180   -0.0005
    0.0004    0.0061   -0.0296    0.0586   -0.3921    0.3650    0.0134
   -0.0018   -0.0155    0.0684   -0.0636    0.5365   -1.4566    0.0298
   -0.0070    0.0592   -0.1486   -0.1055   -2.9693    0.1109   -1.4193
   -0.0017   -0.0424    0.2595   -0.7280    2.4682   -4.1771   -0.5119
    0.0003   -0.0079    0.0273   -0.0205   -0.1224   -0.9982   -0.0710
Frame 5:
    0.0001   -0.0000   -0.0005   -0.0015    0.0804    0.1009    0.0139
   -0.0006    0.0033   -0.0017   -0.0019   -0.5303   -0.5712   -0.0438
   -0.0014   -0.0157    0.0800   -0.1856    0.4182    0.4931    0.0090
    0.0099    0.0522   -0.4140    1.1260    0.6111   -0.0042   -0.1288
    0.0831   -0.4718    0.9591   -2.9510   84.8494   91.3686   10.1751
    0.1599   -0.3194   -0.8962    1.8546  106.1903  117.2751   12.9904
    0.0213   -0.0211   -0.2179    0.4955   12.5323   12.9746    1.3422
Frame 6:
    0.0002   -0.0004   -0.0006    0.0005    0.0945    0.1342    0.0157
   -0.0008    0.0048   -0.0052    0.0013   -0.7012   -0.3668   -0.0231
   -0.0006   -0.0125    0.0347   -0.0396    1.3945   -0.2227   -0.1395
    0.0034    0.0166   -0.0246   -0.0495   -3.2417   -0.3508    0.3284
    0.1373   -0.4804   -0.1436    0.6068  105.5811  101.7766   10.0719
    0.1359   -0.6132    0.8736   -2.8616  121.1074  124.4215   13.7050
    0.0068   -0.0939    0.4312   -1.4152   12.9366   13.1259    1.6024
Frame 7:
    0.0000   -0.0001    0.0006   -0.0024    0.0134    0.0160    0.0014
   -0.0001    0.0006    0.0003   -0.0044   -0.0813   -0.0377   -0.0021
    0.0006    0.0002   -0.0206    0.0816    0.0851   -0.0274   -0.0148
   -0.0028   -0.0074    0.1035   -0.3488    0.0136   -0.1313    0.0288
    0.0177   -0.0358   -0.0968    0.1416   11.4442   11.6279    0.9779
    0.0187   -0.0759    0.0227    0.1041   13.7268   13.3069    1.5629
    0.0002   -0.0164    0.0748   -0.2042    1.6290    1.2827    0.1547

 Reconstruction       b : 
(:,:,1) =
     10     31      4     10     13     15      4      6
     26     24      3     18     17     22     20      5
      6      5      6     11     22     23     23      1
      9     15     18      1     30     24      8      1
     18      4     26     20     31     21      4      6
     25     23     25     14     13      3      3     29
     22     29      7     29     13     31      3     12
     22      3     30      5     10      4      1     19
(:,:,2) =
      1      2     14     31     19     28      6     15
     26     25     25      4      5     15     24      5
      1     29      8     18     22     18     31     23
      8      4     16     21     14      2      2     21
     10      3     14      3     25     10     24     15
      3     16     26     21     16     19     25     27
     28     29      1     20      3     24     31     28
     31     28     14     30     13     29     20      4
(:,:,3) =
     31     26     23      5     22      1     16      8
     21      1     29     10     23     14      9      3
     20     10     11     22     26     31      3     21
      9     24     19      3      4      1     13     29
     18     16      5      6      9     16      8     16
     32     19     32      1      6      4      1     17
     29     29      2     29     27     25     31      6
     28     15     15     22     18      1     18     14
(:,:,4) =
     15      9      4     14     26     10      3     28
     21     24     32     27      1     27      8     16
     10     27     29     15     13      1      5     16
      4      1      8     31     14      6      5     27
      1     19     11     31     12     31     17     26
     27      1     16      6     18      2     17     17
     30      9     15     32     32     29     16      2
      3     11     26      2     23      8     10     31
(:,:,5) =
     12      7      6     12      1     13     30     26
     27     27     20     16     30     28     13     30
     29     15     15      5      1     13     31      2
     31     21     27     30      8      7     11      3
     17      4      6      1      9     25      3     15
     12     18     16      5      9     16      6     13
      3      5     26     30     19     11     32     24
      6     16      7     15     31     10     20     14
(:,:,6) =
     20      7     17     11      4     21     25     17
     18     22     22      6      1      5     15     17
     25     24     16     13     19     16     23     10
      1     31      5     13     11     12      1     18
      1     27      9      5     29     26     23     13
      2     17     17     14     31     21     16      5
     26     21     10     21      9     11      1     15
      8     15     18      4     16      9      3     29
(:,:,7) =
     26      2     30     26      7      4      9      1
     15      2     10     22     16     15      4      3
      4      7     32     27      7      5     17      4
     22     30      6     18     32      2      1     31
     15     19     20     12     10     28     27      3
     26     31     21      2     27     10     22     13
     32      3     27     23      1     11      4     26
      3      1     31     21     27     21     14      9
(:,:,8) =
      2     16     16     23     23      9     27     12
     15     17     20     27      5      4     18     16
     29     32     20      8     14     32     11      4
     28      1     15     19     14      9     30     18
     20      2      8     11     20     24     14      3
     18     15     16      3     23      1     19     31
     32     27     28      9     15     23      9     13
      1     24     30      4     18     11      1     22

function c09ac_example
m  = int64(8);
n  = int64(8);
fr = int64(8);
wavnam = 'DB4';
mode = 'zero';
wtrans = 'Multilevel';
a = zeros(m, n, fr);
a(:, :, 1) = [10, 31, 04, 10, 13, 15, 04, 06;
              26, 24, 03, 18, 17, 22, 20, 05;
              06, 05, 06, 11, 22, 23, 23, 01;
              09, 15, 18, 01, 30, 24, 08, 01;
              18, 04, 26, 20, 31, 21, 04, 06;
              25, 23, 25, 14, 13, 03, 03, 29;
              22, 29, 07, 29, 13, 31, 03, 12;
              22, 03, 30, 05, 10, 04, 01, 19];
a(:, :, 2) = [01, 02, 14, 31, 19, 28, 06, 15;
              26, 25, 25, 04, 05, 15, 24, 05;
              01, 29, 08, 18, 22, 18, 31, 23;
              08, 04, 16, 21, 14, 02, 02, 21;
              10, 03, 14, 03, 25, 10, 24, 15;
              03, 16, 26, 21, 16, 19, 25, 27;
              28, 29, 01, 20, 03, 24, 31, 28;
              31, 28, 14, 30, 13, 29, 20, 04];
a(:, :, 3) = [31, 26, 23, 05, 22, 01, 16, 08;
              21, 01, 29, 10, 23, 14, 09, 03;
              20, 10, 11, 22, 26, 31, 03, 21;
              09, 24, 19, 03, 04, 01, 13, 29;
              18, 16, 05, 06, 09, 16, 08, 16;
              32, 19, 32, 01, 06, 04, 01, 17;
              29, 29, 02, 29, 27, 25, 31, 06;
              28, 15, 15, 22, 18, 01, 18, 14];
a(:, :, 4) = [15, 09, 04, 14, 26, 10, 03, 28;
              21, 24, 32, 27, 01, 27, 08, 16;
              10, 27, 29, 15, 13, 01, 05, 16;
              04, 01, 08, 31, 14, 06, 05, 27;
              01, 19, 11, 31, 12, 31, 17, 26;
              27, 01, 16, 06, 18, 02, 17, 17;
              30, 09, 15, 32, 32, 29, 16, 02;
              03, 11, 26, 02, 23, 08, 10, 31];
a(:, :, 5) = [12, 07, 06, 12, 01, 13, 30, 26;
              27, 27, 20, 16, 30, 28, 13, 30;
              29, 15, 15, 05, 01, 13, 31, 02;
              31, 21, 27, 30, 08, 07, 11, 03;
              17, 04, 06, 01, 09, 25, 03, 15;
              12, 18, 16, 05, 09, 16, 06, 13;
              03, 05, 26, 30, 19, 11, 32, 24;
              06, 16, 07, 15, 31, 10, 20, 14];
a(:, :, 6) = [20, 07, 17, 11, 04, 21, 25, 17;
              18, 22, 22, 06, 01, 05, 15, 17;
              25, 24, 16, 13, 19, 16, 23, 10;
              01, 31, 05, 13, 11, 12, 01, 18;
              01, 27, 09, 05, 29, 26, 23, 13;
              02, 17, 17, 14, 31, 21, 16, 05;
              26, 21, 10, 21, 09, 11, 01, 15;
              08, 15, 18, 04, 16, 09, 03, 29];
a(:, :, 7) = [26, 02, 30, 26, 07, 04, 09, 01;
              15, 02, 10, 22, 16, 15, 04, 03;
              04, 07, 32, 27, 07, 05, 17, 04;
              22, 30, 06, 18, 32, 02, 01, 31;
              15, 19, 20, 12, 10, 28, 27, 03;
              26, 31, 21, 02, 27, 10, 22, 13;
              32, 03, 27, 23, 01, 11, 04, 26;
              03, 01, 31, 21, 27, 21, 14, 09];
a(:, :, 8) = [02, 16, 16, 23, 23, 09, 27, 12;
              15, 17, 20, 27, 05, 04, 18, 16;
              29, 32, 20, 08, 14, 32, 11, 04;
              28, 01, 15, 19, 14, 09, 30, 18;
              20, 02, 08, 11, 20, 24, 14, 03;
              18, 15, 16, 03, 23, 01, 19, 31;
              32, 27, 28, 09, 15, 23, 09, 13;
              01, 24, 30, 04, 18, 11, 01, 22];


% Query wavelet filter dimensions
[lmax, nf, nwct, nwcn, nwcfr, icomm, ifail] = ...
      c09ac(wavnam, wtrans, mode, m, n, fr);

% Transform one less than the max possible number of levels.
nwl = lmax - 1;

% Perform Discrete Wavelet transform
[c, dwtlvm, dwtlvn, dwtlvfr, icomm, ifail] = c09fc(n, fr, a, nwct, nwl, icomm);

% c09ac returns nwct based on max levels, so recalculate.
nwct = sum(7*dwtlvm(1:nwl).*dwtlvn(1:nwl).*dwtlvfr(1:nwl)) + ...
       dwtlvm(1)*dwtlvn(1)*dwtlvfr(1);

fprintf(' Number of Levels :                     %10d\n\n', nwl);
fprintf(' Length of wavelet filter :             %10d\n', nf);
fprintf(' Total number of wavelet coefficients : %10d\n\n', nwct);
fprintf(' Number of coefficients in 1st dimension for each level:\n');
fprintf(' %8d', dwtlvm(1:nwl));
fprintf('\n');
fprintf(' Number of coefficients in 2nd dimension for each level:\n');
fprintf(' %8d', dwtlvn(1:nwl));
fprintf('\n');
fprintf(' Number of coefficients in 3rd dimension for each level:\n');
fprintf(' %8d', dwtlvfr(1:nwl));
fprintf('\n');

% Select the deepest level.
want_level = nwl;

% Select the approximation coefficients.
want_coeffs = 0;

% Identify each set of coefficients in c
for ilevel = nwl:-1:1

  if ilevel ~= want_level
    continue
  end

  nwcm = dwtlvm(nwl-ilevel+1);
  nwcn = dwtlvn(nwl-ilevel+1);
  nwcfr = dwtlvfr(nwl-ilevel+1);

  fprintf('\n--------------------------------\n');
  fprintf(' Level %d output is %d by %d by %d.\n', ilevel, nwcm, nwcn, nwcfr);
  fprintf('--------------------------------\n\n');

  for itype_coeffs = 0:7

    if itype_coeffs ~= want_coeffs
      continue
    end

    % Unless we're looking at the deepest level of nesting, which contains
    % approximation coefficients, advance the pointer on past the preceding
    % levels
    if ilevel == nwl
      locc = 0;
    else
      locc = 8*dwtlvm(1)*dwtlvn(1)*dwtlvfr(1);
      for i = ilevel + 1 : nwl - 1
        locc = locc + 7*dwtlvm(nwl-i+1)*dwtlvn(nwl-i+1)*dwtlvfr(nwl-i+1);
      end
    end

    % Now decide which coefficient type we are considering
    switch (itype_coeffs)
      case {0}
        if (ilevel==nwl)
          fprintf('Approximation coefficients (LLL)\n');
          locc = locc + 1;
        end
      case {1}
        fprintf('Detail coefficients (LLH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients
          locc = locc + nwcm*nwcn*nwcfr + 1;
        else
          locc = locc + 1;
        end
      case {2}
        fprintf('Detail coefficients (LHL)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 1 set of
          % detail coefficients
          locc = locc + 2*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 1 set of detail coefficients
          locc = locc + nwcm*nwcn*nwcfr + 1;
        end
      case {3}
        fprintf('Detail coefficients (LHH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 2 sets of
          % detail coefficients
          locc = locc + 3*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 2 sets of detail coefficients
          locc = locc + 2*nwcm*nwcn*nwcfr + 1;
        end
      case {4}
        fprintf('Detail coefficients (HLL)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 3 sets of
          % detail coefficients
          locc = locc + 4*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 3 sets of detail coefficients
          locc = locc + 3*nwcm*nwcn*nwcfr + 1;
        end
      case {5}
        fprintf('Detail coefficients (HLH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 4 sets of
          % detail coefficients
          locc = locc + 5*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 4 sets of detail coefficients
          locc = locc + 4*nwcm*nwcn*nwcfr + 1;
        end
      case {6}
        fprintf('Detail coefficients (HHL)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 5 sets of
          % detail coefficients
          locc = locc + 6*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 4 sets of detail coefficients
          locc = locc + 5*nwcm*nwcn*nwcfr + 1;
        end
      case {7}
        fprintf('Detail coefficients (HHH)\n');
        if (ilevel==nwl)
          % Advance pointer past approximation coefficients and 6 sets of
          % detail coefficients
          locc = locc + 7*nwcm*nwcn*nwcfr + 1;
        else
          % Advance pointer past 5 sets of detail coefficients
          locc = locc + 6*nwcm*nwcn*nwcfr + 1;
        end
      end

  if itype_coeffs > 0 || ilevel == nwl

    if (itype_coeffs==0)
      % For a multi level transform approx coeffs stored as
      % nwcm x nwcn x nwcfr
      i1 = locc;
      for k = 1:nwcfr
        for j = 1:nwcn
          for i = 1:nwcm
            d(i,j,k) = c(i1);
            i1 = i1 + 1;
          end
        end
      end
    else
      % ... but detail coefficients are stored as ncwfr x nwcm x nwcn
      for k = 1:nwcfr
        for j = 1:nwcn
          for i = 1:nwcm
            i1 = locc - 1 + (j-1)*nwcfr*nwcm + (i-1)*nwcfr + k;
            d(i,j,k) = c(i1);
          end
        end
      end
    end

    % Print out the selected set of coefficients
    fprintf('Level %d, Coefficients %d:\n', ilevel, itype_coeffs);
    format compact;
    for k = 1:nwcfr
      fprintf('Frame %d:\n', k);
      for i = 1:nwcm
        disp(d(i,1:nwcn,k));
      end
    end

  end

  end
end

% Reconstruct original data
[b, ifail] = c09fd(nwl, c, m, n, fr, icomm);

fprintf('\n Reconstruction       b : \n');
% Convert to int16 to get more compact output
disp(int16(b));
 
 Number of Levels :                              2

 Length of wavelet filter :                      8
 Total number of wavelet coefficients :       5145

 Number of coefficients in 1st dimension for each level:
        7        7
 Number of coefficients in 2nd dimension for each level:
        7        7
 Number of coefficients in 3rd dimension for each level:
        7        7

--------------------------------
 Level 2 output is 7 by 7 by 7.
--------------------------------

Approximation coefficients (LLL)
Level 2, Coefficients 0:
Frame 1:
   1.0e-03 *
   -0.0001   -0.0010    0.0024    0.0005    0.1251    0.0186    0.0018
   1.0e-03 *
   -0.0010   -0.0051    0.0376   -0.0781    0.0339   -0.7497   -0.0167
    0.0000    0.0000   -0.0001   -0.0002   -0.0020    0.0036   -0.0002
   -0.0000   -0.0000   -0.0002    0.0021    0.0025   -0.0124    0.0010
    0.0001   -0.0000   -0.0017    0.0009    0.0928    0.1155    0.0004
    0.0002   -0.0007    0.0013   -0.0063    0.1584    0.0931    0.0096
    0.0000   -0.0001    0.0003   -0.0006    0.0123    0.0061    0.0014
Frame 2:
   -0.0000    0.0000    0.0000   -0.0000   -0.0010   -0.0005   -0.0000
    0.0000   -0.0000    0.0001   -0.0006    0.0026    0.0035    0.0004
    0.0001   -0.0000   -0.0008    0.0027    0.0133   -0.0064   -0.0032
   -0.0002    0.0000    0.0032   -0.0067   -0.0708    0.0073    0.0148
   -0.0003    0.0035   -0.0155    0.0406   -0.3676   -0.3434   -0.0682
   -0.0011    0.0004    0.0241   -0.0866   -0.4993   -0.5807   -0.0674
   -0.0002   -0.0003    0.0048   -0.0128   -0.0800   -0.0731   -0.0045
Frame 3:
    0.0000    0.0000   -0.0002    0.0005    0.0006    0.0027    0.0005
   -0.0000    0.0002   -0.0012    0.0037   -0.0224    0.0005   -0.0006
   -0.0002   -0.0011    0.0067   -0.0126    0.0447   -0.0734    0.0068
    0.0008    0.0025   -0.0141   -0.0008    0.0872    0.3261   -0.0494
    0.0012   -0.0173    0.0687   -0.0681    0.5915   -0.1717    0.3943
    0.0016    0.0123   -0.1221    0.4190   -0.5269    1.2295    0.1617
    0.0003    0.0028   -0.0182    0.0396    0.1154    0.2823    0.0102
Frame 4:
   -0.0000   -0.0002    0.0011   -0.0030    0.0059   -0.0102   -0.0026
    0.0000   -0.0010    0.0042   -0.0106    0.0948   -0.0180   -0.0005
    0.0004    0.0061   -0.0296    0.0586   -0.3921    0.3650    0.0134
   -0.0018   -0.0155    0.0684   -0.0636    0.5365   -1.4566    0.0298
   -0.0070    0.0592   -0.1486   -0.1055   -2.9693    0.1109   -1.4193
   -0.0017   -0.0424    0.2595   -0.7280    2.4682   -4.1771   -0.5119
    0.0003   -0.0079    0.0273   -0.0205   -0.1224   -0.9982   -0.0710
Frame 5:
    0.0001   -0.0000   -0.0005   -0.0015    0.0804    0.1009    0.0139
   -0.0006    0.0033   -0.0017   -0.0019   -0.5303   -0.5712   -0.0438
   -0.0014   -0.0157    0.0800   -0.1856    0.4182    0.4931    0.0090
    0.0099    0.0522   -0.4140    1.1260    0.6111   -0.0042   -0.1288
    0.0831   -0.4718    0.9591   -2.9510   84.8494   91.3686   10.1751
    0.1599   -0.3194   -0.8962    1.8546  106.1903  117.2751   12.9904
    0.0213   -0.0211   -0.2179    0.4955   12.5323   12.9746    1.3422
Frame 6:
    0.0002   -0.0004   -0.0006    0.0005    0.0945    0.1342    0.0157
   -0.0008    0.0048   -0.0052    0.0013   -0.7012   -0.3668   -0.0231
   -0.0006   -0.0125    0.0347   -0.0396    1.3945   -0.2227   -0.1395
    0.0034    0.0166   -0.0246   -0.0495   -3.2417   -0.3508    0.3284
    0.1373   -0.4804   -0.1436    0.6068  105.5811  101.7766   10.0719
    0.1359   -0.6132    0.8736   -2.8616  121.1074  124.4215   13.7050
    0.0068   -0.0939    0.4312   -1.4152   12.9366   13.1259    1.6024
Frame 7:
    0.0000   -0.0001    0.0006   -0.0024    0.0134    0.0160    0.0014
   -0.0001    0.0006    0.0003   -0.0044   -0.0813   -0.0377   -0.0021
    0.0006    0.0002   -0.0206    0.0816    0.0851   -0.0274   -0.0148
   -0.0028   -0.0074    0.1035   -0.3488    0.0136   -0.1313    0.0288
    0.0177   -0.0358   -0.0968    0.1416   11.4442   11.6279    0.9779
    0.0187   -0.0759    0.0227    0.1041   13.7268   13.3069    1.5629
    0.0002   -0.0164    0.0748   -0.2042    1.6290    1.2827    0.1547

 Reconstruction       b : 
(:,:,1) =
     10     31      4     10     13     15      4      6
     26     24      3     18     17     22     20      5
      6      5      6     11     22     23     23      1
      9     15     18      1     30     24      8      1
     18      4     26     20     31     21      4      6
     25     23     25     14     13      3      3     29
     22     29      7     29     13     31      3     12
     22      3     30      5     10      4      1     19
(:,:,2) =
      1      2     14     31     19     28      6     15
     26     25     25      4      5     15     24      5
      1     29      8     18     22     18     31     23
      8      4     16     21     14      2      2     21
     10      3     14      3     25     10     24     15
      3     16     26     21     16     19     25     27
     28     29      1     20      3     24     31     28
     31     28     14     30     13     29     20      4
(:,:,3) =
     31     26     23      5     22      1     16      8
     21      1     29     10     23     14      9      3
     20     10     11     22     26     31      3     21
      9     24     19      3      4      1     13     29
     18     16      5      6      9     16      8     16
     32     19     32      1      6      4      1     17
     29     29      2     29     27     25     31      6
     28     15     15     22     18      1     18     14
(:,:,4) =
     15      9      4     14     26     10      3     28
     21     24     32     27      1     27      8     16
     10     27     29     15     13      1      5     16
      4      1      8     31     14      6      5     27
      1     19     11     31     12     31     17     26
     27      1     16      6     18      2     17     17
     30      9     15     32     32     29     16      2
      3     11     26      2     23      8     10     31
(:,:,5) =
     12      7      6     12      1     13     30     26
     27     27     20     16     30     28     13     30
     29     15     15      5      1     13     31      2
     31     21     27     30      8      7     11      3
     17      4      6      1      9     25      3     15
     12     18     16      5      9     16      6     13
      3      5     26     30     19     11     32     24
      6     16      7     15     31     10     20     14
(:,:,6) =
     20      7     17     11      4     21     25     17
     18     22     22      6      1      5     15     17
     25     24     16     13     19     16     23     10
      1     31      5     13     11     12      1     18
      1     27      9      5     29     26     23     13
      2     17     17     14     31     21     16      5
     26     21     10     21      9     11      1     15
      8     15     18      4     16      9      3     29
(:,:,7) =
     26      2     30     26      7      4      9      1
     15      2     10     22     16     15      4      3
      4      7     32     27      7      5     17      4
     22     30      6     18     32      2      1     31
     15     19     20     12     10     28     27      3
     26     31     21      2     27     10     22     13
     32      3     27     23      1     11      4     26
      3      1     31     21     27     21     14      9
(:,:,8) =
      2     16     16     23     23      9     27     12
     15     17     20     27      5      4     18     16
     29     32     20      8     14     32     11      4
     28      1     15     19     14      9     30     18
     20      2      8     11     20     24     14      3
     18     15     16      3     23      1     19     31
     32     27     28      9     15     23      9     13
      1     24     30      4     18     11      1     22


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013