function [R,varx]=bt_psd(x,L,NFFT)

% function [R,r]=bt_psd(x,L)
% Computation of the PSD directly from the
% autocorrelation using Hanning windowing.
% The mean value is removed, i.e., x=x-mean(x).
% The PSD is normalized so that all power is
% located in the positive frequency axis,
% i. e., sum(R)=r(0)=sum(x.^2)/Lx.
% Note that the usual estimator for the
% variance is: var(x)=sum(x.^2)/(Lx-1).
% This normalization plots R(0) 50% down relative
% to the rest of the PSD graph.
% x=input signal
% L=window length
% NFFT=FFT length, NFFT >= 2*L, N is power of 2.
% R(fk) is computed for fk=k/N, k=0,1,...,N/2.
% B=4/(3*L); % Normalized Hanning window bandwidth
% nu=(8*length(x))/(3*L); % Degrees of freedom (Hanning)
% See: Jenkins-Watts, Chapters 5-7.
%
% Programmed by: Dimitris Manolakis, 1994
%
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------



N=2^nextpow2(2*L);
if NFFT<N
disp('NFFT is less than 2L')
break
else
N=NFFT;
end

Lx=length(x);
wd=hanning(Lx);
Uwd=sum(wd.*wd)/Lx;
r=autofft(x.*wd,L);
w=parzen(2*L-1);
rw2=r.*w(L:2*L-1);
f=zeros(N,1);
f(1:L)=rw2;
f(N-L+2:N)=flipud(rw2(2:L));
F=fft(f,N); % F is real beacuse f is even
F2=abs(F(1:N/2+1));
R=F2/Uwd;break
% Callibrate to deal with positive frequencies only

R=F2(1:N/2)+F2(2:N/2+1); R(N/2+1)=2*F2(N/2+1);
R=R/N;

varx=sum(R(1:N/2));