% Eric Thrane
% An example matched filtering calculation
%
clear all;

% trials
T=100;

% sampling frequency
fs = 16384;

% define time
t = 1/fs:1/fs:1;
t = repmat(t, T, 1);

% define signal to be a sine-Gaussian
s = 0.11*exp(-25*(t-0.5).^2) .* sin(2*pi*105*t);
sf = fft(s, [], 2);

% create noise
n = randn(size(t));
nf = fft(n, [], 2);

% create average noise power spectra
Sh = mean(abs(nf).^2,1);
Sh = repmat(Sh, T, 1);

% add signal to noise
h = s + n;
hf = fft(h, [], 2);

% calculate matched filtering SNR
% see Eq. 2.3 in Owen and Sathyaprakash paper
S = real(inner_product(hf, sf, Sh));
N = real(inner_product(nf, sf, Sh));

% display matched filterint SNR denoted rho
fprintf('rho = %1.3e +/- %1.3e\n', mean(S)/std(N), std(S/std(N)));

% compare loudest signal event with loudest noise event
fprintf('max signal and noise in %i trials is: %1.1e, %1.1e\n', ...
  T, max(S), max(N));

% make a plot
figure;
plot(t, n(1,:), 'r');
hold on;
plot(t, s(1,:), 'b');
xlabel('t');
ylabel('h');
title(['\rho = ' num2str(mean(S)/std(N))], 'FontSize', 20);
pretty;
print('-dpng', 'time_series');
