%
%  Satellite Remote Sensing
%  May 10, 2005 - David Sandwell
%
%  DEMONSTRATION OF THE SPEED OF THE CONVOLUTION THEOREM
%  (need file band1c.raw)
%
%  This code performs a high-pass filter of band-1 of the
%  San Diego Landsat image using two methods.
%
%  Method A. convolution by fast fourier transform (fft)
%  Method B. convolution by direct summation
%
%  The size of the convolution filter is increased
%  to show the ratio of direct summation speed
%  speed to fft speed.
%
close all; clear;
%
%  load band-1 of the Landsat data
%
bl = double(fread( fopen( 'band1c.raw', 'r'), [1500 1500], '*uint8'))';
figure(1), colormap('gray'), imagesc(bl), colorbar, title('band 1');
pause
%
%  construct the x and y meshgrids
%
[ny,nx]=size(bl);
x = -nx/2:nx/2-1;
y = -ny/2:ny/2-1;
[X,Y] = meshgrid(x,y);
R2=X.*X+Y.*Y;
%
%  make a Gaussian high-pass filter with a half-width of 10 pixels
%
sig = 10;
fltg = exp(-R2/(2*sig*sig));
fltg = fltg/(sum(sum(fltg)));
flt = zeros(ny,nx);
flt(ny/2+1,nx/2+1) = 1.;
flt = flt-fltg;
aflt=abs(fftshift(fft2(flt)));
figure(2), colormap('gray'), imagesc(aflt), colorbar, title('spectrum of band-pass filter')
pause
%
%  do the convolution using fft and establish elapse time
%
tic
cflt = fft2(flt);
cbl=fft2(bl);
outfft=real(fftshift((ifft2(cbl.*cflt))));
time_fft=toc
figure(3), colormap('gray'), imagesc(outfft,[-64,64]), colorbar, title('convolution by fft');
pause
%
%  do the filter by direct convolution
%
for nc = 4:4:128;
nc2=nc/2;
ncy0=ny/2+1-nc2;
ncyf=ny/2+1+nc2;
ncx0=nx/2+1-nc2;
ncxf=nx/2+1+nc2;
fltc=flt(ncy0:ncyf,ncx0:ncxf);
fltc=fltc-(mean(mean(fltc)));
tic
outconv=conv2(bl,fltc,'same');
time_conv2=toc;
[nc, time_conv2/time_fft]
figure(3), colormap('gray'), imagesc(outconv,[-64,64]), colorbar, title('convolution by direct summation');
%
%  look at the difference between the two methods
%
diffout = outconv - outfft;
figure(4), colormap('gray'), imagesc(diffout,[-32,32]), colorbar, title('difference');
%
pause
end
