test_mw.m

多小波分析源码

% TEST_MW -- test the subroutines in the mw directory

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

disp('testing routines in mw')
disp(' ');

% test wavelet

disp('testing creation of standard wavelets');
disp('this involves a lot of symbolic computation, and will take a while');

% daub
[H,Ht] = wavelet('daub',1);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('daub',1,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''daub'',1)');
end

[H,Ht] = wavelet('daub',2);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('daub',2,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''daub'',2)');
end

[H,Ht] = wavelet('daub',3);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('daub',3,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''daub'',3)');
end

[H,Ht] = wavelet('daub',10);
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''daub'',10)');
end
disp('OK - wavelet(''daub'')');

% coif6
[H,Ht] = wavelet('coif6',[-4,1]);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-4,1],'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-4,1])');
end

[H,Ht] = wavelet('coif6',[-4,1],2);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-4,1],2,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-4,1],2)');
end

[H,Ht] = wavelet('coif6',[-3,2]);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-3,2],'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-3,2])');
end

[H,Ht] = wavelet('coif6',[-3,2],2);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-3,2],2,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-3,2],2)');
end

[H,Ht] = wavelet('coif6',[-2,3]);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-2,3],'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-2,3])');
end

[H,Ht] = wavelet('coif6',[-2,3],2);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-2,3],2,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-2,3],2)');
end

[H,Ht] = wavelet('coif6',[-1,4]);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-1,4],'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-1,4])');
end

[H,Ht] = wavelet('coif6',[-1,4],2);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('coif6',[-1,4],2,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''coif6'',[-1,4],2)');
end

disp('OK - wavelet(''coif6'')');

% cohen
[H,Ht] = wavelet('cohen',1,1);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('cohen',1,1,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''cohen'',1,1)');
end

[H,Ht] = wavelet('cohen',2,4);
P = polyphase(H);
Pt = polyphase(Ht);
[Hs,Hts] = wavelet('cohen',2,4,'symbolic');
Ps = polyphase(Hs);
Pts = polyphase(Hts);
if (~isidentity(P*Pt') | ~isidentity(Ps*Pts'))
    error('failed: wavelet(''cohen'',2,4)');
end

disp('OK - wavelet(''cohen'')');

% daub79
[H,Ht] = wavelet('daub79');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''daub79'')');
end

disp('OK - wavelet(''daub79'')');

% bat
[H,Ht] = wavelet('bat','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''bat'')');
end

disp('OK - wavelet(''bat'')');

% cl2
[H,Ht] = wavelet('cl2','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''cl2'')');
end
[H,Ht] = wavelet('cl2',-3/5,'symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''cl2'',-3/5)');
end

disp('OK - wavelet(''cl2'')');

% cl3
[H,Ht] = wavelet('cl3','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''cl3'')');
end

disp('OK - wavelet(''cl3'')');

% daubbal
[H,Ht] = wavelet('daubbal',2,'symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''daubbal'',2)');
end

disp('OK - wavelet(''daubbal'')');

% dghm
[H,Ht] = wavelet('dghm','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''dghm'')');
end

disp('OK - wavelet(''dghm'')');

% stt
[H,Ht] = wavelet('stt','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''stt'')');
end

disp('OK - wavelet(''stt'')');

% bcl
[H,Ht] = wavelet('bcl','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''bcl'')');
end

disp('OK - wavelet(''bcl'')');

% hc
[H,Ht] = wavelet('hc','dahmen','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''hc'',''dahmen'')');
end

[H,Ht] = wavelet('hc','shortest','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''hc'',''shortest'')');
end

[H,Ht] = wavelet('hc','smoothest','symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''hc'',''smoothest'')');
end

disp('OK - wavelet(''hc'')');

% hm
[H,Ht] = wavelet('hm',0,'symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''hm'',0)');
end

[H,Ht] = wavelet('hm',1/4,'symbolic');
P = polyphase(H);
Pt = polyphase(Ht);
if (~isidentity(P*Pt'))
    error('failed: wavelet(''hm'',1/4)');
end

disp('OK - wavelet(''hm'')');
disp(' ')

% discrete moments

M0 = mpoly(sym('[3,2*sqrt(2);2*sqrt(2),1;2*sqrt(2),-4;0,0]')/5);
M1 = mpoly(sym('[3,0;6*sqrt(2),2;6*sqrt(2),-8;3,-3*sqrt(2)]')/10);
M2 = mpoly(sym('[3,0;9*sqrt(2),-1;9*sqrt(2),-11;9,-6*sqrt(2)]')/10);

H = wavelet('dghm');
M = discrete_moment(H,0:2);
if (~iszero(M(:,:,1) - M0))
    error('failed: 0th discrete moment of DGHM - numerical');
end
if (~iszero(M(:,:,2) - M1))
    error('failed: 1st discrete moment of DGHM - numerical');
end
if (~iszero(M(:,:,3) - M2))
    error('failed: 2nd discrete moment of DGHM - numerical');
end

H = wavelet('dghm','symbolic');
M = discrete_moment(H,0:2);
if (~iszero(M(:,:,1) - M0))
    error('failed: 0th discrete moment of DGHM - symbolic');
end
if (~iszero(M(:,:,2) - M1))
    error('failed: 1st discrete moment of DGHM - symbolic');
end
if (~iszero(M(:,:,3) - M2))
    error('failed: 2nd discrete moment of DGHM - symbolic');
end

disp('OK - discrete_moment');

% continuous moments

mu0 = mpoly(sym('[sqrt(6);sqrt(3);0;0]')/3);
mu1 = mpoly(sym('[sqrt(6);2*sqrt(3);0;0]')/6);
mu2 = mpoly(sym('[16*sqrt(6);52*sqrt(3);-3*sqrt(3);0]')/168);

H = wavelet('dghm');
Mu = continuous_moment(H,2);
if (~iszero(Mu(:,1) - mu0))
    error('failed: 0th continuous moment of DGHM - numerical');