function [pout, fout, Ns] = myPSD(data,fs,Nfft)
% Computes the power spectral density of the data 
% which is recorded with the sample rate fs.  
% The data is cut into blocks of size Nfft. 
% It is assumed that the data is real so that the 
% negative spectrum is redundant and ignored. 
M = size(data,2);

Nh = Nfft/2;
L = length(data);
Ns = floor(L/Nh)-1;
%w = hamming(Nfft);
w = hann(Nfft);
pw = zeros(Nh,M);
for m = 1:M
for n = 0:Ns-1
    Noff = n*Nh;
    H = fft(w.*detrend(data(Noff+(1:Nfft),m)));
    pw(:,m) = pw(:,m) + abs(H(1:Nh)).^2;
end
end
df = fs/Nfft;

fout = df*((1:Nh)-1);
pout = pw/(Ns*Nfft);  
% In the above the Ns is for the mean and the 
% Nfft is for the scaling of the FFT.
end