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: 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 General Public License as published by
0056 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0063 %
0064 %   You can obtain a copy of the GNU General Public License from
0065 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0066 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

Generated by m2html © 2003