function [smooth,nval] = tsg_moveaverage(time, param, timeInterval, stdLim)
%
% Perform a moving average of a time series.
% The average is made over a period equal to 'timeinterval'
% 
% Data exceeding the average over that period by stdLim
% are not taken into account.
%
% Input
% time ............... date in Matlab serial number
% param............... parameter to smooth
% timeInterval ....... time interval used to compute the mean
% stdLim ............. Maximum standard deviation admitted
%
% Output
% smooth ............. time series para smoothed
% nval ............... nb values used to compute the mean in every interval
%
% time, para, smooth, nval are column vector of the same size
%
% Si pas de calcul possible valeur a NaN
%

% Memory allocation
% -----------------
smooth = zeros( size(param) );
nval   = zeros( size(param) );

% Loop over the param time series
% -------------------------------
for i = 1:length(param)
 
	% Select the param data over timeInterval interval
    % ind1 : indices of Param in the timeInterval
    % ind2 : indices of 'good' Param in the timeInterval
    % ---------------------------------------------------
    ind1 = find( time >= time(i) - timeInterval/2 & ...
                 time  <= time(i) + timeInterval/2 );
	ind2 = ind1;
    
    if ~isempty(ind2)
        
		currentStd   = Inf;
        previousStd = 0;
        
        while currentStd > stdLim && currentStd ~= previousStd
			
            previousStd = currentStd; 
            
            % Standard deviation and average over timeInterval
            % ------------------------------------------------
			currentStd = nanstd(  param(ind2) ); 
            meanParam  = nanmean( param(ind2) );
            
			% Indices of 'good' values of Param
            % ---------------------------------
            ind2 = find( param(ind1) >= meanParam - currentStd & ...
                         param(ind1) <= meanParam + currentStd);
            ind2 = ind1( ind2 );
        end
     
        smooth(i) = nanmean( param(ind2) );
    else
        smooth(i) = NaN;
    end
    nval( i ) = length( ind2 );
end