
% KAIST CS891 Special Lecture on Music Information Retrieval (MIR)
% Lab2.1: Spectral low-level features
% Kyogu Lee (August 2008)

clear all;
close all;

% read audio file
%[x,fs] = wavread('../wav/T39-piano.wav'); % violin
%[x,fs] = wavread('../wav/drums_synt.wav'); % drums
%[x,fs] = wavread('../wav/T37-vibraphone.wav'); % vibraphone
[x,fs] = wavread('../wav/T08-violin.wav'); % violin
x = mean(x,2); % make mono

% analysis window size
M = 2^nextpow2(floor(fs*.1)); % 100ms long

% first compute spectrogram
[spec,fi,ti] = specgram(x,M,fs,hann(M),M/2);
nframe = length(ti);

% compute spectral centroid (SC) and spectral spread (SS)
for i=1:nframe
  X = spec(:,i); % X(k) at one frame
  sc(i) = sum(fi.*abs(X).^2)/sum(abs(X).^2); % spectral centroid
  ss(i) = sqrt(sum((fi-sc(i)).^2.*abs(X).^2)/sum(abs(X).^2)); % spectral spread
end


% compute spectral flatness (SF)
% SF is a subband-based feature (SF_b where b is a band)
% first compute filterbank frequencies using bark band (critical band)
bf = [0, 100, 200, 300, 400, 510, 630, 770, 920, 1080, 1270, 1480, 1720, 2000, 2320, 2700, 3150, ...
      3700, 4400, 5300, 6400, 7700, 9500, 12000, 15500, 20000]; % bark band frequencies
nband = length(bf)-1; % # of subbands

for b=1:nband
    bf_idx(b) = find(abs(bf(b+1)-fi)==min(abs(bf(b+1)-fi)));
end

for b=1:nband-1
    n_b(b) = bf_idx(b+1)-bf_idx(b)+1;
end
n_b = [bf_idx(1),n_b];

for b=1:nband
    if b==1, ff{b} = 1:bf_idx(b);
    else ff{b} = bf_idx(b-1)+1:bf_idx(b);
    end
end

% now calculate SF for each band
for i=1:nframe
    X = spec(:,i); % for each frame
    for b=1:nband
        sf(b,i) = (prod(abs(X(ff{b})).^2))^(1/n_b(b))/(sum(abs(X(ff{b})+0.00001).^2)/n_b(b)); % (geometric mean)/(arithmetic mean)
    end
end

sfm = mean(sf,1); % mean SF

% display
subplot(411); plot([1:length(x)]/fs,x); title('signal'); ylabel('x(t)');
subplot(412); plot(ti,sc); title('spectral centroid'); ylabel('SC(t)');
subplot(413); plot(ti,ss); title('spectral spread'); ylabel('SS(t)');
subplot(414); plot(ti,sfm); title('spectral flatness'); ylabel('SF(t)');

figure(2);
imagesc(ti,fi,20*log10(abs(spec))); set(gca,'YDir','normal');
