function out = inner_product(a, b, Sh, deltaF)
% function out = inner_product(a, b, Sh)
% E Thrane: caculate matched filter inner product as defined in Eq. 2.4 of
% "Matched filtering of gravitational waves from inspiraling compact binaries:
% Computational cost and template placement" by Owen and Sathyaprakash
% a = (LxT) fft of strain data
% b = (LxT) fft of template
%   L = length of data
%   T = number of trials
% Sh = noise power spectrum: either (LxT) or (Lx1)
% deltaF (optional) = frequency resolution of FFT; deltaF should not be
%   specified if a and b are calculated using fft_eht
% -------------------------
% If the user supplies arrays of ffts and power spectra, then the dimensions
% should be size(sf) = NTrials x NFreqBins so that size(out) = NTrials x 1

% if deltaF is defined scale FFTs for appropriate comparison with the PSD
try
  a = sqrt(deltaF)*a;
  b = sqrt(deltaF)*b;
catch
  % do nothing
end

% check that array sizes match
if sum(size(Sh)~=size(a)) | sum(size(Sh)~=size(b))
  Sh = repmat(Sh, 1, size(a,2));
end

out = sum( (conj(a).*b)./Sh , 1);

return
