Inverse Problems Lecture 2/2017: first encounter with convolution in Matlab

Показать описание

This is screen capture of Matlab programming I did when teaching my course Inverse Problems at University of Helsinki.

The lecture was given at January 20, 2017.

Here is the final code:

% Simple illustration of convolution in 1D
%
% Samuli Siltanen January 2017

% Parameters for controlling the plot appearance
fsize = 16;
lwidth = 2;

% Build a Point Spread Function (PSF)
M = 13;
psf = ones(1,2*M+1); % This makes sure psf has a unique centre
%psf = psf/sum(psf); % Normalization of the psf

% Construct "unknown signal" f
N = 400;
f = zeros(N,1);
f(1:(end/2)) = 1;

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

% Simulate the "measurement"
m = A*f;

% Simulate "noisy measurement"
sigma = .01;
mn = A*f + sigma*randn(N,1);

% Perform naive inversion
%f0 = inv(A)*m;
%fn = inv(A)*mn;
f0 = pinv(A)*m;
fn = pinv(A)*mn;
%f0 = A\m;

% Take a look
figure(1)
clf
subplot(3,1,1)
plot(f,'k','linewidth',lwidth)
hold on
plot(mn,'r','linewidth',lwidth)
set(gca,'ytick',[0 max(f)],'fontsize',fsize)
subplot(3,1,2)
plot(f,'k','linewidth',lwidth)
hold on
plot(f0,'b','linewidth',lwidth)
set(gca,'ytick',[0 max(f)],'fontsize',fsize)
subplot(3,1,3)
plot(f,'k','linewidth',lwidth)
hold on
plot(fn,'b','linewidth',lwidth)
set(gca,'ytick',[0 max(f)],'fontsize',fsize)