ttfig3.m

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

% ttfig3: BasisPursuit Figure 33 -- Dynamic Range, energy ratio = 1
% Objective:
%	Decomposition Stability
%	When signals consists of features with dynamic range,
%	(1) MP usually misses features with low energy
%	(2) BP can separate features with different energy; thus
%	    does a good job in picking up features with low energy
% Signal:
%	Dynamic-1 = 1   * (cos1 + cos2+ cos3 + cos4) +  dirac1 + dira2
%	Dynamic-3 = 1e3 * (cos1 + cos2+ cos3 + cos4) +  dirac1 + dira2
%	Dynamic-4 = 1e4 * (cos1 + cos2+ cos3 + cos4) +  dirac1 + dira2
%
%	four cosines have close frequencies, not in the dictionary
% Sample Size:	256
% Dictionary:
%	4-overcomplete DCT  +  DIRAC
% Observation:
%	all the plots are transformed by x -> sign(x) * log(abs(x) + 1
%
%	Figure(1), figure(2), figure(3)
%	(a) Signal: Dynamic-1, Dynamic-3, Dynamic-4
%	(b) DCT coef from MP. MP cannot resolve the similarity among
%	    the four close cosines. Instead of giving a representation
%	    within the right frequency range, MP results in a symmetric
%	    decaying structure on all frequencies.
%	(c) DIRAC coef from MP.
%	    Since MP couldn't represent the four close cosines well,
%	    it produce residual messing up the diracs, a residual
%	    which gets bigger and bigger as the dynamic range increases.
%	    Consequently MP fails to reveal the two true diracs in the
%	    signal as dynamic range reaches 10^4.
%		figure(1): EnergyRatio = 1. MP can pick up the two diracs
%	 	figure(2): EnergyRatio = 100, MP can pick up the two diracs,
%			   but with some noise.
%		figure(3): EnergyRatio = 10^4, the two true diracs are
%			   buried in noise.
%	(d) DCT coef from BP. 
%	(e) DIRAC coef from BP.
%	    BP is able to separate the two sets of features with different
%	    energy. Unlike MP, BP finds the right frequency range and
%	    gives a concise representation of the four close cosines
%	    without any residual. Therefore BP is able to pick up the
%	    two DIRACs, no matter what the EnergyRatio is.
%
fig_label=get(0,'Children');
for fi=1:prod(size(fig_label))
   if fi>2
      close (fig_label(fi));
   end
end

help ttfig3

%Sample Size
n = 256;

%Dictionary: 4-overcomplete Discrete Cosines + Diracs
fineness = 4;
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;

%Signal: Dynamic-1, Dynamic-2, Dynamic-4
t = ((1:n)' - .5) / n;freq = pi * ((1:(4*n))' - 1) / (fineness*n);
freq1 = pi * (126.55-1) / (fineness*n);
freq2 = pi * (127.55-1) / (fineness*n);
freq3 = pi * (128.55-1) / (fineness*n);
freq4 = pi * (129.55-1) / (fineness*n);
const = (2/n) ^ .5;
x1 = const * cos(pi * ((126.55  - 1) / fineness) * t);
x2 = const * cos(pi * ((127.55  - 1) / fineness) * t);
x3 = const * cos(pi * ((128.55  - 1) / fineness) * t);
x4 = const * cos(pi * ((129.55  - 1) / fineness) * t);
zerosn = zeros(n,1);
x5 = zerosn; x5(230) = 1;x6 = zerosn; x6(250) = 1;

EnergyRatio = 1;
dynamic_1 = EnergyRatio * (x1 + x2 + x3 + x4) + x5 + x6;
EnergyRatio = 10^3;
dynamic_3 = EnergyRatio * (x1 + x2 + x3 + x4) + x5 + x6;
EnergyRatio = 10^4;
dynamic_4 = EnergyRatio * (x1 + x2 + x3 + x4) + x5 + x6;

%MP
natom = n;
frac = 1e-2;
cMP_1 = MP(dynamic_1, dicts, par1s, par2s, par3s, natom, frac);
natom = n;
frac = 1e-4;
cMP_3 = MP(dynamic_3, dicts, par1s, par2s, par3s, natom, frac);
natom = n;
frac = 1e-5;
cMP_4 = MP(dynamic_4, dicts, par1s, par2s, par3s, natom, frac);


%BP
cBP_1 = BP_Interior(dynamic_1, dicts, par1s, par2s, par3s, 1e-2, 1e-2, 1e-2);
cBP_3 = BP_Interior(dynamic_3, dicts, par1s, par2s, par3s, 1e-4, 1e-4, 1e-4);
cBP_4 = BP_Interior(dynamic_4, dicts, par1s, par2s, par3s, 1e-5, 1e-5, 1e-5);


%Plot for dynamic-1
x = dynamic_1; cMP = cMP_1; cBP = cBP_1;
figure(1);clf
subplot(3,1,1);plot((x > 0) .* log(abs(x)+1) - (x < 0) .* log(abs(x)+1))
title('(a) Dynamic-1 on log scale');
MyAxis(n);

cMP = log(abs(cMP) + 1);
cMPDCT = cMP(1:m1);
cMPDIRAC = cMP((m1+1):m);
subplot(3,2,3);plot(freq, cMPDCT);
axis323 = axis;
axis([freq(100) freq(160) axis323(3) axis323(4)]);
title('(b) MP: DCT Coefs on log scale');
axis323 = axis;
X1 = [freq1 freq1]; Y1 = [axis323(3) axis323(4)];
X2 = [freq2 freq2]; Y2 = [axis323(3) axis323(4)];
X3 = [freq3 freq3]; Y3 = [axis323(3) axis323(4)];
X4 = [freq4 freq4]; Y4 = [axis323(3) axis323(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
plot(X3, Y3, ':');
plot(X4, Y4, ':');
hold off

subplot(3,2,4);plot(cMPDIRAC);
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

cBP = log(abs(cBP) + 1);
cBPDCT = cBP(1:m1);
cBPDIRAC = cBP((m1+1):m);
subplot(3,2,5);plot(freq, cBPDCT);
axis([freq(100) freq(160) axis323(3) axis323(4)]);
title('(d) BP: DCT Coefs on log scale');
axis325 = axis;
X1 = [freq1 freq1]; Y1 = [axis325(3) axis325(4)];
X2 = [freq2 freq2]; Y2 = [axis325(3) axis325(4)];
X3 = [freq3 freq3]; Y3 = [axis325(3) axis325(4)];
X4 = [freq4 freq4]; Y4 = [axis325(3) axis325(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
plot(X3, Y3, ':');
plot(X4, Y4, ':');
hold off

subplot(3,2,6);plot(cBPDIRAC);
MyAxis(m2);
title('(e) BP: DIRAC Coefs on log scale');
axis326 = axis;
X1 = [230 230]; Y1 = [axis326(3) axis326(4)];
X2 = [250 250]; Y2 = [axis326(3) axis326(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

%Plot for dynamic-2
figure(3);clf
x = dynamic_3; cMP = cMP_3; cBP = cBP_3;
subplot(3,1,1);plot((x > 0) .* log(abs(x)+1) - (x < 0) .* log(abs(x)+1))
title('(a) Dynamic-2 on log scale');
MyAxis(n);

cMP = log(abs(cMP) + 1);
cMPDCT = cMP(1:m1);
cMPDIRAC = cMP((m1+1):m);
subplot(3,2,3);plot(freq, cMPDCT);
axis323 = axis;
axis([freq(100) freq(160) axis323(3) axis323(4)]);
title('(b) MP: DCT Coefs on log scale');
axis323 = axis;
X1 = [freq1 freq1]; Y1 = [axis323(3) axis323(4)];
X2 = [freq2 freq2]; Y2 = [axis323(3) axis323(4)];
X3 = [freq3 freq3]; Y3 = [axis323(3) axis323(4)];
X4 = [freq4 freq4]; Y4 = [axis323(3) axis323(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
plot(X3, Y3, ':');
plot(X4, Y4, ':');
hold off

subplot(3,2,4);plot(cMPDIRAC);
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

cBP = log(abs(cBP) + 1);
cBPDCT = cBP(1:m1);
cBPDIRAC = cBP((m1+1):m);
subplot(3,2,5);plot(freq, cBPDCT);
axis([freq(100) freq(160) axis323(3) axis323(4)]);
title('(d) BP: DCT Coefs on log scale');
axis325 = axis;
X1 = [freq1 freq1]; Y1 = [axis325(3) axis325(4)];
X2 = [freq2 freq2]; Y2 = [axis325(3) axis325(4)];
X3 = [freq3 freq3]; Y3 = [axis325(3) axis325(4)];
X4 = [freq4 freq4]; Y4 = [axis325(3) axis325(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
plot(X3, Y3, ':');
plot(X4, Y4, ':');
hold off

subplot(3,2,6);plot(cBPDIRAC);
MyAxis(m2);
title('(e) BP: DIRAC Coefs on log scale');
axis326 = axis;
X1 = [230 230]; Y1 = [axis326(3) axis326(4)];
X2 = [250 250]; Y2 = [axis326(3) axis326(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

%Plot for dynamic-1
figure(4);clf
x = dynamic_4; cMP = cMP_4; cBP = cBP_4;
subplot(3,1,1);plot((x > 0) .* log(abs(x)+1) - (x < 0) .* log(abs(x)+1))
title('(a) Dynamic-4 on log scale');
MyAxis(n);

cMP = log(abs(cMP) + 1);
cMPDCT = cMP(1:m1);
cMPDIRAC = cMP((m1+1):m);
subplot(3,2,3);plot(freq, cMPDCT);
axis323 = axis;
axis([freq(100) freq(160) axis323(3) axis323(4)]);
title('(b) MP: DCT Coefs on log scale');
axis323 = axis;
X1 = [freq1 freq1]; Y1 = [axis323(3) axis323(4)];
X2 = [freq2 freq2]; Y2 = [axis323(3) axis323(4)];
X3 = [freq3 freq3]; Y3 = [axis323(3) axis323(4)];
X4 = [freq4 freq4]; Y4 = [axis323(3) axis323(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
plot(X3, Y3, ':');
plot(X4, Y4, ':');
hold off

subplot(3,2,4);plot(cMPDIRAC);
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

cBP = log(abs(cBP) + 1);
cBPDCT = cBP(1:m1);
cBPDIRAC = cBP((m1+1):m);
subplot(3,2,5);plot(freq, cBPDCT);
axis([freq(100) freq(160) axis323(3) axis323(4)]);
title('(d) BP: DCT Coefs on log scale');
axis325 = axis;
X1 = [freq1 freq1]; Y1 = [axis325(3) axis325(4)];
X2 = [freq2 freq2]; Y2 = [axis325(3) axis325(4)];
X3 = [freq3 freq3]; Y3 = [axis325(3) axis325(4)];
X4 = [freq4 freq4]; Y4 = [axis325(3) axis325(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
plot(X3, Y3, ':');
plot(X4, Y4, ':');
hold off

subplot(3,2,6);plot(cBPDIRAC);
MyAxis(m2);
title('(e) BP: DIRAC Coefs on log scale');
axis326 = axis;
X1 = [230 230]; Y1 = [axis326(3) axis326(4)];
X2 = [250 250]; Y2 = [axis326(3) axis326(4)];
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

help ttfig3