if 1,
    MusPath = '/Users/sdubnov/Sounds/Music/';
    %fname = 'Bach__Prelude_in_G_major_Well_Tempered_Clavier_Book_I.wav'
    %fname = 'wtc1151.wav'
    %fname = 'Schumann__Piano_Quartet_in_E-flat_Major_I.wav'
    %fname = 'BEETHOVEN_SonF-1_Knorre_11025.wav'
    %fname = 'BEETHOVEN_SonF-1_Knorre.mp3'
    fname = 'BEETHOVEN_PianoSonataNo23_1_Steen-Nokleberg.mp3'
    %fname = 'BEETHOVEN_32variationsC-_Adamyan.mp3'

    fname = [MusPath fname];

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

    %[x,fs] = wavread(fname);
    %x = x(25000:end); %silence part in Beethoven SonF-1
    x = x + randn(size(x))*0.0001*max(max(x));

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

    newfs = 11025;
    if fs ~= newfs,
        x = resample(x, newfs, fs);
        fs = newfs;
    end

    close all

    %it seems that using just 60 lower channels of fft (no svd) and long
    %segments (200) worked fine
    %another option is to use cepstrum (32) with svd and shorter ta 20

    Ta = 100; %20; %20 %60; %200? for music (fft 60); %20;
    %SegSize = 3000; %2500; %3000; %1000; %3000;
    SegSize = 7000; %7000; %3000; %3000; %2000; %1000; %500; %Segments shorter the 3000 are too noisy
    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
    hopfactor = 4; %2, %4
    method = 21; %10; %11; %10; %20;

end

%note that CC from DoLamd can not be used directly, since here the
%coefficients CC are averaged into one vector per frame. See if this can be
%done so that CC from DoLamd would be used and here only the averaging would be done.

sndLen = length(x);
SegSamp = floor(SegSize*fs/1000);

hop = floor(SegSamp/hopfactor); %SegStep/4; /2
% segments for processing. each segment is SegSize second
SegLim = [(1:hop:sndLen-SegSamp)' (SegSamp+1:hop:sndLen)'];

disp([' Method ' int2str(method)])


NTa = floor(Ta*fs/1000); %number of samples for Ta msec.
if mod(NTa,2), %Nta odd
    NTa = NTa-1;
end

NSeg = size(SegLim,1);
CCm = zeros(p-1,NSeg);

