test_mw.m

多小波分析源码

end
if (~iszero(Mu(:,2) - mu1))
    error('failed: 1st continuous moment of DGHM - numerical');
end
if (~iszero(Mu(:,3) - mu2))
    error('failed: 2nd continuous moment of DGHM - numerical');
end

H = wavelet('dghm','symbolic');
Mu = continuous_moment(H,2);
if (~iszero(Mu(:,1) - mu0))
    error('failed: 0th continuous moment of DGHM - symbolic');
end
if (~iszero(Mu(:,2) - mu1))
    error('failed: 1st continuous moment of DGHM - symbolic');
end
if (~iszero(Mu(:,3) - mu2))
    error('failed: 2nd continuous moment of DGHM - symbolic');
end

disp('OK - continuous_moment');

% nullspace

v = rand(3,1);
V = v * v';
N = nullspace(V);
if (~iszero(V*N) | (size(N,2) ~= 2))
    error('failed: nullspace');
end

syms a b c real;
v = [a;b;c];
V = v * v';
N = nullspace(V);
if (~iszero(V*N) | (size(N,2) ~= 2))
    error('failed: nullspace');
end

disp('OK - nullspace');

% Refinement matrix

syms h0 h1 h2 h3;
H = mpoly({h0,h1,h2,h3},0,'symbol',2,1) / sqrt(2);
T1 = refinement_matrix(H);
T2 = [h0,0,0,0;
      h2,h1,h0,0;
      0,h3,h2,h1;
      0,0,0,h3] * sqrt(2);
if (~iszero(T1-T2))
    error('failed: refinement_matrix');
end

disp('OK - refinement_matrix');

% Correlation

H1 = wavelet('daub',2,'symbolic');
H2 = wavelet('daub',1,'symbolic');
H = correlation(H1,H2,'-+');
z = mpoly(1,1);
true = (H1 + z*H1)/2;
true = true(1,:);
if (~iszero(H - true))
    error('failed: correlation');
end

disp('OK - correlation');

% Transition matrix
% 
% This routine is just a composition of correlation and refinement matrix
% We test the approximation order: if \phi has approximation order,
% the transition matrix is supposed to have eigenvalues 1, 1/2, 1/4, 1/8

H = wavelet('bat');
T = transition_matrix(H);
s = eig(T);
for lambda = [1,1/2,1/4,1/8]
    index = find(abs(s - lambda) < 1.e-12);
    if (length(index) == 0)
	error('failed: transition_matrix');
    end
end

H = wavelet('bat','symbolic');
T = transition_matrix(H);
s = eig(T);
for lambda = [1,1/2,1/4,1/8]
    index = find(abs(double(s) - lambda) < 1.e-12);
    if (length(index) == 0)
	error('failed: transition_matrix');
    end
end

disp('OK - transition_matrix');

% Point values

H = wavelet('cl2','symbolic');
[xphi,phi,xpsi,psi] = point_values(H,1/2);
true = sym('[0,1/2,1,1/2,0;0,-sqrt(7)/4,0,sqrt(7)/4,0]');
if (~iszero(phi - true))
    error('failed: point_values');
end

disp('OK - point_values');

% Sobolev exponent

H = wavelet('cl2');
s = sobolev_exponent(H);
if (abs(s - 1.0545) > 1.e-4)
    error('failed: sobolev_exponent');
end
s = sobolev_estimate(H);
if (abs(s - 0.6406) > 1.e-4)
    error('failed: sobolev_estimate');
end

disp('OK - sobolev_exponent');

% range

v = [1;2;3;4];
vhat = v/norm(v);
V = v * v';
U = range(V);
if (~iszero(U*U' - vhat*vhat'))
    error('failed: range');
end

v = sym(v);
vhat = v/norm(v);
V = v * v';
U = range(V);
if (~iszero(U*U' - vhat*vhat'))
    error('failed: range');
end

disp('OK - range');

% projection factor

