# v_correlogram

## PURPOSE

V_CORRELOGRAM calculate correlogram [y,ty]=(x,inc,nw,nlag,m,fs)

## SYNOPSIS

function [y,ty]=v_correlogram(x,inc,nw,nlag,m,fs)

## DESCRIPTION

    V_CORRELOGRAM calculate correlogram [y,ty]=(x,inc,nw,nlag,m,fs)
     Usage:
            do_env=1; do_hp=1;                            % flags to control options
            [b,a,fx,bx,gd]=v_gammabank(25,fs,'',[80 5000]); % determine v_filterbank
            y=v_filterbank(b,a,s,gd);                       % filter signal s
            if do_env
                [bl,al]=butter(3,2*800/fs);
                y=filter(bl,al,v_teager(y,1),[],1);           % low pass envelope @ 800 Hz
            end
            if do_hp
                y=filter(fir1(round(16e-3*fs),2*64/fs,'high'),1,y,[],1);  % HP filter @ 64 Hz
            end
            v_correlogram(y,round(10e-3*fs),round(16e-3*fs),round(12e-3*fs),'',fs);

     Inputs:
            x(*,chan)  is the output of a filterbank (e.g. v_filterbank)
                       with one column per filter channel
            inc        frame increment in samples
            nw         window length in samples [or window function]
            nlag       number of lags to calculate
            m          mode string:
                   'h' = Hamming window
            fs         sample freq (affects only plots)

     Outputs:
            y(lag,chan,frame) is v_correlogram. Lags are 1:nlag samples
            ty                is time of the window energy centre for each frame
                                x(n) is at time n

     For each channel, the calculated correlation for frame n comprises
           y(t+1,*,n+1)=(win.*x(n*inc+(1:nw))'*x(n*inc+t+(1:nw))/sqrt(win'*abs(x(n*inc+(1:nw))).^2 * win'*abs(x(n*inc+t+(1:nw))).^2)
     This corresponds to the expression in (1.7) of [1] but incorporating the normalization from (1) of [2].

     Future planned mode options:
           'd' = subtract DC component
           'n' = unnormalized
           'v' = variable analysis window proportional to lag
           'p' = output the peaks only

     Refs:
     [1]    D. Wang and G. J. Brown. Fundamentals of computational auditory scene analysis.
           In D. Wang and G. Brown, editors, Computational Auditory Scene Analysis: Principles,
           Algorithms, and Applications, chapter 1. Wiley, Oct. 2006. doi: 10.1109/9780470043387.ch1
     [2]    S. Granqvist and B. Hammarberg. The correlogram: a visual display of periodicity. J. Acoust. Soc. Amer., 114: 2934, 2003.

## CROSS-REFERENCE INFORMATION

This function calls:

- v_cblabel V_CBLABEL add a label to a colorbar C=(L,H)
- v_colormap V_COLORMAP set and plot color map
- v_voicebox V_VOICEBOX set global parameters for Voicebox functions Y=(FIELD,VAL)
- v_windows V_WINDOWS Generate a standard windowing function (TYPE,N,MODE,P,OV)
- v_xticksi V_XTIXKSI labels the x-axis of a plot using SI multipliers S=(AH)
- v_yticksi V_YTIXKSI labels the y-axis of a plot using SI multipliers S=(AH)

This function is called by:

## SOURCE CODE

    0001 function [y,ty]=v_correlogram(x,inc,nw,nlag,m,fs)
    0002 %V_CORRELOGRAM calculate correlogram [y,ty]=(x,inc,nw,nlag,m,fs)
    0003 % Usage:
    0004 %        do_env=1; do_hp=1;                            % flags to control options
    0005 %        [b,a,fx,bx,gd]=v_gammabank(25,fs,'',[80 5000]); % determine v_filterbank
    0006 %        y=v_filterbank(b,a,s,gd);                       % filter signal s
    0007 %        if do_env
    0008 %            [bl,al]=butter(3,2*800/fs);
    0009 %            y=filter(bl,al,v_teager(y,1),[],1);           % low pass envelope @ 800 Hz
    0010 %        end
    0011 %        if do_hp
    0012 %            y=filter(fir1(round(16e-3*fs),2*64/fs,'high'),1,y,[],1);  % HP filter @ 64 Hz
    0013 %        end
    0014 %        v_correlogram(y,round(10e-3*fs),round(16e-3*fs),round(12e-3*fs),'',fs);
    0015 %
    0016 % Inputs:
    0017 %        x(*,chan)  is the output of a filterbank (e.g. v_filterbank)
    0018 %                   with one column per filter channel
    0019 %        inc        frame increment in samples
    0020 %        nw         window length in samples [or window function]
    0021 %        nlag       number of lags to calculate
    0022 %        m          mode string:
    0023 %               'h' = Hamming window
    0024 %        fs         sample freq (affects only plots)
    0025 %
    0026 % Outputs:
    0027 %        y(lag,chan,frame) is v_correlogram. Lags are 1:nlag samples
    0028 %        ty                is time of the window energy centre for each frame
    0029 %                            x(n) is at time n
    0030 %
    0031 % For each channel, the calculated correlation for frame n comprises
    0032 %       y(t+1,*,n+1)=(win.*x(n*inc+(1:nw))'*x(n*inc+t+(1:nw))/sqrt(win'*abs(x(n*inc+(1:nw))).^2 * win'*abs(x(n*inc+t+(1:nw))).^2)
    0033 % This corresponds to the expression in (1.7) of [1] but incorporating the normalization from (1) of [2].
    0034 %
    0035 % Future planned mode options:
    0036 %       'd' = subtract DC component
    0037 %       'n' = unnormalized
    0038 %       'v' = variable analysis window proportional to lag
    0039 %       'p' = output the peaks only
    0040 %
    0041 % Refs:
    0042 % [1]    D. Wang and G. J. Brown. Fundamentals of computational auditory scene analysis.
    0043 %       In D. Wang and G. Brown, editors, Computational Auditory Scene Analysis: Principles,
    0044 %       Algorithms, and Applications, chapter 1. Wiley, Oct. 2006. doi: 10.1109/9780470043387.ch1
    0045 % [2]    S. Granqvist and B. Hammarberg. The correlogram: a visual display of periodicity. J. Acoust. Soc. Amer., 114: 2934, 2003.
    0046
    0047 %      Copyright (C) Mike Brookes 2011-2018
    0048 %      Version: $Id: v_correlogram.m 10867 2018-09-21 17:35:59Z dmb $
    0049 %
    0050 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0051 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0052 %
    0053 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0054 %   This program is free software; you can redistribute it and/or modify
    0055 %   it under the terms of the GNU Lesser General Public License as published by
    0056 %   the Free Software Foundation; either version 3 of the License, or
    0057 %   (at your option) any later version.
    0058 %
    0059 %   This program is distributed in the hope that it will be useful,
    0060 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0061 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0062 %   GNU Lesser General Public License for more details.
    0063 %
    0064 %   You can obtain a copy of the GNU Lesser General Public License from
    0065 %   https://www.gnu.org/licenses/ .
    0066 %
    0067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0068 memsize=v_voicebox('memsize');     % set memory size to use
    0069 [nx,nc]=size(x);                % number of sampes and channels
    0070 if nargin<6
    0071     fs=1;                       % default sample frequency is 1 Hz
    0072     if nargin<5
    0073         m='h';                  % default analysis window is Hamming
    0074         if nargin<4
    0075             nlag=[];
    0076             if nargin<3
    0077                 nw=[];
    0078                 if nargin<2
    0079                     inc=[];
    0080                 end
    0081             end
    0082         end
    0083     end
    0084 end
    0085 if ~numel(inc)
    0086     inc=128;                    % default frame hop is 128 samples
    0087 end
    0088 if ~numel(nw)
    0089     nw=inc;                     % default analysis window length is the fame increment
    0090 end
    0091 nwin=length(nw);
    0092 if nwin>1                       % nw specifies the window function explicitly
    0093     win=nw(:);
    0094 else                            % nw gives the window length
    0095     nwin=nw;
    0096     if any(m=='h')
    0097         win=v_windows(3,nwin)'; % Hamming window
    0098     else
    0099         win=ones(nwin,1);       % window function
    0100     end
    0101 end
    0102 if ~numel(nlag)
    0103     nlag=nwin;
    0104 end
    0105 nwl=nwin+nlag-1;
    0106 nt=pow2(nextpow2(nwl));         % transform length
    0107 nf=floor((nx-nwl+inc)/inc);     % number of frames
    0108 i1=repmat((1:nwl)',1,nc)+repmat(0:nx:nx*nc-1,nwl,1);
    0109 nb=min(nf,max(1,floor(memsize/(16*nwl*nc))));    % chunk size for calculating
    0110 nl=ceil(nf/nb);                  % number of chunks
    0111 jx0=nf-(nl-1)*nb;                % size of first chunk in frames
    0112 wincg=(1:nwin)*win.^2/(win'*win); % determine window centre of energy
    0113 fwin=conj(fft(win,nt,1));       % conjugate fft of window
    0114 y=zeros(nlag,nc,nf);
    0115 % first do partial chunk
    0116 jx=jx0;
    0117 x2=zeros(nwl,nc*jx);
    0118 x2(:)=x(repmat(i1(:),1,jx)+repmat((0:jx-1)*inc,nwl*nc,1));
    0119 % the next line was previously: v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1));
    0120 v=ifft(conj(fft(x2(1:nwin,:).*repmat(win(:),1,nc*jx),nt,1)).*fft(x2,nt,1));                 % v(1:nlag,:) contains second half of xcorr(x2) output
    0121 w=max(real(ifft(fwin(:,ones(1,nc*jx)).*fft(x2.*conj(x2),nt,1))),0); % v(1:nlag,:) contains second half of xcorr(|x2|^2,win) output
    0122 w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:));
    0123 if isreal(x)
    0124     y(:,:,1:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,jx); % note: some values may be NaN is x=0 throughout the window
    0125 else
    0126     y(:,:,1:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,jx);
    0127 end
    0128 % now do remaining chunks
    0129 x2=zeros(nwl,nc*nb);
    0130 for il=2:nl
    0131     ix=jx+1;            % first frame in this chunk
    0132     jx=jx+nb;           % last frame in this chunk
    0133     x2(:)=x(repmat(i1(:),1,nb)+repmat((ix-1:jx-1)*inc,nwl*nc,1));
    0134     % the next line was previously: v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1));
    0135     v=ifft(conj(fft(x2(1:nwin,:).*repmat(win(:),1,nc*nb),nt,1)).*fft(x2,nt,1));
    0136     w=max(real(ifft(fwin(:,ones(1,nc*nb)).*fft(x2.*conj(x2),nt,1))),0);
    0137     w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:));
    0138     if isreal(x)
    0139         y(:,:,ix:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,nb);
    0140     else
    0141         y(:,:,ix:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,nb);
    0142     end
    0143 end
    0144 ty=(0:nf-1)'*inc+wincg;       % calculate times of window centres
    0145 if ~nargout
    0146     imagesc(ty/fs,(1:nlag)/fs,squeeze(mean(abs(y),2)));
    0147     if nargin<6
    0148         us='samp';
    0149     else
    0150         us='s';
    0151     end
    0152     xlabel(['Time (' v_xticksi us ')']);
    0153     ylabel(['Lag (' v_yticksi us ')']);
    0154     axis 'xy';
    0155     v_colormap('v_thermliny');
    0156     colorbar;
    0157     v_cblabel('Mean Correlation');
    0158     title('Summary Correlogram');
    0159 end
    0160
    0161