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.
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