% plot recovered background (figure 3)
clearvars;
if ~exist('first','var')
    run('custom_settings')
    load([getenv('datadir') 'simulated'],'gendata');
    voigtp = load([getenv('datadir') 'simulated_voigt']);
    nmf = load([getenv('datadir') 'simulated_nmf']);
    first = 1;
end
opt.fig.size = 90*[opt.fig.col2width 2.5];

N = length(gendata.A);
W = length(gendata.B);
Mest = 1000;

V_Bhat = mean(voigtp.Bm(end-Mest:end,:));
N_Ahat = mean(nmf.Am(:,:,end-Mest:end),3);
% normalize
factors = max(N_Ahat);
N_Ahat=N_Ahat./repmat(factors,W,1);

figure
clf
plot(1:W,V_Bhat)
xlim([1 W])
hold on
plot(1:W,N_Ahat(:,2))
ylabel(['Intensity'])
xlabel(['Wavenumber'])
prepfig('',opt.fig)
plot(1:W,gendata.B,'--','linewidth',2)
opt.fig.setProps = false;
prepfig([strFIGdir 'simulated_background'],opt.fig)