% Usage: spectro(f, T, dT, overlap, 'windowtype')
% Given a vector "f" of length N representing a signal of analog time
% length "T" this computes a spectrogram for f using windows of length
% dT with overlap percentage given by "overlap".

function spectro(f, T, dT, overlap, window)

% Convert f to a row vector.
[n,N] = size(f);
if n==1
    f2 = f;
else
    f2 = f';
end
[n,N] = size(f2);

%Set up grayscale color map.
mappy = [(0:255)/256; (0:255)/256; (0:255)/256]';
colormap(mappy);

%Bookkeeping
m = floor(dT/T*N); %Window size
r = floor((100-overlap)*m/100); %Index stepsize through signal
H = 0.5*N/T; %Highest frequency
dt = r/N*T; %Step in terms of time

M = floor((N-m)/r); %Number of blocks we'll do
G = zeros(m,M+1); %Array to hold transforms

for j=0:M
    if strcmp(window,'rect')
        g = f2(r*j+1:r*j+m).*rectwindow(m,1,m);
    elseif strcmp(window,'tri')
        g = f2(r*j+1:r*j+m).*triwindow(m,1,m);
    elseif strcmp(window,'gauss')
        g = f2(r*j+1:r*j+m).*gausswindow(m,1,m,5.0);
    elseif strcmp(window,'hamming')
        g = f2(r*j+1:r*j+m).*hammingwindow(m,1,m);
    else %default to rectangular
        if j==0 display('Window not recognized---default to rectangular'); end;
        g = f2(r*j+1:r*j+m).*rectwindow(m,1,m);
    end; 
    Gg = fft(g);
    G(:,j+1) = Gg';
end;

ma = max(max(log(1+abs(G).^2)));

X = ([0:M]+0.5)*dt; Y = [H:-1:0]; 
image(X,Y,255*log(1+abs(G(floor(m/2):-1:1,:).^2))/ma);
axis xy;
ylabel('Hertz');
xlabel('seconds');