MATLAB Examples

Test E: Reconstruction from Gradients

This script is supplemental material to the paper:

"An Algebraic Framework for the Real-Time Solution of Inverse Problems on Embedded Systems"

by: Matthew Harker, Christoph Gugg, Paul O'Leary and Gerhard Rath

This script require our Discrete orthogonal polynomial toolbox, it is available at:

http://www.mathworks.com/matlabcentral/fileexchange/41250

and the function to convert figures to eps:

http://www.mathworks.com/matlabcentral/fileexchange/42388

Contents

Prepare to execute code.

close all
clear all;
%
% Set some defaults
%
fontSize = 10;
set(0,'DefaultFigureColor',[1,1,1]);
set(0,'DefaultaxesFontName','Times');
set(0,'DefaultaxesFontSize',fontSize);
set(0,'DefaulttextFontName','Times');
set(0,'DefaulttextFontSize',fontSize);
set(0,'DefaultfigurePaperType','A4');
set(0,'DefaultTextInterpreter', 'latex');
%
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.13 0.1 0.86 0.86];
set(0,'DefaultaxesPosition',MyAxesPosition);
%

Load the Test Data

load BVPFitData;
%

Pre-calculations for the Reconstruction

generate the particular and homogeneous solutions.

ls = 11;
%
name = ['BVP',int2str(ls)];
%
% generate a set of constrained basis functions.
%
noBfs = length(x) - 1;
noBfs = 11;
[yc, Bc, S] = dopGenConstrained( x, noBfs, T, ls, ls);
%
D = dopDiffLocal( x, ls, ls);
%
% define the linear differential aperator
%
L = D;
%
% Use the basis functions to generate a regularizing operator
%
M = Bc * pinv( L * Bc );
%
% Compute the particular solution
%
yp = yc - M * L * yc;
%
% Compute the covariancde propagation
%
C = (M * M') ;
%
% The square root of the diagonal of the covariance matrix is the
% unscaled standard deviation.
%
s = sqrt( diag( C ) );
%
Using a local approx for M

Perofrm a Monte Carlo Simulation for the Reconstruction

noSims = 10000;
%
% Define a noise gain in terms of the maximum derivative, in this
% example 1% of the maximum derivative. This corresponds to 1% noise
% in the measurement signal.
%
noiseGain = 0.01;
normDy = max( dy );
stdIn = noiseGain * normDy;
%
% Prepare stoorage for the solutions
%
sols = zeros( length(x), noSims );
G = zeros( length(x), noSims );
%
% Run the Monte Sarlo Simulations
%
for k=1:noSims
    %
    % generate the perturbed forcing function
    %
    gn = stdIn * randn( length(x), 1);
    g = dy + gn;
    G(:,k) = g;
    %
    % Compute the solution
    %
    yrk = M * g + yp;
    sols(:,k) = yrk;
end;

Present the Results of the Monte Carlo Simulation

This fugure shows the analytical derivatives and the perturbations at each node.

FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.15 0.15 0.81 0.8];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
figGs = figure;
errorbar( x, dy, stdIn * ones(size(x)), '-kx');
xlabel('$$x$$');
ylabel('Perturbed $$D y$$');
grid on;
axis([0,1,-1,1]);

Compute the mean reconstruction and plot all the reconstruction points.

FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.15 0.15 0.77 0.8];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
% Compute the mean solution
%
sm = mean(sols')';
%
% Scale the errors to make them visible.
%
scale = 10;
%
figRs = figure;
plot( x, y, 'k-');
hold on;
errorbar( x, sm, scale * stdIn*s, 'kx');
grid on;
xlabel('$$x$$');
ylabel('Reconstructions');
axis([0,1,-0.15,0.25]);
hold on;
[n,m] = size(T);
legend( 'Analitical', 'Reconstruction', 'Location', 'Southwest');
for k=1:n
    if T(k,1) == 0
        % is a value constraint
        plot( T(k,2), T(k,3), 'ko', 'MarkerFaceColor', 'k');
    end;
end;
%

Plot the Mean error of the reconstructions. this verifies if there is a systematic reconstruction error.

FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.15 0.15 0.84 0.8];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
scale = 1E5;
%
figME = figure;
plot( x, scale*(y - sm), '-k.');
xlabel('$$x$$');
ylabel('$$\bar{\epsilon} \, \, \, (\times 10^5)$$');
grid on;

Plot the standard deviation of the reconstructions

Scale the standard deviation to make it visible

stdS = std( sols' )';
scale = 1E3;
%
figVer = figure;
plot( x, scale * stdS, '-k.');
grid on;
%
xlabel('$$x$$');
ylabel('$$\mathrm{std}(\epsilon) \, \, \, (\times 10^3)$$');
hold on;
plot( x, scale*stdIn*s, '-r.');
[n,m] = size(T);
for k=1:n
    if T(k,1) == 0
        % is a value constraint
        plot( T(k,2), 0, 'ko', 'MarkerFaceColor', 'k');
    end;
end;
%

Present the covariance matrix as an image.

figCov = figure;
imagesc( C );
axis image;
colorbar;
xlabel('$$x$$');
ylabel('$$x$$');
%
figMStd = figure;
subplot(2,1,1);
plot( x, scale*(y - sm), '-k.');
ylabel('$$\bar{\epsilon} \, \, \, (\times 10^5)$$');
grid on;
%
hold on;
[n,m] = size(T);
for k=1:n
    if T(k,1) == 0
        % is a value constraint
        plot( T(k,2), 0, 'ko', 'MarkerFaceColor', 'k');
    end;
end;
%
subplot(2,1,2);
plot( x, scale * stdS, '-k.');
grid on;
xlabel('$$x$$');
ylabel('$$\mathrm{std}(\epsilon) \, \, \, (\times 10^3)$$');
hold on;
plot( x, scale*stdIn*s, 'ko');
[n,m] = size(T);
for k=1:n
    if T(k,1) == 0
        % is a value constraint
        plot( T(k,2), 0, 'ko', 'MarkerFaceColor', 'k');
    end;
end;
axis([0,1,0,2]);

legend( 'P','MC','location','northEast');
%

Save the figures as .eps files for possible publication.

h = @figure2eps;
hInfo = functions( h );
if isempty( hInfo.file )
    disp( 'You will need to download the function figure2eps if you');
    disp( 'wish to generate the eps files for the figures' );
    disp( 'The function is available at:');
    disp( 'http://www.mathworks.at/matlabcentral/fileexchange/42388');
else
    %
    figure2eps( figMStd, [name, 'BiasAndStd']);
    figure2eps( figGs, 'BVPGs');
    figure2eps( figVer, [name, 'confid']);
    figure2eps( figRs, [name, 'Resonstrions']);
    figure2eps( figME, [name, 'meanE']);
    figure2eps( figCov, [name, 'CovL']);
end;