clear all
% Create a template of a triangle pointed up
N = 19; 
c = (N+1)/2;
t = zeros(c,N);
for idx = 0:c-1
  [c-idx c+idx]
  t(idx+1,c-idx:c+idx) = 1;
end
[nrt nct] = size(t);

M = 256;
% Insert template into image
f = zeros(M);
f((M/2-nrt/2+1):(M/2+nrt/2),(M/2-c+1):(M/2+c-1)) = t;

normf = sum(sum(f.^2));
SNR = [-10 0];cnt = 1;
for idx = 1:length(SNR)
  sigma = sqrt(N/normf * 10^(-SNR(idx)/10));
  % Ignore the mcidx loop
  for mcidx = 1:1
    % Make noisy data
    y(:,:,cnt) = f+sigma*randn(size(f));
    % Do the matched filtering
    b(:,:,cnt) = conv2(y(:,:,cnt),flipud(fliplr(t)));
    cnt = cnt+1;
  end
end

% Plot the results
figure(1);clf;colormap(gray)
imagesc(f);axis image;title('Noise image')
print -djpeg ps02p3c1.jpg

figure(1);clf;colormap(hsv)
subplot(221);imagesc(y(:,:,1));axis image
title('Signal at -10dB SNR')
subplot(222);mesh(b(:,:,1));
title('Correlation output at -10dB SNR')
subplot(223);imagesc(y(:,:,2));axis image
title('Signal at 0dB SNR')
subplot(224);mesh(b(:,:,2));
title('Correlation output at 0 dB SNR')