function yo = specgram(x,nfft,noverlap)
%SPECGRAM Calculate spectrogram from signal.
%	B = SPECGRAM(A,NFFT,NOVERLAP) calculates the spectrogram
%	for the signal in vector A.  The signal is split into
%	overlapping segments, each of which are windowed and then
%	Fourier transformed to produce an estimate of the
%	short-term frequency content of the signal.
%	The columns of B contain spectrum estimates as a function of time.
%	Time increases linearly across the columns, starting at column 1.
%	Frequency increases linearly down the rows, starting at 0.
%	If A is a length NX real-valued signal, B is a complex matrix
%	with NFFT rows and FIX((NX-NOVERLAP)/(NFFT-NOVERLAP)) columns.  If A
%	is real, only the first NFFT/2 rows are returned.
%	We use the Hanning window, which you can replace with any
%	other window of interest by editing this file.
%	B = SPECGRAM(A) produces the spectrogram of the signal A using a
%	window length which is the minimum of 256 and the length of the 
%	signal, and uses 1/2 the length of the window for the amount of
%	overlap.
%	SPECGRAM(...) plots the absolute value of the spectrogram, using
%	IMAGESC(B),AXIS XY, so the low frequency content of the first portion
%	of the signal is displayed in the lower left corner of the axes.

%       L. Shure 1-1-91
%       Copyright (c) 1991 by the MathWorks, Inc.

x = x(:); % make a column vector for ease later
nx = length(x);
if nargin == 1
	nfft = min(nx,256);
end
if nargin < 3
	noverlap = nfft/2;
end

% Change the following line to change the window used.
win = hanning(nfft);

% figure out number of columns for offsetting the signal
% this may truncate the last portion of the signal since we'd
% rather not append zeros unnecessarily - also makes the fancy
% indexing that follows more difficult.
ncol = fix((nx-noverlap)/(nfft-noverlap));
y = zeros(nfft,ncol);

% now stuff x into columns of y with the proper offset
% should be able to do this with fancy indexing!
colindex = 1 + (0:(ncol-1))*(nfft-noverlap);
rowindex = (1:nfft)';
y(:) = x(rowindex(:,ones(1,ncol))+colindex(ones(nfft,1),:)-1);
%y = x(rowindex(:,ones(1,ncol))+colindex(ones(nfft,1),:)-1);

% Apply the window to the array of offset signal segments.
y(:) = win(:,ones(1,ncol)).*y;

% now fft y which does the columns
y = fft(y);
if ~any(any(imag(x)))
	y = y(1:nfft/2,:);
end
% take abs, and use image to display results
if nargout == 0
	imagesc(abs(y));axis xy
else
	yo = y;
end
