function [spe,f]=fft_data(data,nfft,samp,opt);

% function [spe,f]=fft_data(data,nfft,samp,opt);
%
% Calculate spectrum with FFT
%   for real or complex number data series
%
%   data  : data (may be complex)
%   nfft   : point number
%   samp  : sampling rate [Hz] 
%   opt   : Option for figure plotting (Default: opt=0)
%             0 : Do not plot figure
%             1 : Plot figure
%
%         May 01, 2002, Masaki Ando
%         April 02, 2008, Masaki Ando

if nargin<=3,  opt=0; end

l=length(data);

av=fix(l/nfft);
      
    % Apply the window to the array of offset signal segments.
    window=hanning(nfft);  Be=8/3;
    %window=ones(nfft,1);  Be=1;
    y = reshape(data(1:nfft*av),nfft,av);
    y = window(:,ones(1,av)).*y;
    ys = fft(y,nfft);
    if av==1
      ls = ys.*sqrt(1 ./nfft./samp*2 *Be);  %%% for complex
    else
      ls = mean(ys')'.*sqrt( 1./nfft./samp*2 *Be); %%% for complex
    end;
    sp=ls(1:nfft/2,:);
    fp  = samp*(0:nfft/2-1)/nfft;
     if any(any(imag(data)))    % data are complex
        fm = samp*(-nfft/2+1:0)/nfft;
        sm=ls(nfft/2+1:nfft,:);
        spe=[sm;sp];
        f=[fm,fp]';
    else
        spe=sp;
        f=fp';
    end;
     
if opt==1
    if any(any(imag(data)))    % data are complex
      semilogy(f,spe)
    else
      loglog(f,spe)
    end;
   xlabel('Frequency [Hz]')
   ylabel('Amplitude')
   title([ 'Average : ',num2str(av),'   Band width : ',num2str(f(2)-f(1)),' Hz'])
   grid on
end


