-
Notifications
You must be signed in to change notification settings - Fork 1
/
d2boxstep.m
101 lines (87 loc) · 2.91 KB
/
d2boxstep.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
function varargout=d2boxstep(x,level,dim,cofs)
% f=D2BOXSTEP(x,level,dim,cofs)
%
% Performs one FORWARD iteration (ending at a certain level of
% decomposition) of D2 transform along a certain dimension; putting
% scaling coefficients in front, followed by wavelet coefficients. This
% transform is an isometry. It is equal to the Haar transform.
%
% INPUT:
%
% x The three-dimensional array, sized as a power of two
% level The level we end up at [scalar]
% dim The index identifying the dimension [scalar]
% cofs The wavelet and scaling filter coefficients [D2BOXCOF]
%
% OUTPUT:
%
% f The wavelet transform of x, same dimensions as x
%
% EXAMPLE: My forward is the inverse of my inverse
% d2boxstep('demo1')
%
% SEE ALSO: D2BOXSTEPI, D2BOXCOF
%
% Inspired by Ignace Loris ([email protected]) on 22.06.2009
% Last modified by fjsimons-at-alum.mit.edu, 08/23/2010
if ~isstr(x)
% Initialize output, which you need to always take from the previous step
f=x;
if level==0
% Do nothing
varargout={f};
return
end
% Figure out dimensions
nall=size(x);
if length(nall)==2
nall(3)=1;
end
% To move into 'level' we split 2^(n-level+1) coefficients into two sets
% of k coefficients each, where 2^n is the dimension of the data set
k=size(x,dim)/2^level;
% The LF-tap filter length
LF=length(cofs.H0);
if k<=2^(LF/2); warning('Input signal is not long enough for reconstruction'); end
% Exclude the same number on the left and right, namely none
for i=1:k
% Isolate the DOWNSAMPLED sets of planes in the right dimension
xinside=[x(dindeks(2*i+0,dim,nall))'; ...
x(dindeks(2*i-1,dim,nall))'];
% And put the convolutions in the right spot
% See under ANALYSIS, SN p. 123
% Interior, Lowpass (Scaling Coefficients)
f(dindeks( i,dim,nall))=cofs.H0*xinside;
% Interior, Highpass (Wavelet Coefficients)
f(dindeks(k+i,dim,nall))=cofs.H1*xinside;
end
% All other coefficients (wavelet coeff from previous step) remain
varargout={f};
elseif strcmp(x,'demo1')
cofs=d2boxcof; dim=ceil(rand*3);
% Random-sized array must be at least 2^3 long
n=ceil(rand*7+3);
% Edge treatment only accurate to level 0 if it is 2^3 long
level=ceil(rand*(n-3));
% The random data
if dim==1
x=rand([2^n 1 1 ]);
elseif dim==2
x=rand([1 2^n 1 ]);
elseif dim==3
x=rand([1 1 2^n]);
end
% Initial output
disp(sprintf('\n====== D2BOXSTEP versus D2BOXSTEPI ===== \n'))
disp(sprintf('n = %i ; lev = %i ; dim = %i',n,level,dim))
% The forward transform
xf=d2boxstep(x,level,dim,cofs);
% The inverse of the forward transform
xfi=d2boxstepi(xf,level,dim,cofs);
% The reconstruction error, should be zero
mae1=mean(mean(mean(abs(x-xfi))));
% Further output
% Further output
disp(sprintf('no preconditioning mean(abs(error)) = %8.3e',mae1))
disp(sprintf('\n======================================= \n'))
end