function [dyn] = pitchdynamics(y,sr)
%
% calculating a dynamics of a pitch to detect a music-like segment.
% y is a stretch of signal and sr is a sampling rate of it.
% kslee@ee.columbia.edu, sep. 2007
%

% Parameters
orderlpc = 12; % the order of LPC
nptslpc = 512; % every 32ms segment for LPC analysis

nwincorr = 1600; % 100ms window for autocorrelation
nptscorr = 80; % 5ms segment for autocorrelation
nlagscorr = 200; % maximum lag points for autocorrelation
higherlags = 100; % higher lags to keep for autocorrelation

nwinavgac = 3000; % number of frames for calculating averaged ACF (15-seconds)

% resampling input into 16kHz
sro = 16000;    
if sr ~= sro    
    gg = gcd(sro,sr); 
    y = resample(y,sro/gg,sr/gg);
    sr = sro;
end

% LPC Analysis using every 512 points to fit order of 12 LPC models.
[a,g,e] = lpcfit(y,orderlpc,nptslpc);

% half-wave rectfying and high-pass filtering of a residual signal
e = max(0,e);
e = filter([1 -1], [1 -.99],e);  
    
% autocorrelating it with 1600 points window hopped 80 points, 
% and are return as +/- 200 coefficients of autocorrelation
re = stxcorr(e',e',nwincorr,nptscorr,nlagscorr);
re = re(nlagscorr+1:end,:);


% Keep higher lag-points of AC, and normalize it.
nre = re(higherlags+1:end,:)./repmat(sqrt(sum(re(higherlags+1:end,:).^2)),higherlags+1,1);

% compensate it by long-duration averaged AC to remove periodic stationary noises
avgac = avgacf(nre,nwinavgac);
gainac = sum(nre.*avgac)./sum(avgac.^2);
nre = nre - avgac.*repmat(gainac,size(nre,1), 1);

% estimate the dyanmics of a pitch with cosine similarity measure
dyn = cosim(nre(:,2:end),nre(:,1:end-1));


% ----------------------------------------------------------------------------------------------------------
function [a] = avgacf(cor,w)
%
% average input with w-frames window hopped every frame
% kslee@ee.columbia.edu

[nlag,nfrm] = size(cor);

if nargin < 2
    w = nfrm;
end

if nfrm > w
    a = zeros(size(cor));
    w = min(w,nfrm-1);
    hw = floor(w/2);
    s = zeros(1,nfrm);

    for i = 1:nfrm
%        iw = i + [-hw:hw];
        iw = [[i-hw:i-1],[i+1:i+hw]];
        iw = iw(find(iw>=1 & iw <= nfrm));
        niw = length(iw);
        if i == 1
            a(:,i) = sum(cor(:,iw),2);
        else
            if iw(1) == 1 
                a(:,i) = a(:,i-1) + cor(:,iw(end));
            elseif iw(end) == nfrm
                a(:,i) = a(:,i-1) - cor(:,iw(1)-1);
            else
                a(:,i) = a(:,i-1) - cor(:,iw(1)-1) + cor(:,iw(end));
            end
  
        end
        s(1,i) = niw;
    end
    a = a./repmat(s,nlag,1);
    
else
    
    a = repmat(mean(cor,2),1,nfrm);
end


%----------------------------------------------------------------------------------------------------------
function s = cosim(x,y)
% calculating the cosine similarity between x and y
% kslee@ee.columbia.edu
%

[ndim,ndat] = size(x);
for nd = 1:ndat
    xx = x(:,nd);
    yy = y(:,nd);
    s(:,nd) = (xx'*yy)./sqrt(sum(xx.^2).*sum(yy.^2))';
end