% spectrogram.m - compute fft spectrogram of X
 
W = hamming(512);        % 32 ms window at 16kHz SR
l = round(length(X)/256);
S = zeros(256,l);
Ph= zeros(256,l);
 
for i = 0:l-2
Y = fft(W .* X(256*i+1:256*i+512),512);
Pyy = (Y .* conj(Y)/512)';
S(:,i+1) = log(Pyy(1:256));
Ph(:,i+1)=angle(Y)(1:256);
end
 
f = 16000 * (0:255)/512;
surf(S)
 
xlabel('t (x 32 ms)')
ylim([1 256])
ylabel('f (x 256 Hz)')
zlabel('Power')