% convolution/correlation demonstration
%
% create some 1D and 2D signals, convolution kernels, and demonstrate
%  the results of the convolutions
%

% set up some test signals

% get some images
im1 = zeros(32);
im1(8:23,8:23) = 1;

im2 = double(pnmread( 'washington2.pgm' ));



subplot(2,2,1), surf(im1);
subplot(2,2,2), image(im2+1), axis image off, colormap(gray(256))


% convolve them with an averaging filter

g = ones(7) / 49;

h1 = conv2(im1, g, 'same');
h2 = conv2(im2, g, 'same');

subplot(2,2,3), surf(h1)
subplot(2,2,4), image(h2+1), axis image off, colormap(gray(256))

pause



% do the process in the Fourier domain

im1 = im1 + 0.1 * randn(size(im1));
im2 = im2 + 20 * randn(size(im2));

IM1 = fft2( im1 );
IM2 = fft2( im2 );

g1 = zeros(size(im1));
g1(1:4, 1:4) = 1;
g1(30:32, 1:4) =1;
g1(1:4, 30:32) =1;
g1(30:32, 30:32) =1;
g1 = g1 / 49;

[nr,nc] = size(im2);
g2 = zeros(size(im2));
g2(1:4, 1:4) = 1;
g2(1:4, nc-2:nc) = 1;
g2(nr-2:nr, 1:4) = 1;
g2(nr-2:nr, nc-2:nc) = 1;
g2 = g2 / 49;

G1 = fft2(g1);
G2 = fft2(g2);

H1 = G1 .* IM1;
H2 = G2 .* IM2;

h1 = ifft2(H1);
h2 = ifft2(H2);

subplot(2,2,1), surf(im1)
subplot(2,2,2), image(im2+1), axis image off, colormap(gray(256))
subplot(2,2,3), surf(real(h1))
subplot(2,2,4), image(real(h2)+1), axis image off, colormap(gray(256))

