Calculating the SNR of a bandpass delta-sigma modulator in Matlab

% This calculates the SNR of a bandpass delta-sigma modulator
% Alper Ucar, ucar {at} ieee {dot} org
% -------------------------------------------------------
% [snrv,PxtdB,PxnbdB] = snr(data,tone,fftlen,osr,offset)
% data: output stream in time-domain
% tone: input stream in time-domain
% fftlen: FFT length
% osr: oversampling ratio
% offset: number of samples to be removed from both sides
% of the spectral content for the measurement
% PxtdB: input power
% PxnbdB: noise power
% snrv: SNR
function[snrv,PxtdB,PxnbdB] = snr(data,tone,fftlen,osr,offset)

clc
close all
N = 2^floor(log2(fftlen));
data = squeeze(data(1:N));
tone = squeeze(tone(1:N));
% hann-windowing
W = hann(N);
% coherent gain factor
CGF = 8/3;
x = data.*W;
xt = tone.*W;
% Normalized spectral content from 0 to Fs
x = abs(fft(x))/N;
x = x(1:N/2);
xt = abs(fft(xt))/N;
xt = xt(1:N/2);
% Search the max signal amplitude
[val pos] = max(xt(2:end));
% remove the input from the spectral content
% with coherent sampling and Hann windowing the tone consists
% of only 3 samples
xn = x;
xn((pos-offset):(pos+offset))= min(x);
% Extract in-band noise and calculate its power
xnb = xn(floor(pos-(N/(4*osr)):(floor(pos+(N/(4*osr))))));
Pxnb = 2*CGF*xnb.^2;
Pxnb = sum(Pxnb);
PxnbdB = 10*log10(Pxnb);
% input power
Pxt = 2*CGF*xt.^2;
Pxt = sum(Pxt);
PxtdB = 10*log10(Pxt);
% disp(['input signal power = ',num2str(PxtdB),'']);
% calculate SNR
snrv = PxtdB-PxnbdB;
% Plot PSD
Px = 2*CGF*x.^2;
PxdB = 10*log10(Px);
figure;
subplot 211
plot(PxdB)
axis([0 0.5*fftlen min(PxdB) 0])
xlabel('Frequency bin')
ylabel('PSD [dBFS]')
title('Output PSD')
text(50,-20,sprintf('SNR: %4.2fdB\nOSR: %d',snrv,osr),...
'BackgroundColor',[1 1 1]);
grid on
subplot 212
Pxn = 2*xn.^2;
PxndB = 10*log10(Pxn);
plot(PxndB)
axis([0 0.5*fftlen min(PxndB) 0])
xlabel('Frequency bin')
ylabel('PSD [dBFS]')
title('Output PSD')
text(50,-20,sprintf('Noise Power: %4.2f\nInput Power: %4.2f',...
PxnbdB,PxtdB),'BackgroundColor',[1 1 1]);
grid on