thfig53.m

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

% thfig53: BP Thesis Figure 5.3 -- BP DeNoising Super Resolution
%
%       Basis Pursuit is able to achieve super-resolution, while 
%	Matching Pursuit sometimes cannot
% Signal: TwinSine-2, SNR = 10
%	x = cos(pi * w1 * t) + cos(pi * w2 * t)
% Signal Length: 256
% Dictionary:
%	4-fold overcomplete discrete cosine dictionary
% Observation:
%	(a) The signal: TwinSine-2
%	(b) The inner product b/w the signal and the basis functions. clearly
%	    MP is going to select the basis function with the highest
%	    inner product, which is neither w1, nor w2, but between. Then
%	    the residule messes up the whole thing
%	(c) MP couldn't resolve w1 and w2
%	(d) The representation from Basis Pursuit has two picks, correponding
%	    to w1 and w2

help thfig53

n = 256;
fineness = 4;
m = n * fineness;
t = ((1:n)' - .5) / n;
const = (2/n) ^ .5;
x1 = const * cos(pi * (125.2  - 1) / fineness * t);
x2 = const * cos(pi * (127.2  - 1) / fineness * t);
%x1 = const * cos(pi * (126  - 1) / fineness * t);
%x2 = const * cos(pi * (128  - 1) / fineness * t);

x = x1 + x2;
SNR = 10;
%noise level = 1
randn('seed', 123456789);
[x y] = NoiseMaker(x, SNR);

%MOFDeNoise
time = cputime;
[xrecMOF cMOF]= MOFDN(y, 'DCT', fineness, 0, 0);
cputime1 = cputime - time;

%MPDeNoise
time = cputime;
[xrecMP cMP] = MPDN(y, 'DCT', fineness, 0, 0);
cputime3 = cputime - time;

%BPDeNoise
time = cputime;
[xrecBP cBP] = BPDN_Interior(y, 'DCT', fineness, 0, 0, 1e-2, 1e-2, 1e-2);
cputime4 = cputime - time;

fprintf('\n')
fprintf('CPU Running Time of MOFDN = %8.4e\n', cputime1);
fprintf('CPU Running Time of  MPDN = %8.4e\n', cputime3);
fprintf('CPU Running Time of  BPDN = %8.4e\n', cputime4);

%Plots
freqs = (0:(m-1))'/m;
figure(1);clf;
subplot(3,2,1);plot(t,x);title('(a) TwinSine-2')
titlestr = sprintf('(b) Noised TwinSine, SNR = %g', SNR);
subplot(3,2,3);plot(t, y);title(titlestr);
subplot(3,2,5);plot(freqs, FastDCTAnalysis(x, fineness));
title('(c) DCT transform')
xlabel('Frequency/Nyquist')

subplot(3,2,2); plot(freqs, cMOF);
title('(d) MOF Coefs')
xlabel('Frequency/Nyquist')

subplot(3,2,4); plot(freqs, cMP);
title('(e) MP Coefs')
xlabel('Frequency/Nyquist')

subplot(3,2,6); plot(freqs, cBP);
title('(f) BP Coefs')
xlabel('Frequency/Nyquist')

subplot(3,2,5);
AXIS = axis;
axis([(101-1)/(fineness*n) (150-1)/(fineness*n) AXIS(3) AXIS(4)]);
AXIS = axis;
X1 = [(125.2-1)/(fineness * n) (125.2-1)/(fineness * n)];
Y1 = [AXIS(3) AXIS(4)];
X2 = [(127.2-1)/(fineness * n) (127.2-1)/(fineness * n)];
Y2 = Y1;
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

subplot(3,2,2);
axis(AXIS);
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

subplot(3,2,4);
axis(AXIS);
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

subplot(3,2,6)
axis(AXIS);
hold on
plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off