function SNR=computeSNR(clean,noisy)
% function SNR=computeSNR(clean,noisy)
% This function computes the SNR of a noisy signal given its (noiseless) reference
%
% INPUT: 
% noisy..........signal for which the SNR is computed
% clean..........reference signal (same length as clean)
%
% OUTPUT:
% SNR............signal-to-noise-ratio of the noisy signal
%
% Bernhard Geiger, SPSC, 2012-06-01

% make column vectors
clean=clean(:);
noisy=noisy(:);

% find possible time shift between signal and reference by correlation
[coeffs,lags]=xcorr(clean,noisy);
[~,maxlag]=max(coeffs);
maxlag=lags(maxlag);

% compute SNR (after correcting possible shift)
if maxlag<0
    SNR=10*log10(var(clean(1:end+maxlag))/var(clean(1:end+maxlag)-noisy(1-maxlag:end)));
else
    SNR=10*log10(var(clean(1+maxlag:end))/var(clean(1+maxlag:end)-noisy(1:end-maxlag)));
end