% Program to generate signal with noise and background with noise, fit them
% via a least squares fitting function provided by matlab, and plot the
% progress of the fit with number of scans.
% The function is a constant background with amplitude = 25, a peak signal
% with height = 2, centred at x = 250 and with a SD value of 16.

clear;
close all;
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
vid = VideoWriter('snr.mp4','MPEG-4');
vid.FrameRate = 60;    % Default 30
vid.Quality = 100;    % Default 75
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);

philBlue = [0.3 0.28 0.55];

numScans = 2000; % Number of scans
scanLength = 500; % Length of any given scan
sigPos = scanLength/2; % Place peak signal ideally at centre of scan
sigma = 16; % SD in Gaussian description of peak signal
xx = 1:scanLength; % x-coordinate of scan
bckav = 25; % Average background signal
sigpk = 2; % Average signal peak above background
backSig = zeros(1,scanLength); % Initialized background signal
sig = zeros(1,scanLength); % Initialized peak signal

% Initialize plots of progress in fitting
sigPkArray = zeros(1,numScans);
sigPosArray = zeros(1,numScans);
sigSigmaArray = zeros(1,numScans);
bckArray = zeros(1,numScans);

for i = 1:numScans
    subplot(4,5,[1,2,3,4,6,7,8,9,11,12,13,14,16,17,18,19])
    newplot
    hold off
    % Simulated background with Poisson noise
    backSig = backSig+poissrnd(bckav*ones(1,scanLength));
    % Idealized peak signal
    perfSig = sigpk*exp(-(xx - scanLength/2).^2/(2*sigma^2));
    % Simulated peak signal with Poisson noise
    sig = sig+poissrnd(perfSig);
    % Add up simulated noisy signals and normalized by dividing by i
    totSig = (sig + backSig)/i;
    plot(xx,totSig,'color',philBlue,'LineWidth', 2);
    set(gca,'linewidth',2,'fontsize',18);
    ymin = 0.98*min(totSig);
    ymax = 1.02*max(totSig);
    ylim([ymin ymax]);
    fun = @(x,xx)x(1)*exp(-(xx-x(2)).^2/(2*x(3)^2))+x(4);
    x0 = [10,scanLength/2.2,sigma*2,10]; % Random starting guess
    x = lsqcurvefit(fun,x0,xx,totSig);
    hold on
    fity = x(1)*exp(-(xx-x(2)).^2/(2*x(3)^2))+x(4);
    plot(xx,fity,'color','red','LineWidth', 2);
    set(gca,'linewidth',2,'fontsize',14);
    str1 = 'Scan ';
    str2 = num2str(i,'% 4.0f');
    strTot = [str1,str2];
    hTextYpos = ymin+(ymax-ymin)*0.88;
    hText = text(0.88*scanLength, hTextYpos, strTot, 'FontSize',18);
    sigPkArray(i) = x(1);
    sigPosArray(i) = x(2);
    sigSigmaArray(i) = x(3);
    bckArray(i) = x(4);
    
    
    %hold on
    subplot(4,5,[5])
    plot(sigPkArray(1:i),'color','red','LineWidth', 2);
    set(gca,'linewidth',2,'fontsize',14);
    xlim([1 numScans]);
    ylim([sigpk*0.5 sigpk*1.5]);
    hText = text(numScans*0.8, sigpk*1.25, 'a_0', 'FontSize',14);
    
    subplot(4,5,[10])
    plot(sigPosArray(1:i),'color','red','LineWidth', 2);
    set(gca,'linewidth',2,'fontsize',14);
    xlim([1 numScans]);
    ylim([sigPos*0.98 sigPos*1.02]);
    hText = text(numScans*0.8, sigPos*1.01, 'a_1', 'FontSize',14);
    
    
    subplot(4,5,[15])
    plot(sigSigmaArray(1:i),'color','red','LineWidth', 2);
    set(gca,'linewidth',2,'fontsize',14);
    xlim([1 numScans]);
    ylim([sigma*0.5 sigma*1.5]);
    hText = text(numScans*0.8, sigma*1.25, 'a_2', 'FontSize',14);
    
    subplot(4,5,[20])
    plot(bckArray(1:i),'color','red','LineWidth', 2);
    set(gca,'linewidth',2,'fontsize',14);
    xlim([1 numScans]);
    ylim([bckav*0.98 bckav*1.02]);
    hText = text(numScans*0.8, bckav*1.01, 'a_3', 'FontSize',14);
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mpg file
close(vid);
