% Exercise 2, M2
clear all; close all; clc;

% a) Choose n and evaluate the continuous convolution
% Choose the number of data points
n = 64;

% Construct data evaluation points
xvec = (1/n)*[0:(n-1)];

% Parameter determining the width of the PSF
a = .04;

% Create integration points for the interval [-a,a].
% Here we just use the midpoint rule.
Ny = 2000;
y  = linspace(-a,a,Ny);
Dy = y(2)-y(1);
psf1=PSF(y,a);

% Initialize the result
m_tilde = zeros(size(xvec));

% Loop over x values
for iii = 1:length(xvec)
    % Evaluate the integrand with the current x
    x = xvec(iii);
    intgnd = psf1.*targetf(x-y);
    
    % Compute the integral
    m_tilde(iii) = Dy*sum(intgnd);
end

m_tilde=reshape(m_tilde,length(m_tilde),1); % force the result to be column vector

% Plot
figure(1)
finex = linspace(0,1,2000);
plot(finex,targetf(finex),'k','linewidth',1) % plot target function
hold on
plot(xvec,m_tilde,'b.','markersize',3)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% b) Compute discrete convolution
% Define the PSF
% Choose nu to be the smallest integer such that
% nu >= na

M=ceil(n*a); % in fact, since psi(x) = 0 when x=a or x=-a, we can also
% choose M = ceil(n*a)-1;
x=(1/n)*[-M:1:M];
psf=PSF(x,a);


% Construct the convolution matrix
A = convmtx(psf,n);
A = (1/n)*A(:,(M+1):(end-M));

% Compute the discrete convolution
m = A*targetf(xvec).';

% Plot
plot(xvec,m,'r.','markersize',3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% c) Relative square norm error
err_rel=norm(m_tilde-m)/norm(m_tilde);

% find (roughly) n so that the relative error between m and m_tilde less
% than 1%. Result: n is about 1024
% Error < 0.1%, n is about 8182
n=64;
a = .04;

Ny = 2000;
y  = linspace(-a,a,Ny);
Dy = y(2)-y(1);
psf1=PSF(y,a);

xvec = (1/n)*[0:(n-1)];

% Initialize the result
m_tilde = zeros(size(xvec));

% Loop over x values
for iii = 1:length(xvec)
    % Evaluate the integrand with the current x
    x = xvec(iii);
    intgnd = psf1.*targetf(x-y);
    
    % Compute the integral
    m_tilde(iii) = Dy*sum(intgnd);
end

m_tilde=reshape(m_tilde,length(m_tilde),1);

M=ceil(n*a);
x=(1/n)*[-M:1:M];
psf=PSF(x,a);

% Construct the convolution matrix
A = convmtx(psf,n);
A = (1/n)*A(:,(M+1):(end-M));

% Compute the discrete convolution
m = A*targetf(xvec).';

err_rel=norm(m_tilde-m)/norm(m_tilde);

while err_rel>0.001
    n=n*2;
    xvec = (1/n)*[0:(n-1)];
    
    % Initialize the result
    m_tilde = zeros(size(xvec));

    % Loop over x values
    for iii = 1:length(xvec)
        % Evaluate the integrand with the current x
        x = xvec(iii);
        intgnd = psf1.*targetf(x-y);
        
        % Compute the integral
        m_tilde(iii) = Dy*sum(intgnd);
    end
    
    m_tilde=reshape(m_tilde,length(m_tilde),1);
    
    M=ceil(n*a);
    x=(1/n)*[-M:1:M];
    psf=PSF(x,a);
    
    % Construct the convolution matrix
    A = convmtx(psf,n);
    A = (1/n)*A(:,(M+1):(end-M));
    
    % Compute the discrete convolution
    m = A*targetf(xvec).';
    
    err_rel=norm(m_tilde-m)/norm(m_tilde);
end

% Plot
figure(2)
finex = linspace(0,1,2000);
plot(finex,targetf(finex),'k','linewidth',1) % plot target function
hold on
plot(xvec,m_tilde,'b.','markersize',3)
plot(xvec,m,'r.','markersize',3)
