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

clear

% create some functions to filter

step = zeros(1,100);
step(25:74) = 1;

t = 0:100;
sinusoid =     cos(2*pi*0.025*t) + ...
   0.5*cos(2*pi*0.05*t+pi/2) + ...
   0.25*cos(2*pi*0.1*t+2*pi/3);

subplot(4,2,1), plot(step)
subplot(4,2,2), plot(sinusoid)


STEP = fft(step);
SINUSOID = fft(sinusoid);

subplot(4,2,3), plot( fftshift(abs(STEP)) )
subplot(4,2,4), plot( fftshift(abs(SINUSOID)) )

% convolve them with an averaging filter
f1 = step;
f2 = sinusoid;

g = ones(1,15) / 15;

h1 = conv(g, f1);
h2 = conv(g, f2);

subplot(4,2,5), plot(h1)
subplot(4,2,6), plot(h2)

% look at Fourier transforms of the filtered signals
H1 = fft(h1);
H2 = fft(h2);

subplot(4,2,7), plot( fftshift(abs(H1)) )
subplot(4,2,8), plot(fftshift(abs(H2)))

pause


% repeat the process with some noise
% convolve them with an averaging filter
f1 = step + 0.25 * randn(size(step));
f2 = sinusoid + 0.25 * randn(size(sinusoid));

subplot(4,2,1), plot(f1);
subplot(4,2,2), plot(f2);

F1 = fft(f1);
F2 = fft(f2);

subplot(4,2,3), plot(fftshift(abs(H1)))
subplot(4,2,4), plot(fftshift(abs(H2)))

g = ones(1,15) / 15;

h1 = conv(g, f1);
h2 = conv(g, f2);

subplot(4,2,5), plot(h1)
subplot(4,2,6), plot(h2)

% look at Fourier transforms of the filtered signals
H1 = fft(h1);
H2 = fft(h2);

subplot(4,2,7), plot(fftshift(abs(H1)))
subplot(4,2,8), plot(fftshift(abs(H2)))

pause


% do the process in the Fourier domain
g1 = zeros(1,100);
g1(1:8) = 1;
g1(94:100) =1;
g1 = g1 / 15;
g2 = zeros(1,101);
g2(1:8) = 1;
g2(end-6:end) =1;
g2 = g2 / 15;

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

H1 = G1 .* F1;
H2 = G2 .* F2;

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

subplot(6,2,1), plot(f1)
subplot(6,2,2), plot(f2)
subplot(6,2,3), plot(fftshift(abs(F1)))
subplot(6,2,4), plot(fftshift(abs(F2)))
subplot(6,2,5), plot(g1)
subplot(6,2,6), plot(g2)
subplot(6,2,7), plot(fftshift(abs(G1)))
subplot(6,2,8), plot(fftshift(abs(G2)))
subplot(6,2,9), plot(fftshift(abs(H1)))
subplot(6,2,10), plot(fftshift(abs(H2)))
subplot(6,2,11), plot(real(h1))
subplot(6,2,12), plot(real(h2))

