function [yo,fo,to] = specg_data(x,nfft,Fs,window,noverlap)
%SPECGRAM Calculate spectrogram from signal.
%   [yo,fo,to] = specg_data(x,nfft,Fs,window,noverlap) 
%   calculates the spectrogram for 
%   the signal in vector x.  SPECGRAM splits the signal into overlapping 
%   segments, windows each with the WINDOW vector and forms the columns of
%   B with their zero-padded, length NFFT discrete Fourier transforms.  Thus
%   each column of B contains an estimate of the short-term, time-localized
%   frequency content of the signal A.  Time increases linearly across the 
%   columns of B, from left to right.  Frequency increases linearly down 
%   the rows, starting at 0.  If A is a length NX complex signal, B is a 
%   complex matrix with NFFT rows and 
%
%   Author(s): L. Shure, 1-1-91
%              T. Krauss, 4-2-93, updated
%   Copyright 1988-2001 The MathWorks, Inc.
%   $Revision: 1.7 $  $Date: 2001/04/02 20:20:51 $
%
%    reshaped by Masaki Ando
%                       2002/04/26

nx = length(x);
nwind = length(window);
if nx < nwind    % zero-pad x if it has length less than the window length
    x(nwind)=0;  nx=nwind;
end
x = x(:); % make a column vector for ease later
window = window(:); % be consistent with data set

nuni = nwind-noverlap;
ncol = fix((nx-noverlap)/nuni);
colindex = 1 + (0:(ncol-1))*nuni;

if nwind==nuni*8,
    y1=reshape(x(1:nuni*(ncol+7)),nuni,ncol+7);
    y2=[y1;[y1(:,2:ncol+7),zeros(nuni,1)]];
    y3=[y2;[y2(:,3:ncol+7),zeros(nuni*2,2)]];
    y4=[y3;[y3(:,5:ncol+7),zeros(nuni*4,4)]];
    y=y4(:,1:ncol);
elseif noverlap==0,   
    y=reshape(x(1:nuni*ncol),nuni,ncol);
elseif nwind==nuni*2,
    y1=reshape(x(1:nuni*ncol),nuni,ncol);
    y=[y1;[y1(:,2:ncol),x(nuni*ncol+1:nuni*(ncol+1))]];
elseif nwind==nuni*4,
    y1=reshape(x(1:nuni*ncol),nuni,ncol);
    y2=[y1;[y1(:,2:ncol),x(nuni*ncol+1:nuni*(ncol+1))]];
    y=[y2;[y2(:,3:ncol),x(nuni*ncol+1:nuni*(ncol+2)),x(nuni*(ncol+1)+1:nuni*(ncol+3))]];
else   
    rowindex = (1:nwind)';
    y = zeros(nwind,ncol);
    y(:) = x(rowindex(:,ones(1,ncol))+colindex(ones(nwind,1),:)-1);  
end;
    
    % Apply the window to the array of offset signal segments.
    y = window(:,ones(1,ncol)).*y;

    ys = fft(y,nfft);
    ls = ys.*conj(ys)./nfft./Fs*2;
    sp=ls(1:nfft/2,:);
    fp  = Fs*(0:nfft/2-1)/nfft;
     if any(any(imag(x)))    % x complex
        fm = Fs*(-nfft/2+1:0)/nfft;
        sm=ls(nfft/2+1:nfft,:);
        spe=[sm;sp];
        f=[fm,fp];
    else
        spe=sp;
        f=fp;
    end;
     
t = (colindex-1)'/Fs;

% take abs, and use image to display results
if nargout == 0
    newplot;
    if length(t)==1
 %       imagesc([0 1/f(2)],f,20*log10(abs(spe)+eps));axis xy; colormap(jet)
    else
  %      imagesc(t,f,20*log10(abs(spe)+eps));axis xy; colormap(jet)
    end
  %  xlabel('Time')
  %  ylabel('Frequency')
elseif nargout == 1,
    yo = spe;
elseif nargout == 2,
    yo = spe;
    fo = f;
elseif nargout == 3,
    yo = spe;
    fo = f;
    to = t;
end
