thfig55.m

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

% thfig55: BP Thesis Figure 5.5 -- TV DeNoise and Dictionary Merger
%	
%	BP with Heaviside dictionaries is equivalent to TV DeNoising
% Signal:	HeaviSine
% Sample Size:	1024
% Dictionary:	
%	(1) Heavisides
%	(2) Waves: Orthgonal Wavelets with qmf=(Symmlet, 8)
%	(3) Waves+Jumps: Orthgonal Wavelets with qmf=(Symmlet, 8) normalized 
%		by TV + Jumps with decay rate .9
% Observations:
%	(a) Signal HeaviSine
%	(b) Noisy HeaviSine
%	(c) Plots of sorted coefficients (compression rate in the three
%	    different dictionaries):
%	(d) BP Denoising reconstruction using Heavisides (TV DeNoising)
%	(e) BP Denoising reconstruction using Waves
%	(f) BP Denoising reconstruction using Waves+Jumps
%
%	Notice that the quality of reconstruction depends on 
%	how well the dictionary compresses the signal. Since
%	Waves+Jumps give the best compression for HeaviSine, it
%	delivers the best reconstruction.
%

help thfig55

n = 1024; D = log2(n); par1 = 3;
t = (1:n)'/n;
qmf1 = MakeONFilter('Symmlet', 8);
qmf2 = MakeONFilter('Haar');
x = InputSignal('HeaviSine', n);

	%%Noiseless case
	%cHS = BP_Interior(x,'HS',0,0,0,1e-2,1e-2,1e-3);
	cHS = [x(1);diff(x)];
	cSymm = IWT_PO(x,par1,qmf1);

	dicts = MakeList('POTV', 'JUMP');
	par1s = MakeList(par1, .9);
	par2s = MakeList(qmf1, 0);
	par3s = MakeList(0, 0);
	cMerg = BP_Interior(x, dicts, par1s, par2s, par3s);


	%%DeNoising
	randn('seed',931316785);
	SNR = 7;

	%%HeaviSide Dictionary
	[x y] = NoiseMaker(x,SNR);
	[xrecHS crecHS] = BPDN_Interior(y, 'HS', 0, 0, 0, 1e-1, 1e-1, 1e-3);

	%%Symmlet
	cclean = SoftThresh(FWT_PO(y,par1, qmf1), (2 * log(n)) ^ .5); 
	xrecSymm = IWT_PO(cclean, par1, qmf1);

	%%Merge
	dicts = MakeList('POTV', 'JUMP');
	par1s = MakeList(par1, .9);
	par2s = MakeList(qmf1, 0);
	par3s = MakeList(0, 0);
	[xrecMerg crecMerg] = BPDN_Interior(y, dicts, par1s, par2s, par3s,1e-1,1e-1,1e-2);

%%Plots	
	
	figure(1);clf
	subplot(3,2,1); plot(t,x); title('(a) Signal: HeaviSine');
	titlestr = sprintf('(b) SNR = %g', SNR');
	subplot(3,2,2); plot(t,y); title(titlestr);

	subplot(3,2,4); plot(t,xrecHS); title('(d) BPDeNoise: Heaviside');
	subplot(3,2,5); plot(t,xrecSymm); title('(e) Wavelet Shrinkage: Symmlet');
	subplot(3,2,6); plot(t, xrecMerg); title('(f) BPDeNoise: Jump+Wave');
	subplot(3,2,1); AXIS = axis;
	subplot(3,2,1); axis(AXIS);
	subplot(3,2,2); axis(AXIS);
	subplot(3,2,4); axis(AXIS);
	subplot(3,2,5); axis(AXIS);
	subplot(3,2,6); axis(AXIS);

	subplot(3,2,3); 
	semilogy(reverse(sort(abs(cMerg)))/max(abs(cMerg)));
	axis([0 n 1e-6 1])
	title('(c) Sorted Coefs')
	hold on;
	xlabel('Order')
	ylabel('Amplitude')
	semilogy(reverse(sort(abs(cHS)))/max(abs(cHS)), ':');
	semilogy(reverse(sort(abs(cSymm)))/max(abs(cSymm)), '--');
	MyAxis(200);
	axes('Position',[.24 .44 .9 .9],'Visible','off');
	h1 = text(0, 0, 'Dotted: Heaviside');
	set(h1, 'FontSize', 6);
	h2 = text(0, .03, 'Dashed: Wave'); 
	set(h2, 'FontSize', 6);
	h3 = text(0, .06, 'Solid: Jump+Wave');
	set(h3, 'FontSize', 6);
	hold off