	function [xpsd,freq,check,nwindows]=nonparametric_Spectragram(X,block,norm,window,deltat,errorcheck);
%
% function [xpsd,freq,check]=Nonparametric_Spectragram(X,block,norm,window,deltat,errorcheck);
%
% Here, E(X^2)=sum(psd*deltaf) is used as the normalization.
%
% This function determines the nonparametric spectrogram for a signal.  The signal
% is segmented into subwindows each block data points long.  A Fourier Transform is
% taken for each subwindow and the result is divided by each subwindows variance.  
% The power is multiplied by frequency.  Reduced frequency is not employeed.
%
% Inputs:
%
% X is a vector containing a real signal.   
%
% block is the size of desired subblock. block should be even.
%
% norm=1 will cause normalization of the psd at each subwindow.  NOT THE SAME AS NORMALIZING WITH 
% variance.
%
% window  When window=1, a cosine window is used.  When window=0,
%   no window is used.
%
% deltat is the sampling interval in seconds
%
% errorcheck  is used to test for errors  (no errorchecking is
%   currently implemented)
%
% Outputs:
%
% xpsd Power spectral density at frequencies in vector freq.
%
% Note, due to negative frequency info, that
%
%    sigma^2=2*sum(xpsd)*deltaf
%           =2*sum(xpsd)/deltat/length(signal)
%
% freq Frequencies for the xpsd vector.
%


if window==1;  % use cosine window    
  Ntau=block/2;
  windtemp=ones([2*Ntau 1]); 
  wind=ones([2*Ntau 1]);
  windtemp((floor(Ntau)-3):(floor(Ntau)+3))=zeros([7 1]);
  windtemp(floor(.8*Ntau):floor(Ntau))= ...
     cos(pi/(.4*Ntau)*(floor(.8*Ntau):floor(Ntau)));
  windtemp(floor(Ntau):floor(1.2*Ntau))= ...
     sin(pi/(.4*Ntau)*(floor(Ntau):floor(1.2*Ntau))-pi/2);
  
  wind(1:Ntau)=windtemp(Ntau:(block-1));
  wind((Ntau+1):block)=windtemp(1:Ntau);
else  % use no window
   wind=ones([block 1]);
end

% Determine the frequencies
freq=(0:(block/2-1))*(1/deltat)/block;

%Determine the number of subwindows employing a 50% offset approach
nwindows=floor(length(X)/(block/2))-1;

%Create a matrix full of zeroes of size (block X nwindows)to accomdate positive and negative frequencies
xpsd=zeros([block nwindows]);

%Iterate through the subwindows, calculate Fourier transform for each subwindow and normalize it with the variance.
for kkk=1:nwindows  
   Xwindow=X((1+(kkk-1)*block/2) : ((kkk-1)*block/2+block) );
   temp=fft(wind.*( Xwindow-mean(Xwindow) )  );
   xpsd(:,kkk)=conj(temp).*temp/block/var(Xwindow);
end


% Multiply by two accounting for on%ly using positive frequencies and normalizing, limit to positive frequencies
xpsd=2*xpsd(1:block/2,:)/block;
xpsd=xpsd*deltat*block;

%Check that the area under the individual PSD's equal 1
check=sum(xpsd(1:block/2,:)*(freq(2)-freq(1)));

%Multiply power by frequency
 for iii=1:(block/2)
   xpsd2(iii,:)=xpsd(iii,:)*freq(iii);
end

% if norm==1  % normalize psds to allow easier display
   for kk=1:nwindows; 
       xpsd(kk,:)=xpsd(kk,:)/sum(xpsd(kk,:))/(freq(2)-freq(1)); 
   end   
% end 

%Plotting routine
%figure(1);
%contour(1:273,log10(freq(5:625)),xpsd(5:625,1:273),15);
%set(gca, 'xtick',9:48:273);
%set(gca,'xticklabel',{'10:00','14:00','18:00','22:00','02:00','06:00'})
% [li col]=size(xpsd);
% t=0:(block/2)*0.1:length(X)*0.1-block*0.1;
% t=ones(li,1)*t;
% freq=freq';
% freq=freq*ones(1,col);

figure(1);
contourf(xpsd,15);

%set(gca, 'xtick',9:48:273);
%set(gca,'xticklabel',{'10:00','14:00','18:00','22:00','02:00','06:00'})