u = [1;2;3;4];
F = projection_factor(u,2,2);
if (~isidentity(F * F'))
    error('failed - projection_factor');
end

v = [3;1;4;1];
[F,Ft] = projection_factor(u,v,2,2);
if (~isidentity(F * Ft'))
    error('failed - projection_factor');
end

u = sym(u);
v = sym(v);
F = projection_factor(u,2,2);
if (~isidentity(F * F'))
    error('failed - projection_factor');
end
[F,Ft] = projection_factor(u,v,2,2);
if (~isidentity(F * Ft'))
    error('failed - projection_factor');
end

disp('OK - projection_factor');

% projection shift

[H,Ht] = wavelet('hm',0);
P = polyphase(H);
Pt = polyphase(Ht);
[Pnew,Ptnew,F] = projection_shift(P,Pt);
if (~iszero(P - Pnew * F{1}))
    error('failed: projection_shift');
end
if (~iszero(Pt - Ptnew * F{1}))
    error('failed: projection_shift');
end
if (Pt.min ~= Ptnew.min+1)
    error('failed: projection_shift');
end

[H,Ht] = wavelet('hm',0,'symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
[Pnew,Ptnew,F] = projection_shift(P,Pt);
if (~iszero(P - Pnew * F{1}))
    error('failed: projection_shift');
end
if (~iszero(Pt - Ptnew * F{1}))
    error('failed: projection_shift');
end
if (Pt.min ~= Ptnew.min+1)
    error('failed: projection_shift');
end

disp('OK - projection_shift');

% projection factorization

H = wavelet('dghm');
P = polyphase(H);
factors = projection_factorization(P);
product = factors{1};
for i = 2:length(factors)
    product = product * factors{i};
end
if (~iszero(product - P))
    error('failed - projection_factorization orthogonal numeric');
end

H = wavelet('dghm','symbolic');
P = polyphase(H);
factors = projection_factorization(P);
product = factors{1};
for i = 2:length(factors)
    product = product * factors{i};
end
if (~iszero(product - P))
    error('failed - projection_factorization orthogonal symbolic');
end

[H,Ht] = wavelet('hm',0);
P  = polyphase(H);
Pt = polyphase(Ht);
[factors,factorst] = projection_factorization(P,Pt);
product = factors{1};
for i = 2:length(factors)
    product = product * factors{i};
end
productt = factorst{1};
for i = 2:length(factorst)
    productt = productt * factorst{i};
end
if (~iszero(product - P))
    error('failed - projection_factorization biorthogonal numeric');
end
if (~iszero(productt - Pt))
    error('failed - projection_factorization biorthogonal numeric');
end

[H,Ht] = wavelet('hm',0,'symbolic');
P  = polyphase(H);
Pt = polyphase(Ht);
[factors,factorst] = projection_factorization(P,Pt);
product = factors{1};
for i = 2:length(factors)
    product = product * factors{i};
end
productt = factorst{1};
for i = 2:length(factorst)
    productt = productt * factorst{i};
end
if (~iszero(product - P))
    error('failed - projection_factorization biorthogonal symbolic');
end
if (~iszero(productt - Pt))
    error('failed - projection_factorization biorthogonal symbolic');
end

disp('OK - projection_factorization');


% TST

H = wavelet('dghm');
Ht = H;
[Hn,Htn] = tst(H,Ht,'lr');
if (approximation_order(Hn) ~= 3 | approximation_order(Htn) ~= 1)
    error('failed - tst');
end

[Hn,Htn] = itst(H,Ht,'lr');
if (approximation_order(Hn) ~= 1 | approximation_order(Htn) ~= 3)
    error('failed - tst');
end

H = wavelet('dghm','symbolic');
Ht = H;
[Hn,Htn] = tst(H,Ht,'lr');
if (approximation_order(Hn) ~= 3 | approximation_order(Htn) ~= 1)
    error('failed - tst');
end

[Hn,Htn] = itst(H,Ht,'lr');
if (approximation_order(Hn) ~= 1 | approximation_order(Htn) ~= 3)
    error('failed - tst');
end

disp('OK - tst/itst');

% prefilter and postfilter

s = rand(64,1);
sp = prefilter(s,'dghm');
s2 = postfilter(sp,'dghm');
if (~iszero(s - s2) | ~iszero(norm(s) - norm(sp)))
    error('failed - prefilter / postfilter');
end

% the roundoff error is bigger here than usual, since this
% prefilter is only given numerically.
sp = prefilter(s,'dghm','ap2');
s2 = postfilter(sp,'dghm','ap2');
if (~iszero(s - s2) | ~iszero(norm(s) - norm(sp)))
    error('failed - prefilter / postfilter');
end

disp('OK - prefilter/postfilter')

% DWT - IDWT

H = wavelet('cl3');
s = rand(256,1);
sd = dwt(s,H,3);
s2 = idwt(sd,H,3);
if (~iszero(s - s2))
    error('failed - DWT / IDWT');
end

[H,Ht] = wavelet('hm',0);
sd = dwt(s,H,3);
s2 = idwt(sd,Ht,3);
if (~iszero(s - s2))
    error('failed - DWT / IDWT');
end

disp('OK - DWT/IDWT');