
% KAIST CS891 Special Lecture on Music Information Retrieval (MIR)
% Lab1.4: Short-Time Fourier Transform (STFT) and time-frequency representation
% Kyogu Lee (August 2008)

clear all;
close all;

% Generate a sinusoid
fs = 44100; % sampling rate
amp = 1; % amplitude
freq = 440; % frequency in Hz
T = 1/fs; % sampling period
dur = .2; % duration in seconds
phi = 0; % phase in radian

% time-domain signal
x = amp*sin(2*pi*freq*[0:T:dur-T]+phi);

% block-processing (STFT)
M = 128; % block (window) size
hopsize = M/2; % hop size = M - overlap
nframe = floor(length(x)/hopsize); % # of frames
spec = []; % spectrogram

win = hann(M)'; % window
for i=1:nframe
  tt = (i-1)*hopsize+1:(i-1)*hopsize+M; % block index
  if tt(end) > length(x)
    xt = [x(tt(1):end),zeros(1,M-(length(x)-tt(1)+1))];
  else
    xt = x(tt);
  end
  xw = win.*xt; % windowing
  Xw = fft(xw);
  %Xw = fft(fftshift(xw));
  spec = [spec,Xw'];
end

[F,N] = size(spec); % # of freq./time frames
fi = 0:fs/F:fs-fs/F; % freq. resolution = fs/F
ti = 0:dur/N:dur-dur/N; % time resoultion = dur/N;

% display 2-d representation
subplot(211);
imagesc(ti,fi,20*log10(abs(spec))); ylim([0 fs/2]); colorbar;
set(gca,'YDir','normal');
title('2-D spectrogram');
ylabel('Frequency');

% Matlab function specgram() does all this for you
subplot(212);
specgram(x,M,fs,hann(M),M/2);