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

## 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 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;