function [psd_final,t,f,f1,aver_FFT,check,u_mean,u_var]=fft_analysis(x,Fs,w,percent,av_freq,av_time)


% %The power spectral density, PSD, describes how the power (or variance) of a time series is distributed 
% %with frequency. Mathematically, it is defined as the Fourier Transform of the autocorrelation sequence 
% %of the time series. An equivalent definition of PSD is the squared modulus of the Fourier transform 
% %of the time series, scaled by a proper constant term
% 
% % Use next highest power of 2 greater than or equal to length(x) to calculate FFT. 
% NFFT= 2^(nextpow2(length(x))); 
% 
% 
% % Take fft, padding with zeros so that length(FFTX) is equal to 
% % NFFT 
% 
% FFTX = fft(x,NFFT); 
% 
% % Calculate the numberof unique points 
% NumUniquePts = ceil((NFFT+1)/2); 
% 
% % FFT is symmetric, throw away second half 
% FFTX = FFTX(1:NumUniquePts); 
% 
% % Take the magnitude of fft of x 
% MX = abs(FFTX); 
% 
% % Scale the fft so that it is not a function of the 
% % length of x 
% MX = MX/length(x); 
% 
% % Take the square of the magnitude of fft of x. 
% MX = MX.^2; 
% 
% % Multiply by 2 because you 
% % threw out second half of FFTX above 
% MX = MX*2; 
% 
% % DC Component should be unique. 
% MX(1) = MX(1)/2; 
% 
% %Note that running the above code fragment generates a simple mean-square (power) spectrum plot, not a power-spectral-density plot. 
% 
% % Nyquist component should also be unique.
% if ~rem(NFFT,2) 
%    % Here NFFT is even; therefore,Nyquist point is included. 
%    MX(end) = MX(end)/2; 
% end 
% 
% % This is an evenly spaced frequency vector with 
% % NumUniquePts points. 
% f = (0:NumUniquePts-1)*Fs/NFFT; 
% 
% % Generate the plot, title and labels. 
% plot(f,MX);
% %title('Power Spectrum of a 200Hz Sine Wave'); 
% xlabel('Frequency (Hz)'); 
% ylabel('Mean-squared Power'); 
%removing NaN from the data
u_NaN = find(~isnan(x));
x=x(u_NaN,:);
u_mean=[];
u_var=[];
% w=rectwin(length(x));
 NFFT=length(w);%2^(nextpow2(length(w)));

% if NFFT < 100
%     %disp('Data are not good');
%     psd_final=[];
%     t=[];
%     f=[];
%     f1=NaN;
%     aver_FFT=NaN;
%     check=[];
% else
%    NFFT=1024;  

w_squared=abs(w).^2;
lwin=length(w);
numoverlap=lwin*percent;
non_re_use=lwin-numoverlap;

k = fix((length(x)-numoverlap)/(length(w)-numoverlap))
FFT=[];
for i=1:k
     x_detrend=detrend(x(1+(i-1)*non_re_use :(i-1)*non_re_use+lwin),'linear');
     x_detrend=detrend(x_detrend,'constant');
     u_mean=[u_mean;mean(x(1+(i-1)*non_re_use :(i-1)*non_re_use+lwin))];
     u_var=[u_var;var(x(1+(i-1)*non_re_use :(i-1)*non_re_use+lwin))];
    xwin=w.*x_detrend;
    %xwin=w.*x(1+(i-1)*non_re_use :(i-1)*non_re_use+lwin );
    F=fft(xwin,NFFT);
    F=abs(F);
    F=F.^2;
    F=F./var(x_detrend);
    FFT=[F FFT];
end

aver_FFT=FFT;
 NumUniquePts=ceil((NFFT+1)/2);
 FFT=FFT(1:NumUniquePts,:); 
 f1=(0:NumUniquePts-1)*Fs/NFFT;
 f1=f1';
% for iii=1:601
%     FFT(iii,:)=FFT(iii,:)*f1(iii);
% end
 U=sum(w_squared)/lwin;
aver_FFT=mean(FFT,2);
  aver_FFT=2.*aver_FFT/(lwin*U);
  aver_FFT=aver_FFT./Fs;
  aver_FFT=aver_FFT.*(f1);

 aver_FFT(1)=aver_FFT(1)/2;
% semilogx(f1,aver_FFT);
% %title('Power Spectrum of a 200Hz Sine Wave'); 
% xlabel('Frequency (Hz)'); 
% ylabel('Power Spectral Density');
% %-----added for reviewers to remove noise-----%
 aver_FFT2=aver_FFT;
f2=f1;
aver_FFT=[];
f1=[];

L=5;
no_reused=1;
k=fix((length(f2)-(L-no_reused))/(L-(L-no_reused)));
for i=1:k
aver_FFT=[aver_FFT;mean(aver_FFT2(1+(i-1)*no_reused :(i-1)*no_reused+L))];
f1=[f1;mean(f2(1+(i-1)*no_reused :(i-1)*no_reused+L))];
end

%----end f reviewers stuff----%

%%%% for spectro
psd2=2.*FFT./(lwin*U);
psd2=psd2./Fs;
check=sum(psd2(1:NumUniquePts,:)*(f1(2)-f1(1)));
[li col]=size(psd2);
%psd2=psd2.*(f1*ones(1,col));


 t1=0:col-1;
if av_freq==0
     psd_final=psd2;
     f=f1;
    
else
    m=fix(li/av_freq);
    r=rem(li,av_freq);
    
    
    for i=1:m
        temp1=f1(1+av_freq*(i-1):av_freq*i,1);
        temp2=psd2(1+av_freq*(i-1):av_freq*i,:);
        temp1=mean(temp1);
        temp2=mean(temp2);
        f(i,1)=temp1;
        psd_final(i,:)=temp2;
    end
    
    if r~=0
        f(m+1,1)=mean(f1(m*av_freq+1:m*av_freq+r,1)); 
        psd_final(m+1,:)=mean(psd2(m*av_freq+1:m*av_freq+r,:));
    end
end



%moyennage horizontale
if av_time==0
     t=t1;
    
else 
    m=fix(col/av_time); 
    r=rem(col,av_time);
    for i=1:m
    temp1=t1(1,1+av_time*(i-1):av_time*i);
        temp2=psd_final(:,1+av_time*(i-1):av_time*i);
        temp1=mean(temp1);
        temp2=mean(temp2,2);
        t(1,i)=temp1;
        psd_final2(:,i)=temp2;
    end
    if r~=0
        t(1,m+1)=mean(t1(1,m*av_time+1:m*av_time+r)); 
        psd_final2(:,m+1)=mean(psd_final(:,m*av_time+1:m*av_time+r),2);
    end
    psd_final=psd_final2;
end
    [li col]=size(psd_final);
    t=ones(li,1)*t;
    f=f*ones(1,col);
    
    %end  
    
    
    % xwin=(x(1+(k-1)*numoverlap:length(x)));
% reste=lwin-length(xwin);
% z=zeros(reste,1);
% xwin=[xwin;z];
% xwin=w.*xwin;
% F=fft(xwin,NFFT);
% F=abs(F);
% F=F.^2;
% FFT=[FFT F];

