thfig35.m

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

% thfig35: BP Thesis Figure 3.5	-- Dynamic Range, energy ratio = 1
%
%	  Matching Pursuit behaves badly when the signal is made up
%	  atoms within dynamic range
%	  Basis Pursuit has handles dynamic range problem pretty well
% Signal: Dynamic-0
%	x = EnergyRatio * (cos1 + cos2+ cos3 + cos4) +  dirac1 + dira2
%	four cosines have close frequencies.
%	energy ratio: EnergyRatio = 1
% Signal Length: 256
% Dictionary:
%	[cosines diracs], includes all the components of the signal.
% Observation:
%	(all the plots are transformed by x -> sign(x) * log(abs(x) + 1)
%	(a) the signal
%	(b) (c) Mathing Pursuit produces worse result as the energy ratio
%		increase from 1 to 10^6
%	  	thfig35: EnergyRatio = 1, it can pick up the two diracs
%	 	thfig36: EnergyRatio = 1e2, it can pick up the two diracs,
%			 but with some noise on the other diracs
%	  	thfig37: EnergyRatio = 1e4, it compleltely misses the two
%			 dirac.
%	(d) (e) BasisPursuit finds the exact decomposition, not affected by 
%	        the dynamic range problem
%

help thfig35
n = 256;
EnergyRatio = 1;

fineness = 4;
t = ((1:n)' - .5) / n;
const = (2/n) ^ .5;
x1 = const * cos(pi * (.5 + (126  - 1) / fineness) * t);
x2 = const * cos(pi * (.5 + (127  - 1) / fineness) * t);
x3 = const * cos(pi * (.5 + (128  - 1) / fineness) * t);
x4 = const * cos(pi * (.5 + (129  - 1) / fineness) * t);
zerosn = zeros(n,1);
x5 = zerosn; x5(230) = 1;
x6 = zerosn; x6(250) = 1;
x = EnergyRatio * (x1 + x2 + x3 + x4) + x5 + x6;

dicts = MakeList('DCT', 'DIRAC');
par1s = MakeList(fineness, 0);
par2s = MakeList(0, 0);
par3s = MakeList(0, 0);
m1 = SizeOfDict(n, 'DCT', fineness, 0, 0);
m2 = SizeOfDict(n, 'DIRAC', 0, 0, 0);
m = m1 + m2;

%MP
natom = n;
frac = 1e-5;
time = cputime;
cMP = MP(x, dicts, par1s, par2s, par3s, natom, frac);
cputime3 = cputime - time;

%BP
time = cputime;
cBP = BP_Interior(x, dicts, par1s, par2s, par3s, 1e-3, 1e-3, 1e-3, 1e-5, 1e-5);
cputime4 = cputime - time;

fprintf('\n')
fprintf('CPU Running Time of  MP = %8.4e\n', cputime3);
fprintf('CPU Running Time of  BP = %8.4e\n', cputime4);

%Plot
figure(1);clf
subplot(3,1,1);plot((x > 0) .* log(abs(x)+1) - (x < 0) .* log(abs(x)+1))
title('(a) Signal on log scale: 4 close cosines & 2 dirac; Energy Ratio = 1e0')
MyAxis(n);

temp = log(abs(cMP) + 1);
subplot(3,2,3);plot(temp(1:m1));
MyAxis(m1);
title('(b) MP: DCT Coefs on log scale');
axis323 = axis;

subplot(3,2,4);plot(temp((m1+1):m));
MyAxis(m2);
title('(c) MP: DIRAC Coefs on log scale');
axis324 = axis;
X1 = [230 230]; Y1 = [axis324(3) axis324(4)];
X2 = [250 250]; Y2 = [axis324(3) axis324(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off


temp = log(abs(cBP) + 1);
subplot(3,2,5);plot(temp(1:m1));
title('(d) BP: DCT Coefs on log scale');
axis(axis323)

subplot(3,2,6);plot(temp((m1+1):m));
title('(e) BP: DIRAC Coefs on log scale');
axis(axis324)
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off