function [meanSE,kn,fi,fN,Serr_low,Serr_high]=espectro(VAR,nu)

%% ----------------------------- Spectra ----------------------------------

nA=length(VAR); % use length of shortest timeseries
% nu=10; % degrees of freedom
ns=nu/2; % number of segments
nn=nA/ns; % length of segments

m=1;
p=floor(nn); 

N=p; % number of data points in each segment
dt=1;  %1 hour, mean sampling time interval
dk=1/(dt*N);   % sampling wavenumber
n=1:(N/2);
kn=n./((N-1)*dt);

% Nyquist frequency
fN=1/(2*dt);

% Lowest resolved frequency
fi=1/(N*dt);

    % Hanning window
    wind=hann(N);
    M=floor(N/2);
    
% Detrend Data
WSPDis=VAR';

for j=1:ns % ns: number of segments
      
    am(:,j)=fft(WSPDis(m:p).*wind.*sqrt(8/3)); % windspeed
    SE(:,j)=2*(abs(am(2:(N/2)+1,j)).^2)/dk/N^2;
    varcheck(j,1)=var(WSPDis);
    varcheck(j,2)=sum(SE(:,j)*dk); 

    m=m+floor(nn);
    p=p+floor(nn);
end

% Compute Spectral Estimates (averaging over segments)
meanSE=mean(SE,2); % Mean Spectral Estimate

% Compute Spectral error bar for nu d.0.f.
alpha=0.05; % for 95% confidence
Serr_low=nu/chi2inv(alpha/2,nu);
Serr_high=nu/chi2inv(1-alpha/2,nu);

% If using the hanning window, nu is nu effective
g2=sum(wind.^2);
nueff=(8/3)*nA/M; % from Emery & Thompson Table 5.6.4, or 2/g2 where g2 is the sum of the squared weights
alpha=0.05; % for 95% confidence
Serr_low=nueff/chi2inv(alpha/2,nueff);
Serr_high=nueff/chi2inv(1-alpha/2,nueff);

% %% ------------------------------ Plots ----------------------------------
% figure(1)
% axes('Parent',gcf,'FontSize',18,'CLim',[1 2]);
% maximize;
% loglog(kn,meanSE,'k','linewidth',2)
% ylabel('Spectral Density [(\circ C)^2/cph]','FontSize',16)
% xlabel('Frequency [cph]','FontSize',16)
% hold on
% xlim([fi fN])
% set(gca,'XGrid','on');%,'XTick',[0.005 0.01 0.016 0.02 0.025 0.033 0.05 0.0833 0.1 0.11 0.125 0.14 0.16 0.2 0.25 0.33 0.5 1]);
% [hx,hy]=format_ticks_logscale(gca,{'$\frac{1}{200}$','$\frac{1}{144}$','$\frac{1}{100}$','$\frac{1}{72}$','$\frac{1}{60}$','$\frac{1}{48}$','$\frac{1}{40}$','$\frac{1}{30}$','$\frac{1}{24}$','$\frac{1}{20}$','$\frac{1}{12}$','$\frac{1}{10}$','$\frac{1}{9}$','$\frac{1}{8}$','$\frac{1}{7}$','$\frac{1}{6}$','$\frac{1}{5}$','$\frac{1}{4}$','$\frac{1}{3}$','$\frac{1}{2}$'},[],...
%     [0.005 0.0069 0.01 0.0139 0.016 0.0208 0.025 0.033 0.0417 0.05 0.0833 0.1 0.11 0.125 0.14 0.16 0.2 0.25 0.33 0.5]);
% 
% semilogy([0.33 0.33],[Serr_low Serr_high]*2,'color',[0.4 0.4 0.4],'linewidth',3)
% title([num2str(nu),' d.o.f, 95% confidence '],'fontsize',22)
% % legend('SST');
% grid on
% set(gcf,'PaperPositionMode','auto')
% set(get(gca,'XLabel'),'Position',[0.05 0.02 0])
% 
% hold on 

% print(['SpectraSST_nu',num2str(nu),'_SOTS_2012.eps'],'-depsc2')

