function[peaks_pos, peaks_mag, peaks_pha]=stft(sig,winsize,winspace,peakspace) 
% [peaks_pos,peaks_mag,peaks_pha]=stft(sig,winsize,winspace,peakspace) -
% given a speech signal sampled at 8KHz, this function computes the stft at
% overlapping intervals.  These intervals are then windowed and the dft
% taken.  The peaks of each dft are computed, and these peaks and their
% positions (frequency in radians) are returned in the matrices peaks_pos &
% peaks_mag.  The phase of each peak is is returned in the matrix peaks_pha.
% The argument winspace is the spacing - in samples - between the peaks of
% sucessive windows.  The argument winsize is the size of the hamming window
% - in samples.  the argument peakspace is the size of the window in which
% the peaks function looks for local maxima. The columns of the returned
% matrices represent the frames which we manipulate. One frame corresponds to
% a column of amplitudes, phases and frequencies.

% the signal has some DC bias, so get rid of it by subtracting the mean
sig = sig - mean(sig);

% get the length of the signal
len = length(sig);
% set the hamming window size
% get the hamming window
hamwin = hamming(winsize);
% normalize the window so that the sum = 1.
hamwin = hamwin./sum(hamwin);
% set up the window counter; this also indexes the peaks matrices
windows = 1;
%set up a position counter
pos = 1;

% window the data and find all of the peaks 
% make sure we don't go past the end of the data
while (pos + winsize) < len,
	% window the data
	data = sig(pos:pos+winsize-1).*hamwin;
	% take the 512 point fft of the data
	data = fft(data,512);
	% put th magnitudes in DB representation; this makes it easier to
	% find the peaks
	data_mag = 20*log10(abs(data));
	data_pha = angle(data);
	% find the peaks of the windowed data
	% p, m, & z are column matrices
	[p,m,z]=peaks(data_mag(1:256),data_pha(1:256),peakspace);
	% convert the peak positions to frequecy in radians
	p = p.*(pi/256);
	% convert the magnitudes back to normal values
	m = 10.^(m/20);
	% put the positions and magnitudes of the peaks into 
	% their respective matrices
	peaks_pos(1:length(p),windows) = p;
	peaks_mag(1:length(m),windows) = m;
	peaks_pha(1:length(z),windows) = z;
	% increment the windows counter
	windows = windows + 1;	
	% adjust the position counter accroding to the overlap argument
	pos = pos + winspace;
end % of the while loop