function [pav,f,tout,Pn]=fftpsdk(V,N,dT,OVLP,i)
%[P,f,t, Pn]=ftpsd(X,N,dT,OVLP,W);
%
% By G. Gonzalez
%
% Calculates and averages power spectrum P as a function of frequency F 
% from data in a vector X, taking N(=2^n) points for each record. 
% and overlaping them OVLP% . dT is the smapling time. 
% The program will use a Hanning window.  
% P=P(f,t) is a matrix with N/2 x T*100/OVLP elements.
% f is a frequency vector in Hertz with N/2 elements, and 
% t is a time vector in seconds with length(X)*dT*100/OVLP elements
% Pn is a power spectrum nromalized by the mean in each frequency bin
%
% Example : if  X is a vector of size 16384 x 1, with 8 seconds of a channel
%  with sampling frequency fs = 2048 (2048 Hz; 2^11 samples = 1sec), 
%  T = 2*fs = 2sec for each average 
%  dt = 1/fs;
% [P,f,t,Pn]=ftpsd(X,T,dt,50);
% produces:
%   a time vector t=[0 1 2 ...6]sec (Dt=1 sec = 50% of 2*fs samples) 
%   a frequency vector f = [0 0.5 1 ... 800] Hz (Df = 0.5 Hz = 1/(2*fs))
%   a matrix P(f,t) = 1601 x 7

last=max(size(V));
N1=round(N*800/2048)+1;
t=[0:N-1]*dT;
T=N*dT;


%Hanning window
w=2*sin(pi*t/T).^2;
bw=1.5/T;

if (size(w)~=size(V(1:N))), w=w'; end;

fs=N/T;
f=fs*[0:N1-1]/N;

inc=round((1-OVLP/100)*N);
k=N;
navg=0;

while k<=last
  y=fft(V(k-N+1:k).*w);
  navg=navg+1;
  pav(:,navg)=4*abs(y(1:N1)).^2/(2*bw*N^2);
  k=k+inc;
end;

f=f';

tout=[0:size(pav,2)-1]*inc/fs;

Pn = pav./(mean(pav')'*ones(1,size(pav,2)));

h(i) = figure('PaperPositionMode', 'auto');

subplot(221)
%avg psd
loglog(f,sqrt(mean(pav'))) % ASD
%loglog(f,mean(pav'))
xlabel('Hz')
ylabel('asd')
grid on
xlim([0.1 10000])

subplot(222)
%spectrogram
imagesc(tout,f,log10(pav))
axis xy
colorbar
xlabel('sec')
title('log10(P(f,t))')
ylabel('Hz')

subplot(223)
semilogx(f,std(pav')./mean(pav'))
grid on
ylabel('std/mean')
xlabel('Hz')
xlim([0.1 10000])

subplot(224)
imagesc(tout,f,log10(Pn),[-2 2])
axis xy
colorbar
title('log10(P/mean(P))')
xlabel('sec')
ylabel('Hz')

shg