thfig38.m

％ Atomizer Main Directory, Version .802 里面信号含有分解去噪合成过程的代码 ％---------------------------------------

% thfig38: BP Thesis Figure 38 -- Phase plane evolution at
%	BP-Interior iterations.
%
%	BP-Interior iteratively sparsifies the MOF solution.
%
% Signal:	FM-Cosine
% Signal Length: 1024
% Dictionary: Cosine Packets with D = log2(n)-4
% Observations:
%	(a) Initial phase plane, i.e. the MOF phase plane.
%	(b) Phase plane at the 1st BP-Interior iteration.
%	(c) Phase plane at the 2st BP-Interior iteration.
%	(d) Phase plane at the 3st BP-Interior iteration.
%	(e) Phase plane at the 4st BP-Interior iteration.
%	(f) Phase plane at the BP-Interior iteration teramination.
%

help thfig38

n = 1024; D = log2(n); qmf = MakeONFilter('Symmlet', 8);bell = 'Sine';
x = InputSignal('FM', n);
cBPMovie= zeros(n * (D-4+1), 6);

NameOfDict = 'CP'; par1 = log2(n) - 4; par2 = 0; par3 = 0;
MaxBPIter = 30;
CGAccuracy = 1e-1;
FeaTol = 1e-1;
OptTol = 1e-1;

%*** Size of the problem
b = x(:); n = length(b); 
[m L] = SizeOfDict(n, NameOfDict, par1, par2, par3);

%*** Setting the controlling parameters
delta = OptTol;
gamma = OptTol;
sigma = .995;
MaxCGIter = m * log2(m);

%*** Initializing (x, y, z), make sure x and z are INTERIOR 
%*** Initializing by MOF representation
fprintf('\nBP-Interior:\n\n');
disp('Initializing BP-Interior by MOF ...')
x = MOF(b, NameOfDict, par1, par2, par3,1e-5);
x = [ x .* (x > 0) ; (-x) .* (x < 0) ];
sg = sign(x);
y = FastSS(sg, n, NameOfDict, par1, par2, par3);
ap = FastAA(y,NameOfDict,par1,par2, par3);ap = ap(1:m);
mp  = 1.1 * max(max(abs(ap)));
ap = ap ./ mp;
y = y ./ mp;
z1 = 1- ap(:);
z2 = 1+ ap(:);
z = [z1 ; z2];
x = x + .1;
x0 = x; y0 = y; z0 = z;

%*** Initialization for BP-Interior Iterations
mu0 = .01;
f = [z .* x; b - FastSS(x, n, NameOfDict, par1, par2, par3) - delta * y; ...
	- FastAA(y, NameOfDict, par1, par2, par3) - z + gamma * x + 1];
mu = mu0 * norm(f) / sqrt(2 * m);
k = -1;
gamma = gamma ^ 2;
delta = delta ^ 2;
zerosn = zeros(n,1);


%*** Iteration
disp(' ')
disp('BP-Interior Iterations ...')
head = '       mItn	PrimalFea   DualFea     DualGap         Obj           CGStep';
disp(head)

