N = 256;
M = 32;

f = zeros(N, N);
f(N/2 - M : N/2 + M, N/2 - M : N/2 + M) = 200;

figure(1);
colormap(gray(256));
image(f);

imwrite(uint8(f), 'f.jpg');
    
for sigma = 2:3:15;
    h = compute_psf(N, sigma);
    
    figure(2);
    colormap(gray(256));
    imagesc(h);
    
    g = filter2(h, f);

    figure(3);
    colormap(gray(256));
    image(g);  
    
    hfile = sprintf('h_%d.jpg', sigma);
    gfile = sprintf('g_%d.jpg', sigma);
    
    imwrite(uint8(h), hfile);
    imwrite(uint8(g), gfile);

end