for segnum = 1:NSeg, %size(SegLim,1),

    fprintf('Processing segment %d...out of %d\n',segnum,size(SegLim,1));
    % calculate Cepstral Coeff
    % ************************

    X = buffer(x(SegLim(segnum,1):SegLim(segnum,2)),NTa,NTa/2); %50 percent overlap between analysis frames
    X = X(:,1:end-1);

    if method == 10 || method == 11,
        FV = log(abs(fft(X))); %log(abs(fft(X))+0.0001); %log(abs(fft(X)));
        FV = FV(1:end/2+1,:);
    elseif method == 20 || method == 21,
        FV = rceps(X); %Do all cepstral analysis at once
    end
    I = find(isfinite(sum(FV)));     % find frames that are finite. (not infinite)
    if p > size(FV,1),
        p = size(FV,1);
    end
    FV = FV(:,I); %This causes time misalignment (missing frames)


    if method == 10, %fft
        % In case of fft, SVD is done on all data and then
        % reduced to p dimensions.        %
        [U,S,V] = svd(FV');
        us = FV'*V(:,1:p-1);  %using p dimensions (data reduction)
    elseif method == 20, %cepst (data reduction is by liftering)
        [U,S,V] = svd(FV(2:p,:)'); %no energy
        us = FV(2:p,:)'*V;  %should be equal to U*S (?)
    elseif method == 11 || method == 21, %no SVD
        us = FV(2:p,:)';
    end
    CCm(:,segnum) = mean(us)';
end
    fprintf('Done processing\n');


%Spectral Clustering: Singular Value Analysis of Markov matrix
%Clustering and Recurrence analysis is done by Eig. of the Markov Matrix.

% Note: the dot product distance works for cepstral coefficients, not
% log-Spectral-magnitude. If using methods 10, 11, probably need to use
% distpf

Smat = CCm'*CCm;
figure
imagesc([0:size(Smat,1)-1]*SegSize,[0:size(Smat,1)-1]*SegSize,Smat)
title('Similarity Matrix')

figure
imagesc([0:size(Smat,1)-1]*SegSize,[0:size(Smat,1)-1]*SegSize,(Smat*diag(1./sum(Smat))'))

%Note that here the normalization is "on the right", which means that the
%transition Markov matrix is oriented columnwise (columns sum to 1). So the
%probability vectors are "on the right" as well. This is the reason why we
%take matrix V as the probablity eigenvector.
%Note that in the presentation (Computer Audition and Ted's paper) this
%might be transposed.

%Similarity Matrix
invsqrtSmat = diag(1./sqrt(sum(Smat)));
CCmn = CCm * invsqrtSmat;
% Normalization by "node degree"
[U,L,V] = svd(CCmn);

%Zinv = diag(1./sum(Smat));
%Pmat = Smat*Zinv;
%[Vp,D] = eig(Pmat');
%-Vp(:,2) is equal to V1
% (check the math: is the svd trick above a correct way to estimate Pmat' ? 
% note that in sumulation Vp of both Pmat', Pmat and V1 seem very close)

Vav = conv(V(:,2)-mean(V(:,2)),ones(1,10)*0.1);
Vav = Vav(5:end-5);
Vav = Vav / std(Vav);

V1 = V(:,2);
%V1 = -V1;
V1 = V1/std(V1);
V1 = V1 - min(V1);

V2 = V(:,3);
V2 = V2/std(V2);
V2 = V2 - min(V2);


figure
plot(SegLim(:,1)/fs,Vav)
axis([0 SegLim(end,1)/fs -2.5 2.5])
legend('Spectral Clustering: 1st eigenvector')
xlabel('Seconds')

figure
specgram(x,[],fs)
hold on
plot(SegLim(:,1)/fs,V1*fs/2*0.2)
plot(SegLim(:,1)/fs,V2*fs/2*0.2,'k--')
legend('Spectral Clustering: 1st & 2nd eigenvector')
xlabel('Seconds')

break
% Here are the plots for the paper
%load BeethSon1
%load BeethSon23
%save 'BeethSon23.mat' Lant Ldis Sant Sdis V V1 V2 Smat fs SegLim SaT Tm Te E En ET fname
%load BeethSonF-1seg12
%save 'BeethSonF-1seg12.mat' Lant Ldis Sant Sdis V V1 V2 Smat fs SegLim SaT Tm Te E En ET fname
%load BeethSonF-1seg7Ta100.mat
%save 'BeethSonF-1seg7Ta100.mat' Lant Ldis Sant Sdis V V1 V2 Smat fs SegLim SaT Tm Te E En ET fname
%save 'Beeth32variationsC-seg7Ta100.mat' Lant Ldis Sant Sdis V V1 V2 Smat fs SegLim SaT Tm Te E En ET fname
%save 'BeethSon23seg7Ta100.mat' Lant Ldis Sant Sdis V V1 V2 Smat fs SegLim SaT Tm Te E En ET fname

figure
specgram(x,[],fs)
hold on
plot(SegSize/2/1000+ SegLim(:,1)/fs,V1*fs/2*0.15)
plot(Tm,Sant*fs/2*0.8,'-kd')
%plot(Tm-SegSize/hopfactor/2/1000,Sant*fs/2*0.8,'-kd')
%plot(SegLim(:,1)/fs,V2*fs/2*0.2,'k--')
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.');
TH = title('Beethoven, Sonata No.23 "Appassionata", 1st Mov.');
set(TH,'FontSize',18)
set(get(get(gcf,'CurrentAxes'),'XLabel'),'FontSize',14)
set(get(get(gcf,'CurrentAxes'),'YLabel'),'FontSize',14)


[B,A] = butter(4,0.3);
Santlp = filter(B,A,Sant);
Santlp = filter(B,A,Santlp(end:-1:1));
Santlp = Santlp(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
% end

len = min(length(V1lp),length(Santlp));
figure
plot3([1:len],V1lp(1:len),Santlp(1:len))

figure
specgram(x,[],fs)
hold on
plot(SegSize/2/1000+ SegLim(:,1)/fs,V1lp*fs/2*0.15)
plot(Tm,Santlp*fs/2*0.8,'-kd')

gplay(x,fs,fs,3)


