function [fftR,fftAxe] = FFTR(data,timeSamp,fftLength,fftDim);

% [fftR,fftAxe] = FFTR(data,timeSamp,fftLength,fftDim)
% Computes and displays (default) the one-sided
% FFT of a REAL 'data', plus the frequency axis.
% 'timeSamp' = 1 (default)
% 'fftLength' equals data length (default) or specified
% length or the next power of two (fftLength = 'n')
% 'fftDim' applies the FFT operation across the
%  dimension DIM.
%
%   Author: Laurent C. Duval
%   Institution: IFP, Technology Department
%   Created: 05/07/2002
%   Modified: 05/05/2003


if ~isreal(data) 
   error('Input data is complex')
end
if nargin < 4
   fftDim = 1;
end

if nargin < 2
   timeSamp = 1;
end
freqSamp = 1/timeSamp;
freqNyq = freqSamp/2;
[dataLength,dataNb] = size(data);
% if (isempty(fftLength))
%        fftLength = dataLength;
% end
% Choose FFT length
% If not given, simply dataLength, else 'n' for the next
% power of two
if (nargin < 3) 
   fftLength = dataLength;
else
   if ~isnumeric(fftLength)
      switch lower(fftLength)
      case 'n'
         fftLength = 2^nextpow2(dataLength);
      otherwise
         error('Ambiguous ''fftLength'' argument!')
      end
   end
end
nbFreq = ceil((fftLength+1)/2);
dataFft = fft(data,fftLength,fftDim);
fftR = 2*abs(dataFft(1:nbFreq,:));
% Compensate for DC 2-factor
fftR(1,:) = fftR(1,:)/2;
% Compensate for last frequency 2-factor 
if ~rem(fftLength,2),
   fftR(nbFreq,:) = fftR(nbFreq,:)/2;
end
% Data length independency
fftR = fftR/dataLength;
fftAxe = (0:nbFreq-1)'*2*freqNyq/fftLength;
if nargout == 0
      semilogx(fftAxe,20*log10(fftR));axis tight;grid on
      xlabel('Freq. (Hz)')
      ylabel('Magnitude (dB)')
end