% thfig25: BP Thesis Figure 2.5 -- Super-resolving two close frequencies
%		within the Rayleigh distance
%
%       BP is able to achieve super-resolution, while
%	MOF and MP sometimes cannot
%
% Signal:	TwinSine-1 (two frequencies in the dictionary)
%		cos(pi * w1 * t) + cos(pi * w2 * t)
% Signal Length: 256
% Dictionary:	factor 4-overcomplete discrete cosine dictionary
% Observation:
%	(a) The signal: TwinSine-1
%	(b) MOF results a Dirichlet-like oscillatory structure, consisting
%	    not of two but of many frequencies.
%	(c) MP fails to superresolve. MP initially selects the frequency
%	    in between the two frequencies making up the signal. Because of
%	    this mistake, MP is forced to make a series of alternating
%	    corrrections which suggest a highly complex and organized
%	    strucuture.
%	(d) BP exactly picks up the doublet structure.
%

help thfig25

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

%MOF
time = cputime;
cMOF = MOF(x, 'DCT', fineness, 0, 0);
cputime1 = cputime - time;

%MP
natom = n;
frac = 1e-2;
time = cputime;
cMP = MP(x, 'DCT', fineness, 0, 0, n, frac);
cputime3 = cputime - time;

%BP
time = cputime;
cBP = BP_Interior(x, 'DCT', fineness, 0, 0, 1e-2, 1e-2, 1e-2);
cputime4 = cputime - time;

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

%True frequencies
P11 = [(126-1)/(fineness * n) (126-1)/(fineness * n)];
P12 = [-.5 2];
P21 = [(128-1)/(fineness * n) (128-1)/(fineness * n)];
P22 = P12;


%Plots
figure(1);clf
subplot(2,2,1);plot(t, x);title('(a) Signal: TwinSine-1')
subplot(2,2,2);
PlotCoeff(cMOF, n, 'DCT', fineness, 0, 0);
hold on;
axis([(101-1)/(fineness*n) (150-1)/(fineness*n) -.5 2])
plot(P11, P12, ':')
plot(P21, P22, ':')
hold off
xlabel('Frequency/Nyquist')
ylabel('Amplitude')
title('(b) MOF coefs')

subplot(2,2,3); 
PlotCoeff(cMP, n, 'DCT', fineness, 0, 0);
title('(c) MP Coefs');
hold on;
axis([(101-1)/(fineness*n) (150-1)/(fineness*n) -.5 2])
plot(P11, P12, ':')
plot(P21, P22, ':')
hold off
xlabel('Frequency/Nyquist')
ylabel('Amplitude')

subplot(2,2,4);
PlotCoeff(cBP, n, 'DCT', fineness, 0, 0);
title('(d) BP Coefs');
hold on;
axis([(101-1)/(fineness*n) (150-1)/(fineness*n) -.5 2])
plot(P11, P12, ':')
plot(P21, P22, ':')
hold off
xlabel('Frequency/Nyquist')
ylabel('Amplitude')