MATRIX_IWT Discrete Inverse N-D Wavelet Transform.
IWTND(WX,RH,RG,SCALES) calculates the N-D inverse wavelet transform of
vector WX, which should be a SCALES-scales direct wavelet transform. By
default, N, the number of dimensions to transform is the same as the
number of the dimensions of the input WX. The second argument RH is the
synthesis lowpass filter and the third argument RG the synthesis highpass
filter.
Notes 1: Dimensions to transform
--------------------------------
If the passed matrix WX has N dimensions, then the default behaviour is to
transform all N dimensions. It is possible to select a subset of
dimensions to transform, using the P_DIMS input (see below). For example
we could choose to do a 2D transform on a 3D matrix.
Notes 2: Output dimensions from the iwt
---------------------------------------
The original wavelet transform may have added rows columns etc to the
output, because it needed dimensions divisible by 2 to do the
transform. In doing the inverse transform, we need to know or work out
what the original dimensions were. IWTND will calculate the size of the
reconstructed matrix dimensions to be as large as possible (maybe 1 point
larger than the original) unless the size is given using
IWTND(WX,RH,RG,SCALES,SIZ). A value of 0 for SIZ is the same as omitting
it. SIZ should either be scalar, in which case it it the default size for
any reconstructed dimension, or should be a vector with the same number of
values as reconstructed dimensions (see above, and P_DIM); each value
specifies an output dimension size.
Note that iwtnd differs a little from UviWave iwt2d in its understanding
of dimensions. For iwtnd, the first value in the SIZ vector refers to the
number of rows (M) and the second value refers to the number of columns
(N). UviWave iwt2d (see help iwt2d) takes the first passed size (SIZX) as
the width (number of columns N) and the second passed size (SIZY) as the
height (number of rows, M). In this sense iwtnd follows the matlab
convention (number of rows first, then number of columns).
IWTND can be used to perform a single process of multiresolution
analysis. The way to do it is by selecting the scales whose highpass bands
(detail signals) should be ignored for reconstruction.
This can be done Using IWTND(WX,RH,RG,SCALES,SIZ,SC_LEVELS) where
SC_LEVELS is a SCALES-sized vector,1's or 0's. An i-th coefficient of 0
means that the i-th scale detail (starting from the deepest) should be
ignored. SC_LEVELS vector can be replaced by a single number for selecting
just only the SC_LEVELS deepest scales.
An all-ones vector, or an empty values, or a single number equal to
SCALES, is the same as the normal inverse transform.
IWTND(WX,RH,RG,SCALES,SIZ,SC_LEVELS,DEL1,DEL2) calculates the inverse
transform or performs the multiresolution analysis, but allowing the user
to change the alignment of the outputs with respect to the input
signal. This effect is achieved by setting to DEL1 and DEL2 the analysis
delays of H and G respectively, and calculating the complementary delays
for synthesis filters RH and RG. The default values of DEL1 and DEL2 are
calculated using the function WTCENTER.
Note 3: Using P_DIMS
--------------------
IWTND(WX,RH,RG,SCALES,SIZ,SC_LEVELS,DEL1,DEL2,P_DIMS) adds the ability to
iwt over a subset of the matrix dimensions. P_DIMS should be a matrix of
flag values, one value per dimension of X; each dimension with a flag
value of 1 is iwt'ed, dimensions of size 1 or with a 0 flag value are
ignored. Thus, if we have a 3D matrix, the default value for P_DIMS will
be [1 1 1]; if you want to transform just over the rows of this matrix,
you could override the default by passing a value of P_DIMS of [0 1 0].
P_DIMS can also be a scalar, in which case it specifies that the first
P_DIMS dimensions should be processed; for example, for a 3D matrix, a
P_DIMS value of 2 will expand to a P_DIMS vector of [1 1 0].
Note 4: Different dimensions, different filters
-----------------------------------------------
Each dimension can have its own wavelet filter; the input filters rh and
rg can be given as cell arrays, each cell containing a filter for the
matching processed dimension. For example, let use imagine you are passing
a 3D matrix, but only want to transform the first 2 dimensions. You want
a different filter for the first and second dimensions. You would use a
P_DIMS of [1 1 0] (see above), and pass a two element cell array in h and
in g, with the rh filter for dimension 2 in h{2}, etc. Similarly the
delays DEL1 and DEL2 can be vectors, with one delay for each processed
dimension.
See also: WT, WTND, WTCENTER, WTMETHOD
Based on iwt.m from UviWave 3.0, with thanks - see below
$Id: iwtnd.m,v 1.1 2004/09/26 04:00:24 matthewbrett Exp $0001 function y = iwtnd(wx, rh, rg, scales, o_sz, sc_levels, del1, del2, p_dims) 0002 % MATRIX_IWT Discrete Inverse N-D Wavelet Transform. 0003 % 0004 % IWTND(WX,RH,RG,SCALES) calculates the N-D inverse wavelet transform of 0005 % vector WX, which should be a SCALES-scales direct wavelet transform. By 0006 % default, N, the number of dimensions to transform is the same as the 0007 % number of the dimensions of the input WX. The second argument RH is the 0008 % synthesis lowpass filter and the third argument RG the synthesis highpass 0009 % filter. 0010 % 0011 % Notes 1: Dimensions to transform 0012 % -------------------------------- 0013 % 0014 % If the passed matrix WX has N dimensions, then the default behaviour is to 0015 % transform all N dimensions. It is possible to select a subset of 0016 % dimensions to transform, using the P_DIMS input (see below). For example 0017 % we could choose to do a 2D transform on a 3D matrix. 0018 % 0019 % Notes 2: Output dimensions from the iwt 0020 % --------------------------------------- 0021 % 0022 % The original wavelet transform may have added rows columns etc to the 0023 % output, because it needed dimensions divisible by 2 to do the 0024 % transform. In doing the inverse transform, we need to know or work out 0025 % what the original dimensions were. IWTND will calculate the size of the 0026 % reconstructed matrix dimensions to be as large as possible (maybe 1 point 0027 % larger than the original) unless the size is given using 0028 % IWTND(WX,RH,RG,SCALES,SIZ). A value of 0 for SIZ is the same as omitting 0029 % it. SIZ should either be scalar, in which case it it the default size for 0030 % any reconstructed dimension, or should be a vector with the same number of 0031 % values as reconstructed dimensions (see above, and P_DIM); each value 0032 % specifies an output dimension size. 0033 % 0034 % Note that iwtnd differs a little from UviWave iwt2d in its understanding 0035 % of dimensions. For iwtnd, the first value in the SIZ vector refers to the 0036 % number of rows (M) and the second value refers to the number of columns 0037 % (N). UviWave iwt2d (see help iwt2d) takes the first passed size (SIZX) as 0038 % the width (number of columns N) and the second passed size (SIZY) as the 0039 % height (number of rows, M). In this sense iwtnd follows the matlab 0040 % convention (number of rows first, then number of columns). 0041 % 0042 % IWTND can be used to perform a single process of multiresolution 0043 % analysis. The way to do it is by selecting the scales whose highpass bands 0044 % (detail signals) should be ignored for reconstruction. 0045 % 0046 % This can be done Using IWTND(WX,RH,RG,SCALES,SIZ,SC_LEVELS) where 0047 % SC_LEVELS is a SCALES-sized vector,1's or 0's. An i-th coefficient of 0 0048 % means that the i-th scale detail (starting from the deepest) should be 0049 % ignored. SC_LEVELS vector can be replaced by a single number for selecting 0050 % just only the SC_LEVELS deepest scales. 0051 % 0052 % An all-ones vector, or an empty values, or a single number equal to 0053 % SCALES, is the same as the normal inverse transform. 0054 % 0055 % IWTND(WX,RH,RG,SCALES,SIZ,SC_LEVELS,DEL1,DEL2) calculates the inverse 0056 % transform or performs the multiresolution analysis, but allowing the user 0057 % to change the alignment of the outputs with respect to the input 0058 % signal. This effect is achieved by setting to DEL1 and DEL2 the analysis 0059 % delays of H and G respectively, and calculating the complementary delays 0060 % for synthesis filters RH and RG. The default values of DEL1 and DEL2 are 0061 % calculated using the function WTCENTER. 0062 % 0063 % Note 3: Using P_DIMS 0064 % -------------------- 0065 % 0066 % IWTND(WX,RH,RG,SCALES,SIZ,SC_LEVELS,DEL1,DEL2,P_DIMS) adds the ability to 0067 % iwt over a subset of the matrix dimensions. P_DIMS should be a matrix of 0068 % flag values, one value per dimension of X; each dimension with a flag 0069 % value of 1 is iwt'ed, dimensions of size 1 or with a 0 flag value are 0070 % ignored. Thus, if we have a 3D matrix, the default value for P_DIMS will 0071 % be [1 1 1]; if you want to transform just over the rows of this matrix, 0072 % you could override the default by passing a value of P_DIMS of [0 1 0]. 0073 % P_DIMS can also be a scalar, in which case it specifies that the first 0074 % P_DIMS dimensions should be processed; for example, for a 3D matrix, a 0075 % P_DIMS value of 2 will expand to a P_DIMS vector of [1 1 0]. 0076 % 0077 % Note 4: Different dimensions, different filters 0078 % ----------------------------------------------- 0079 % 0080 % Each dimension can have its own wavelet filter; the input filters rh and 0081 % rg can be given as cell arrays, each cell containing a filter for the 0082 % matching processed dimension. For example, let use imagine you are passing 0083 % a 3D matrix, but only want to transform the first 2 dimensions. You want 0084 % a different filter for the first and second dimensions. You would use a 0085 % P_DIMS of [1 1 0] (see above), and pass a two element cell array in h and 0086 % in g, with the rh filter for dimension 2 in h{2}, etc. Similarly the 0087 % delays DEL1 and DEL2 can be vectors, with one delay for each processed 0088 % dimension. 0089 % 0090 % See also: WT, WTND, WTCENTER, WTMETHOD 0091 % 0092 % Based on iwt.m from UviWave 3.0, with thanks - see below 0093 % 0094 % $Id: iwtnd.m,v 1.1 2004/09/26 04:00:24 matthewbrett Exp $ 0095 0096 % Restrictions: 0097 % 0098 % - Synthesis filters from the same set as in analysis must be used. If 0099 % forced delays were used in analysis, the same delays should be forced in 0100 % synthesis (iwtnd calculates the complementary ones), so as to get a 0101 % perfect vector reconstruction (or the equivalent in the case of not full 0102 % reconstruction). 0103 % 0104 % - The number of scales indicated in SCALES, must match the number of 0105 % scales specified in the analysis process. Otherwise, the reconstruction 0106 % will be completely incorrect. The same applies to the original vector size 0107 % SIZ, but if it's not given, or it's set to zero, IWTND will give the 0108 % reconstructed vector the largest size possible. 0109 0110 %-------------------------------------------------------- 0111 % iwt: Copyright (C) 1994, 1995, 1996, by Universidad de Vigo 0112 % 0113 % Uvi_Wave is free software; you can redistribute it and/or modify it 0114 % under the terms of the GNU General Public License as published by the 0115 % Free Software Foundation; either version 2, or (at your option) any 0116 % later version. 0117 % 0118 % Uvi_Wave is distributed in the hope that it will be useful, but WITHOUT 0119 % ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 0120 % FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 0121 % for more details. 0122 % 0123 % You should have received a copy of the GNU General Public License 0124 % along with Uvi_Wave; see the file COPYING. If not, write to the Free 0125 % Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 0126 % 0127 % Authors: Sergio J. Garcia Galan 0128 % Cristina Sanchez Cabanelas 0129 % e-mail: Uvi_Wave@tsc.uvigo.es 0130 %-------------------------------------------------------- 0131 0132 % ----------------------------------- 0133 % CHECK PARAMETERS AND OPTIONS 0134 % ----------------------------------- 0135 0136 if nargin < 4 0137 error('Need at least matrix to iwt, iwt filters and scales'); 0138 end 0139 if nargin < 5 0140 o_sz = []; 0141 end 0142 if nargin < 6, 0143 sc_levels=[]; 0144 end; 0145 if nargin < 7 0146 del1 = []; 0147 end 0148 if nargin < 8 0149 del2 = []; 0150 end 0151 if xor(isempty(del1), isempty(del2)) 0152 error('You cannot specify only one of the two delays'); 0153 end 0154 if nargin < 9 0155 p_dims = []; 0156 end 0157 0158 % Process scales levels: if SC_LEVELS specifies the number of scales then 0159 % build the SC_LEVELS vector with SC_LEVELS ones and SCALES-SC_LEVELS zeros. 0160 sc_l_l = prod(size(sc_levels)); 0161 if sc_l_l == 0 % empty 0162 % SC_LEVELS not given means reconstructing all bands. 0163 sc_levels=ones(1, scales); 0164 elseif sc_l_l == 1 0165 sc_levels=[ones(1, sc_levels) ... 0166 zeros(1, scales - sc_levels)]; 0167 elseif sc_l_l ~= scales 0168 error(['SC_LEVELS should be a single number (<=SCALES) ' ... 0169 'or a vector with SCALES elements']); 0170 else 0171 % And make sure that all nonzero elements in SC_LEVELS are ones 0172 sc_levels(sc_levels~=0) = 1; 0173 end 0174 0175 % input dimensions 0176 sz = size(wx); 0177 n_dims = length(sz); 0178 0179 % Process p_dims 0180 p_d_l = prod(size(p_dims)); 0181 if p_d_l == 0 % empty 0182 p_dims = ones(1, n_dims); 0183 elseif p_d_l == 1 0184 p_dims = [ones(1, p_dims) zeros(1, n_dims-p_dims)]; 0185 elseif p_d_l < n_dims 0186 p_dims = [p_dims(:)' ones(1, n_dims-p_d_l)]; 0187 elseif p_d_l > n_dims 0188 error('p_dims array has too many values for passed matrix'); 0189 end 0190 p_dims = p_dims & sz > 1; 0191 if ~any(p_dims), y = wx; return, end 0192 n_p_dims = sum(p_dims); 0193 0194 % make filters into cells, and expand to N-D if necessary 0195 if ~iscell(rh), rh = {rh}; end 0196 if prod(size(rh))==1, rh = repmat(rh, 1, n_p_dims); end 0197 if ~iscell(rg), rg = {rg}; end 0198 if prod(size(rg))==1, rg = repmat(rg, 1, n_p_dims); end 0199 0200 % ---------------------------------- 0201 % CHECK THE ORIGINAL LENGTH 0202 % ---------------------------------- 0203 0204 % Process input sizes 0205 def_o_sz = zeros(1, n_dims); 0206 def_o_sz(~p_dims) = sz(~p_dims); 0207 switch prod(size(o_sz)); 0208 case 0 % empty 0209 o_sz = def_o_sz; 0210 case {1, n_p_dims} 0211 def_o_sz(p_dims) = o_sz(:)'; 0212 o_sz = def_o_sz; 0213 case n_dims 0214 % Dimensions to be transformed will not change 0215 tmp = ~p_dims & o_sz==0; 0216 o_sz(tmp) = sz(tmp); 0217 if any(o_sz(~p_dims) ~= sz(~p_dims)) 0218 error('Dimensions not being transformed should not change size'); 0219 end 0220 otherwise 0221 error(['Original size vector should have one element for each ' ... 0222 'dimension of the input matrix, or one element for each '... 0223 'processed dimension, or be a scalar']); 0224 end 0225 0226 % Get maximum size for any missing input sizes 0227 missing_dims = (o_sz == 0); 0228 for d = find(o_sz==0) 0229 o_sz(d) = maxrsize(sz(d), scales); 0230 end 0231 if any(o_sz == 0) 0232 error(sprintf(['Can''t determine the original length for' ... 0233 'dimension %d; SCALES might be wrong\n'], ... 0234 find(o_sz == 0))); 0235 end 0236 0237 % Estimated input / output sizes for current original sizes 0238 [wt_sz sc_o_sz sc_wt_sz] = wt_dims(o_sz, scales, p_dims); 0239 0240 % Check as yet unchecked dims to see if they work 0241 dims_to_check = p_dims & ~missing_dims; 0242 dims_to_redo = find(wt_sz(dims_to_check) ~= sz(dims_to_check)); 0243 if ~isempty(dims_to_redo) 0244 for r = 1:length(dims_to_redo) 0245 new_o_sz(r) = maxrsize(sz(dims_to_redo(r)), scales); 0246 end 0247 fprintf('Size %d for dim %d is not correct, trying default %d: \n', ... 0248 [o_sz(dims_to_redo); dims_to_redo; new_o_sz]); 0249 if any(new_o_sz == 0) 0250 error(sprintf(['No default size found for dim %d. ' ... 0251 'SCALES might be wrong\n'], ... 0252 dims_to_redo(new_o_sz==0))); 0253 end 0254 o_sz(dims_to_redo) = new_o_sz; 0255 [wt_sz sc_o_sz sc_wt_sz] = wt_dims(o_sz, scales, p_dims); 0256 0257 wrong_sz = wt_sz ~= sz; 0258 if any(wrong_sz) 0259 error(sprintf(['Default failed for dim %d. ' ... 0260 'SCALES might be wrong\n'], ... 0261 find(wrong_sz))); 0262 end 0263 fprintf(['(Check the resulting size(s); may be wrong. ' ... 0264 'If so, check SCALES)\n']); 0265 end 0266 0267 %-------------------------- 0268 % DELAY CALCULATION 0269 %-------------------------- 0270 0271 for d = 1:n_p_dims 0272 llp(d) = length(rh{d}); % Length of the lowpass filters. 0273 lhp(d) = length(rg{d}); % Length of the highpass filters. 0274 end 0275 0276 % The total delay of the analysis-synthesis process must match the sum of 0277 % the analysis delay plus the synthesis delay. SUML holds this total delay, 0278 % which is different depending on the kind of filters. 0279 suml = llp+lhp-2; 0280 difl = abs(lhp-llp); 0281 tmp = rem(difl,2)==0; 0282 suml(tmp) = suml(tmp)/2; 0283 0284 if isempty(del1) 0285 % Calculate analysis delays as the reciprocal M. C. 0286 for d = 1:n_p_dims 0287 dlpa(d) = wtcenter( rg{d} ); 0288 dhpa(d) = wtcenter( rh{d} ); 0289 end 0290 % difference between them must be even 0291 dhpa = dhpa + (rem(dhpa-dlpa,2)~=0); 0292 else 0293 % Other experimental filter delays can be forced from the arguments 0294 o = ones(1, n_p_dims); 0295 dlpa=del1; if prod(size(dlpa))==1, dlpa = o * dlpa; end 0296 dhpa=del2; if prod(size(dhpa))==1, dhpa = o * dhpa; end 0297 end 0298 0299 % Synthesis delays are the total minus the analysis delays. 0300 dlp = suml - dlpa; 0301 dhp = suml - dhpa; 0302 0303 %------------------------------ 0304 % START THE ALGORITHM 0305 %------------------------------ 0306 0307 % For doing tricksy permute thing in loop below (see help wt_dim_shifts) 0308 [shifts end_shift] = wt_dim_shifts(p_dims, sz); 0309 dims = [1:n_dims]'; 0310 0311 % upend scale size matrix, scales has opposite meaning compared to wt 0312 sc_wt_sz = flipud(sc_wt_sz); 0313 0314 % Flags for truncation (scales by dims) 0315 truncations = sc_wt_sz - flipud(sc_o_sz); 0316 0317 % offsets caused by extra 0's put in by wt 0318 offsets = zeros(1, n_dims); 0319 0320 for sc=1:scales 0321 0322 % get new data block containing lowpass + highpass bits of vector 0323 % following assumes UviWave ordering 0324 0325 % do simplest case simply for speed 0326 if sc==scales & ~any(offsets) 0327 t = wx; 0328 else % general case 0329 in_sz = sc_wt_sz(sc,:); 0330 op1 = offsets + 1; 0331 switch n_dims 0332 case 2 0333 t = wx(op1(1):in_sz(1)+offsets(1), ... 0334 op1(2):in_sz(2)+offsets(2)); 0335 case 3 0336 t = wx(op1(1):in_sz(1)+offsets(1), ... 0337 op1(2):in_sz(2)+offsets(2), ... 0338 op1(3):in_sz(3)+offsets(3)); 0339 case 4 0340 t = wx(op1(1):in_sz(1)+offsets(1), ... 0341 op1(2):in_sz(2)+offsets(2), ... 0342 op1(3):in_sz(3)+offsets(3), ... 0343 op1(4):in_sz(4)+offsets(4)); 0344 otherwise 0345 error('Not implemented'); 0346 end 0347 end 0348 0349 % if available put the result from the last pass into data block 0350 if sc > 1 0351 t_sz = size(last_t); 0352 switch n_dims 0353 case 2 0354 t(1:t_sz(1), 1:t_sz(2)) = last_t; 0355 case 3 0356 t(1:t_sz(1), 1:t_sz(2), 1:t_sz(3)) = last_t; 0357 case 4 0358 t(1:t_sz(1), 1:t_sz(2), 1:t_sz(3), 1:t_sz(4)) = last_t; 0359 otherwise 0360 error('Not implemented'); 0361 end 0362 end 0363 0364 p_truncs = truncations(sc, p_dims); 0365 0366 for d = 1:n_p_dims 0367 0368 % move next dimension to X (columns) 0369 s = shifts(d); 0370 if s, t = permute(t, circshift(dims, s)); end 0371 0372 % hard work is in C 0373 t = do_iwtx(t, rh{d}, rg{d}, dlp(d), dhp(d), ... 0374 sc_levels(sc), ... 0375 p_truncs(d)); 0376 0377 end % for dim 0378 0379 % and back to original shape 0380 if end_shift 0381 t = permute(t, circshift(dims, end_shift)); 0382 end 0383 0384 % set up for next pass 0385 offsets = offsets + truncations(sc, :); 0386 last_t = t; 0387 0388 end % for scale 0389 0390 y = t; 0391 0392 %------------------------------ 0393 % END OF THE ALGORITHM 0394 %------------------------------ 0395