example_mw.m

几个关于多小波的程序

% Example calls to subroutines in mw

% 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

format compact;

disp(' ');
disp('## This Matlab script file illustrates the use of the multiwavelet toolbox.');
disp('## A description of each step, and the command used to produced it,');
disp('## are shown AFTER each step. Scroll back to see the output.');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
H = wavelet('dghm','symbolic')
disp(' ');
disp('## Create the DGHM wavelet in symbolic form');
disp('## H = wavelet(''dghm'',''symbolic'')');
disp(' ');
dummy = input('Hit <return> to continue');

disp(' ');
P = polyphase(H)
disp(' ');
disp('## Convert it to polyphase form');
disp('## P = polyphase(H)');
disp(' ');
dummy = input('Hit <return> to continue');

disp(' ');
H = symbol(P)
disp(' ');
disp('## Convert it back to a symbol');
disp('## H = symbol(P)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
[p,Y] = approximation_order(H)
disp(' ');
disp('## Find its approximation order and approximation vectors');
disp('## [p,Y] = approximation_order(H)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
F = projection_factorization(H)
disp(' ');
disp('## Factor the polyphase matrix into projection factors');
disp('## F = projection_factorization(H)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
F{1}
disp(' ');
disp('## F{1} is a constant orthogonal matrix');
disp('## F{1}');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
F{2}
disp(' ');
disp('## F{2} is a projection factor');
disp('## F{2}');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
[H2,Ht2] = tst(H,H,'lr')
disp(' ');
disp('## Move one approximation order to the dual side, by TST');
disp('## [H2,Ht2] = tst(H,H,''lr'')');
disp(' ');
dummy = input('hit <return> to continue');


disp(' ');
p = approximation_order(H2)
disp(' ');
pt = approximation_order(Ht2)
disp(' ');
disp('## Check: p = approximation_order(H2), pt = approximation_order(Ht2)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
[H3,Ht3] = itst(H2,Ht2,'lr')
disp(' ');
disp('## Move the approximation order back, by ITST (this gives a new multiwavelet)');
disp('## [H3,Ht3] = itst(H2,Ht2,''lr'')');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
p = approximation_order(H3)
disp(' ');
pt = approximation_order(Ht3)
disp(' ');
disp('## Check: p = approximation_order(H3), pt = approximation_order(Ht3)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
H = double(H)
disp(' ');
disp('## The remaining computations need to be done numerically;');
disp('## They would be prohibitively slow or impossible in symbolic form');
disp(' ');
disp('## Convert H to numerical form');
disp('## H = double(H)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
s = sobolev_estimate(H)
disp(' ');
disp('## Quick and dirty estimate of Sobolev exponent');
disp('## s = sobolev_estimate(H)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
s = sobolev_exponent(H)
disp(' ');
disp('## Better estimate of Sobolev exponent');
disp('## s = sobolev_exponent(H)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
[xphi,phi,xpsi,psi] = point_values(H,1/64)
disp(' ');
disp('## Calculate point values, and plot them')
disp('## (If you have your plot window iconized, bring it up now and leave it up)');
disp('## [xphi,phi,xpsi,psi] = point_values(H,1/64)');
disp(' ');
figure(1)
clf
subplot(2,2,1)
plot(xphi,phi(1,:));
title('\phi_1');
subplot(2,2,2)
plot(xphi,phi(2,:));
title('\phi_2');
subplot(2,2,3)
plot(xpsi,psi(1,:));
title('\psi_1');
subplot(2,2,4)
plot(xpsi,psi(2,:));
title('\psi_2');
dummy = input('hit <return> to continue');

s = [(1:128)/128,cos((1:128)*pi/128)]';
disp(' ');
disp('## The next several steps illustrate the discrete wavelet transform');
disp(' ');
disp('## Plot the original signal (straight line, followed by cosine)');
clf
plot(s);
axis([0,256,-2,2]);
title('original signal');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
sp = prefilter(s,'dghm','ap2')
disp(' ');
disp('## Do prefiltering, and plot the results');
disp('## sp = prefilter(s,''dghm'',''ap2'')');
plot(1:128,reshape(sp,2,128));
axis([0,128,-2,2]);
title('components of coefficient vector after prefiltering');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ');
sd = dwt(sp,H,2)
disp(' ');
disp('## Perform DWT through 2 levels');
disp('## sd = dwt(sp,H,2)');
disp(' ');
disp('Sort out the pieces, and plot them');
clf
subplot(4,1,1)
plot(1:128,reshape(sp,2,128))
axis([0,128,-2,2]);
title('preprocessed original signal s_2')
subplot(4,1,2)
plot(1:32,reshape(sd(1:64),2,32))
axis([0,32,-5,5]);
title('coarse approximation s_0')
subplot(4,1,3)
plot(1:32,reshape(sd(65:128),2,32))
axis([0,32,-.003,.003]);
title('fine detail d_0')
subplot(4,1,4)
plot(1:64,reshape(sd(129:256),2,64))
axis([0,64,-.005,.005]);
title('fine detail d_1')
disp(' ');
dummy = input('hit <return> to continue');

disp(' ')
sp2 = idwt(sd,H,2)
disp(' ');
norm(sp-sp2)
disp(' ');
disp('## Do reconstruction')
disp('## sp2 = idwt(sd,H,2)');
disp(' ');
disp('## Check accuracy: norm(sp - sp2)');
disp(' ');
dummy = input('hit <return> to continue');

disp(' ')
s2 = postfilter(sp2,'dghm','ap2');
disp(' ');
norm(s-s2)
disp(' ');
disp('## Do Postfiltering')
disp('## s2 = postfilter(sp2,''dghm'',''ap2'')');
disp(' ');
disp('## Check accuracy: norm(s - s2)');