function showspectrogram(y,fs)
% usage
%        showspectrogram(y,fs)
% 
% Shows the spectrogram and the time domain plot of the
% signal y.  If the sampling rate fs is not shows, fs=8192 is used.
% 
% R Kakarala, PhD
% U C Berkeley Extesion
% last rev: 5 Oct 2003
%

if nargin<2
    fs = 8192;
end;

% compute the spectrogram
% assume a 256 point sliding Hamming window, with an overlap of 128 samples
% is applied to the data.

Ly = length(y);
if (Ly<256)  % if the length < 256, add zeros to increase length to match window
    y(Ly+1:256)=0;
end;

L = 256;
Nfft = 256;
Noverlap = 256;  % this, the max overlap, is adjusted below
num_segs = 1025; 
while num_segs > 1024 
    Noverlap = Noverlap-10;
    shift = abs(L-Noverlap);
    num_segs = 1 + fix( (Ly-L)/shift );
end;

window = 0.54-0.46*cos(2*pi*(0:L-1)/(L-1));  % Hamming window
window = window(:);  y = y(:);  % make both column vectors to multiply term by term

% the following idea is from The DSP FIRST Toolbox, by McClellan, Shafer
% and Yoder. Basically, it divides the input into segments, each of length
% L and takes the windowed fft of each one
% 
B = zeros( Nfft/2+1, num_segs );     %- Pre-allocate the matrix
iseg = 0;
while( iseg<num_segs ) 
    nstart = 1 + iseg*shift;
    ysegw = window .* y( nstart:nstart+L-1);
    YF = fft( ysegw, Nfft );
    iseg = iseg + 1;
    B(:,iseg) = YF(1:Nfft/2+1);    
end

% use imagesc -- image scaled to show the histogram in dB
F = (0:(Nfft/2))/Nfft * fs;
T = ( L/2 + shift*(0:num_segs-1) ) / fs;
subplot(2,1,1);
imagesc(T,F,20*log10(abs(B)+eps)); axis xy; colormap(jet);
ylabel('frequency (Hz)'); xlabel('time (sec)'); title('Spectrogram (dB)');

% plot time signal, trying to align time axes for both plots
maxsampleused = nstart+L-2;
subplot(2,1,2);
plot((0:maxsampleused)/fs,y(1:maxsampleused+1)); 
axis([min(T),max(T),min(y),max(y)]); 
xlabel('time (sec)'); ylabel('signal');