for i = 1:MaxBPIter,

	if i >= 1 & i <= 5,
		cBPMovie(:,i) = x(1:m)-x((m+1):(2*m));
	end

	% Calculate the Residules
	v1 = mu - z .* x;
	v2 = b - FastSS(x, n, NameOfDict, par1, par2, par3);
	v3 = - FastAA(y, NameOfDict, par1, par2, par3) - z + 1;
	
	% Test Termination
	PrimalFeas = norm(v2) / (1 + norm(x));
	DualFeas   = norm(v3) / (1 + norm(y));
	PDGap = z' * x / (1 + abs(sum(x)));
	Obj = sum(abs(x));
	fprintf('BPStep %4g:	%8.2e    %8.2e    %8.2e    %14.7e    %3g\n', i-1, PrimalFeas, DualFeas, PDGap, Obj, k)
	Feasible = PrimalFeas < FeaTol & DualFeas < FeaTol;
	Converged = PDGap < OptTol;
	if Feasible & Converged,
		break;
	end

	% Solve Newton's Direction by CG: (SHA + delta I) * dy = h
	H = 1 ./ (z ./ x + gamma);
	h = v2 - FastSS(H .* (v1 ./ x - v3), n, NameOfDict, par1, par2, par3);
	thresh = CGAccuracy/i * norm(h); thresh = thresh ^ 2;
	
	dy = zerosn;
	r = h;
	rho_1 = norm(r) ^ 2;
	for k = 1:MaxCGIter,
		%fprintf('	CG Step: %3g:	%g\n', k, rho_1);
		if k == 1,
			p = r;
		else
			beta = rho_1 / rho_0;
			p = r + beta * p;
		end
		w = FastSS(H .* FastAA(p, NameOfDict, par1, par2, par3), n, ...
			NameOfDict, par1, par2, par3) + delta * p;
		alpha = rho_1 / (p' * w);
		r = r - alpha * w;
		rho_0 = rho_1;
		rho_1 = norm(r) ^2;
		dy = dy + alpha * p;
		if rho_1  < thresh,
			break;
		end
	end

	if k == MaxCGIter,
		if imp < 0,
			x = xold;
		end
		fprintf('Too many CG iterations!\n')
		break;
	end


	dx=H .* (v1 ./ x - v3) + H .* FastAA(dy, NameOfDict, par1, par2, par3);
	dz = v1 ./ x - dx .* z ./ x;

	negdx = dx < -1e-14;
        negdz = dz < -1e-14;
        Rp = min(x(negdx) ./ abs(dx(negdx)));
        Rd = min(z(negdz) ./ abs(dz(negdz)));

	% Update (x y z mu)
	xold = x;
	x = x + sigma * Rp * dx; 
	y = y + sigma * Rd * dy; 
	z = z + sigma * Rd * dz;

	mu = (1 -  min([Rp Rd sigma])) * mu;
end


%% Feasibility:
%% Let dx = x - x0
%% Project dx onto the null space of \Phi by Conjugate Gradients:
%%	dx = (I - \Phi^T * (\Phi * \Phi^T)^-1 * \Phi) dx
%% x = x0 + dx
%%
dx = x - x0;
r = FastSS(dx, n, NameOfDict, par1, par2, par3);
thresh = 1e-16 * norm(r); thresh = thresh ^ 2;
dd = zerosn;
	rho_1 = norm(r) ^ 2;
	for k = 1:MaxCGIter,
		%fprintf('	CG Step: %3g:	%g\n', k, rho_1);
		if k == 1,
			p = r;
		else
			beta = rho_1 / rho_0;
			p = r + beta * p;
		end
		w = FastSS(FastAA(p, NameOfDict, par1, par2, par3), n, ...
			NameOfDict, par1, par2, par3);
		alpha = rho_1 / (p' * w);
		r = r - alpha * w;
		rho_0 = rho_1;
		rho_1 = norm(r) ^2;
		dd = dd + alpha * p;
		if rho_1  < thresh,
			break;
		end
	end
dx = dx - FastAA(dd, NameOfDict, par1, par2, par3);

%dx = dx - FastAA(FastSS(dx, n, NameOfDict, par1, par2, par3), ...
%	NameOfDict, par1, par2, par3) / (2 * L);
x = x0 + dx;
c = x(1:m) - x((m+1):(2*m));
cBPMovie(:,6) = c;
fprintf('Obj = %14.7e\n', sum(abs(c)));

TFScale = 20;
figure(1);clf
for i = 1:6,
	subplot(3,2,i);
	PhasePlane(cBPMovie(:, i), 'CP', n, 0, 256, TFScale);
end
subplot(3,2,1);title('Phase Plane: BP Iteration = 0');
subplot(3,2,2);title('Phase Plane: BP Iteration = 1');
subplot(3,2,3);title('Phase Plane: BP Iteration = 2');
subplot(3,2,4);title('Phase Plane: BP Iteration = 3');
subplot(3,2,5);title('Phase Plane: BP Iteration = 4');
subplot(3,2,6);title('Phase Plane: BP Termination');