function [IR,E,T,CCm] = IR_Analysis(x,fs,Ta,SegSize,p,q,method)
% IR = IR_Analysis(x,fs,Ta,SegSize,p,q)
% Perform and information rate analysis of a sound
% OUTPUT:
% IR = cumulative IR values for 1 to (p-1) principal components
% E = energy (for cepstrum), DC in case of fft (unused).
% T = times of the IR and E analysis points in sec. (middle of macro-frame) 
% CCm = mean feature vector of every segment (fft or cepstrum), without the energy E.
% INPUT:
% x = input sound vector
% fs = sampling frequency
% Ta = lenght in msec. of the  "timbre atom" (size of the cepstral analysis frame)
% SegSize = length in msec. of the IR analysis frame (macro-frame)
% p = smoothness of the spectral envelope (number of cepstral coefficients)
% or data reduction. See explanation below.
% q = power spectral density resolution (PSD) for IR estimation
% method = 10,11 - magnitude spectrum (dB) with or without SVD; 20,21 - cepstrum with or without SVD (respectively). 
%
% Usage of the p parameter depends on the method:
% In case of fft, SVD is done on all fft data and then reduced to p dimensions. 
% In case of cepstrum, p is the liftering coefficient and SVD is done on all p coefficients, with no further data reduction.
% The reason for this method is that using p cepstral coefficient is usually sufficient for data reduction.
%
% Copyright (c) Shlomo Dubnov sdubnov@ucsd.edu

PLT = 0;

if nargin == 2,
    Ta = 20; %200 %analysis frame (analysis step is Ta/3).
    SegSize = 1000;
    p = 32; %23; %120 %Adiel: 23;%120; %30; %10; %number of coeff.
    q = 64; %psd resolution
    method = 20;
elseif nargin == 3,
    SegSize = 1000;
    p = 32; %23; %120 %Adiel: 23;%120; %30; %10; %number of coeff.
    q = 64; %psd resolution
    method = 20;
elseif nargin == 4,
    p = 32; %23; %120 %Adiel: 23;%120; %30; %10; %number of coeff.
    q = 64; %psd resolution
    method = 20;
elseif nargin == 5,
    q = 64; %psd resolution
    method = 20;
elseif nargin == 6,
    method = 20;
end

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

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

CCm = zeros(p-1,size(Seg,1));

disp([' Method ' int2str(method)])
% Loop over all segments of the piece and extract sfm for each segment
segtot = size(Seg,1);
progressbar(0);
for segnum = 1:segtot,
    
    progressbar(segnum/segtot)
    %fprintf('Processing segment %d...out of %d\n',segnum,size(Seg,1));
    % calculate Cepstral Coeff
    % ************************
    
    
    NTa = floor(Ta*fs/1000); %number of samples for Ta msec.
    if mod(NTa,2), %Nta odd
        NTa = NTa-1;
    end
    
    X = buffer(x(Seg(segnum,1):Seg(segnum,2)),NTa,NTa/2); %50 percent overlap between analysis frames
    X = X(:,1:end-1);
    if method == 10 || method == 11,
        CC = abs(fft(X)); %log(abs(fft(X))+0.0001); %log(abs(fft(X)));
        CC = CC(1:end/2+1,:);
    elseif method == 20 || method == 21,
        CC = rceps(X); %Do all cepstral analysis at once
    end
    I = find(isfinite(sum(CC)));     % find frames that are finite. (not infinite)
    if p > size(CC,1),
        p = size(CC,1);
    end
    CC = CC(:,I);
    
    %fprintf('Processing ...')
    if method == 10, %fft
        % In case of fft, SVD is done on all data and then
        % reduced to p dimensions.        %
        [U,S,V] = svd(CC');
        us = CC'*V(:,1:p-1);  %using p dimensions (data reduction)
    elseif method == 20, %cepst (data reduction is by liftering)
        [U,S,V] = svd(CC(2:p,:)'); %no energy
        us = CC(2:p,:)'*V;  %should be equal to U*S (?)
    elseif method == 11 || method == 21, %no SVD
        us = CC(2:p,:)';
    end
    us = us - ones(size(us,1),1)*mean(us);
    if PLT,
        figure(1)
        imagesc(us'), axis xy
    end

    E(segnum) = mean(CC(1,:)); %energy
    CCm(:,segnum) = mean(CC(2:p,:)');
    
    % Spectral flatness
    % fprintf('calculate sfm ...\n');
    for i = 1:size(us,2), %29,
            spf(segnum,i) =sfm(us(:,i)',q); %gsfm(us(:,i)',q); - gsfm does not work since filter transient is too long
    end
    if PLT,
        figure(2)
        plot(-0.5*log(spf(segnum,:)))
        title('IR (-0.5*log(SFM))')
        pause     
    end
    
    %spf is t*nCepstracCoeff 
    % cumulative sum of spf
   
end

IR = cumsum(-0.5*log(sort(spf,2))'); %./([1:p]'*ones(1,size(spf,1)));
T = mean(Seg')/fs;
