function [snrdB,ptotdB,psigdB,pnoisedB] = calcSNR_lu(vout,Fin,Fs,N,R)

% vout: input size is  data(1,N)

% Fin : signal 's frequency in analog domain

% Fs : sampling rate for current stage

% N:    samples number for current stage

% R:    total oversampling  factor for current stage after previous stage's

% decimation

% snrdB:     SNR in dB

% ptotdB:    Bit-stream power spectral density (vector)

% psigdB:    Extracted signal power spectral density (vector)

% pnoisedB:  Noise power spectral density (vector)

%

N=floor(N);

%w=hann(N)';         % windowing vector

w=blackman(N,'periodic')'; 

f=Fin/Fs; % Normalized signal frequency，Fin=Fs*nper/N;=> f=nper/N % Normalized signal frequency，Fin=Fs*nper/N;=> f=nper/N

fB=N/2/R; % Base-band frequency bins,always keep same，bw=Fs/2/R;=>fB=N/2/R,R is osr in current stage

fB=ceil(fB);

signal=(N/sum(w))*sinusx(vout(1:N).*w,f,N); % Extracts sinusoidal signal

noise=vout(1:N)-signal;             % Extracts noise components

stot=((abs(fft((vout(1:N).*w)'))).^2); % Bit-stream PSD

ssignal=(abs(fft((signal(1:N).*w)'))).^2; % Signal PSD

snoise=(abs(fft((noise(1:N).*w)'))).^2; % Noise PSD

pwsignal=sum(ssignal(1:fB));                % Signal power

pwnoise=sum(snoise(5:fB));             % Noise power,skip DC

snr=pwsignal/pwnoise;

snrdB=dbp(snr);

norm=sum(stot)/sum(vout(1:N).^2)*N; % PSD normalization

if nargout > 1

ptot=stot/norm;

ptotdB=dbp(ptot);

end

if nargout > 2

psig=ssignal/norm;

psigdB=dbp(psig);

end

if nargout > 3

pnoise=snoise/norm;

pnoisedB=dbp(pnoise);

end

function SNR = snr_qian_lvkun(s_y,frq,sfrq,bw)

% s_y: data

% frq: Input frequency

% sfrq: sampling frequency

% bw: bandwitdh of the signal

% frq/fres= frq/sfrq*Ns=int number,=>Ns=int_number*sfrq/frq

format long

Ns=20*sfrq/frq;

s_y=s_y(1:Ns);

%Ns <=length(s_y)

w = hanning(Ns, 'periodic');

z = 1;

[a, fre] = pwelch(s_y, w, [], Ns, sfrq);

fres = sfrq / Ns;

bx = round(Ns * bw/sfrq + 1);

bbpwr = fres * sum(a(z+2:bx));

sigbin = 1+ round(frq/fres);

sigbinmin = max([(sigbin-z) 1]);

sigbinmax = sigbin+z;

sigpwr = fres * sum(a(sigbinmin:sigbinmax));

SNR = 20*log10(sqrt(sigpwr/(bbpwr-sigpwr)));

figure

semilogx(fre, pow2db(a), 'r');

xlabel('Frequency(Hz)', 'FontWeight', 'Bold', 'FontSize',10);

ylabel('PSD(dB/Hz)', 'FontWeight', 'Bold', 'FontSize',10);

title('CIFB3L7B output PSD');

txt = ['\leftarrow SNR=', num2str(SNR)];

text(20e3, pow2db(max(a))-10, txt);

grid on;

hold on;

plot([fres 20e3], [pow2db(max(a)) pow2db(max(a))], 'b--', 'LineWidth', 1);

plot([20e3 20e3], [pow2db(min(a)) pow2db(max(a))], 'b--', 'LineWidth', 1);

plot([20e3 0.5*sfrq], [pow2db(min(a)) pow2db(min(a))], 'b--', 'LineWidth', 1);