%Do Information Rate Analysis on a sound file

Path = '../Sounds/';
%fname = 'Bach__Prelude_in_G_major_Well_Tempered_Clavier_Book_I.mp3';
%fname = 'wtc1151.wav';
fname = 'BEETHOVEN_SonF-1_Knorre.mp3'

fname

if strcmp(fname(end-2:end),'mp3'),
    [x,fs] = mp3read([Path fname],[],1,4); %read whole file, 'mono', downsample by 4
    %note that this uses my modified function 
else
    [x,fs] = wavread([Path fname]);
end

if min(size(x)) == 2,
    x = sum(x')';
end

%remove silence from both ends
xmat = buffer(x,1024,512);
xrms = std(xmat);
snr = 60;
abovesnr = find(xrms > 10^(-snr/20));
from = (min(abovesnr)-1)*512 + 1;
to = max(abovesnr)*512;
%BEETHOVEN_SonF-1: 

x = x(from:to);
close all


Ta = 20; %60; %20; %50; %200? for music (fft 60); %20;
%SegSize = 3000; %2500; %3000; %1000; %3000;
SegSize = 3000; %7000; %3000 %600 jungle; %1000; %500;
p = 30; %30; %60; %20; %32; %5; %32; 60; %5 or 32 fpr cepst %60 at least for fft
%p = Ta*fs/1000; %32 fpr cepst %60 at least for fft
q = 8; %16; %8; %spectral resolution (burg) used for IR estimation;
method = 20; %20; %21; %10; %20;

x = x - mean(x);
snr = 60;
x = x + randn(size(x))*10^(-snr/20)*max(abs(x));

tic
[IR,E,TT,CC] = IR_Analysis(x,fs,Ta,SegSize,p,q,method);
toc
tic
[P1,V] = Recur_Analysis(CC,method);
toc

Sant = IR(end,:);

Sant = Sant(2:end-1);
Sant = Sant - min(Sant);
Sant = Sant / max(Sant);

En = E(2:end-1);
En = En - min(En);
En = En / max(En);

Z1 = P1;
Z1 = Z1(2:end-1);
Z1 = Z1 - min(Z1);
Z1 = Z1/max(Z1(10:end-10));


T = TT(2:end-1);

paramtxt = ...
    sprintf('Parameters: Ta=%d, SegSize=%d, p=%d, q=%d, method=%d',Ta,SegSize,p,q,method);


% U1 = -log(P1);
% U1 = U1(2:end-1);
% U1 = U1 - min(U1);
% U1 = U1/max(U1(10:end-10));
% 
% V1 = V(:,1); 
% %note that this is the first eigenvector since we are using a probability (Markov) matrix
% V1 = V1(2:end-1);
% V1 = V1 - min(V1);
% V1 = V1/max(V1(10:end-10));
% 
% 
% figure
% specgram(x,[],fs)
% hold on
% plot(T,U1*fs/2,'--')
% plot(T,V1*fs/2)
% plot(T,Sant*fs/2*0.8,'-kd')
% LH = legend('Spectral Clustering 1st eigenvector','V-IR Spectral Flatness');
% xlabel('Seconds')
% set(LH,'FontSize',14)
% %TH = title('Beethoven, Sonata Fmin No.1, 1st Mov.,IR2');
% %TH = title('Beethoven, Sonata No.23 "Appassionata", 1st Mov.');
% %TH = title('Beethoven, Sonata No.15 "Pastorale", 1st Mov.');
% TH = title(fname(1:end-4));
% set(TH,'FontSize',18)
% xlabel(paramtxt)
% set(get(get(gcf,'CurrentAxes'),'XLabel'),'FontSize',14)
% set(get(get(gcf,'CurrentAxes'),'YLabel'),'FontSize',14)


figure
specgram(x,[],fs)
hold on
plot(T,Z1*fs/2*0.8)
plot(T,Sant*fs/2*0.8,'-kd')
LH = legend('Spectral Clustering 1st eigenvector (Stationary Distribution)','V-IR Spectral Flatness');
xlabel('Seconds')
set(LH,'FontSize',14)
%TH = title('Beethoven, Sonata Fmin No.1, 1st Mov.,IR2');
%TH = title('Beethoven, Sonata No.23 "Appassionata", 1st Mov.');
%TH = title('Beethoven, Sonata No.15 "Pastorale", 1st Mov.');
TH = title(fname(1:end-4));
set(TH,'FontSize',18)
xlabel(paramtxt)
set(get(get(gcf,'CurrentAxes'),'XLabel'),'FontSize',14)
set(get(get(gcf,'CurrentAxes'),'YLabel'),'FontSize',14)

[B,A] = butter(8,0.5);
Santlp = filter(B,A,Sant);
Santlp = filter(B,A,Santlp(end:-1:1));
Santlp = Santlp(end:-1:1);
Z1lp = filter(B,A,Z1);
Z1lp = filter(B,A,Z1lp(end:-1:1));
Z1lp = Z1lp(end:-1:1);
% V1lp = filter(B,A,V1);
% V1lp = filter(B,A,V1lp(end:-1:1));
% V1lp = V1lp(end:-1:1);

% figure
% for i = 1:length(Santlp),
%     plot(V1lp(1:i),Santlp(1:i))
%     pause          v

len = min(length(Z1lp),length(Santlp));
figure
plot3(T(1:len-4),Z1lp(1:len-4),Santlp(1:len-4))
%plot3(T(2:110),Z1lp(2:110),Santlp(2:110))
grid

figure
specgram(x,[],fs)
hold on
%plot(T,V1lp*fs/2*0.8)
plot(T,Z1lp*fs/2*0.8)
plot(T,Santlp*fs/2*0.8,'-kd')
TH = title(fname(1:end-4));
set(TH,'FontSize',18)
LH = legend('Signal Recurrence','V-IR using SFM');
xlabel(paramtxt)
set(get(get(gcf,'CurrentAxes'),'XLabel'),'FontSize',14)
set(get(get(gcf,'CurrentAxes'),'YLabel'),'FontSize',14)
