correlogram

 make correlogram,
 Usage:
        do_env=1; do_hp=1;                            % flags to control options
        [b,a,fx,bx,gd]=gammabank(25,fs,'',[80 5000]); % determine filterbank
        y=filterbank(b,a,s,gd);                       % filter signal s
        if do_env
            [bl,al]=butter(3,2*800/fs);
            y=filter(bl,al,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
        correlogram(y,round(10e-3*fs),round(16e-3*fs),round(12e-3*fs),'',fs);

 Inputs:
        x(*,chan)  is the output of a 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
              [d = subtract DC component]
              [n = unnormalized]
              [v = variable analysis window proportional to lag]
              [p = output the peaks only]
               'h' = Hamming window
        fs         sample freq (affects only plots)

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

0001 function [y,ty]=correlogram(x,inc,nw,nlag,m,fs)
0002 % make correlogram,
0003 % Usage:
0004 %        do_env=1; do_hp=1;                            % flags to control options
0005 %        [b,a,fx,bx,gd]=gammabank(25,fs,'',[80 5000]); % determine filterbank
0006 %        y=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,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 %        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
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 %              [d = subtract DC component]
0024 %              [n = unnormalized]
0025 %              [v = variable analysis window proportional to lag]
0026 %              [p = output the peaks only]
0027 %               'h' = Hamming window
0028 %        fs         sample freq (affects only plots)
0029 %
0030 % Outputs:
0031 %        y(lag,chan,frame) is correlogram. Lags are 1:nlag samples
0032 %        ty                is time of the window energy centre for each frame
0033 %                            x(n) is at time n
0034 
0035 memsize=voicebox('memsize');    % set memory size to use
0036 [nx,nc]=size(x);  % number of sampes and channels
0037 if nargin<6
0038     fs=1;
0039     if nargin<5
0040         m='h';
0041         if nargin<4
0042             nlag=[];
0043             if nargin<3
0044                 nw=[];
0045                 if nargin<2
0046                     inc=[];
0047                 end
0048             end
0049         end
0050     end
0051 end
0052 if ~numel(inc)
0053     inc=128;
0054 end
0055 if ~numel(nw)
0056     nw=inc;
0057 end
0058 nwin=length(nw);
0059 if nwin>1
0060     win=nw(:);
0061 else
0062     nwin=nw;
0063     if any(m=='h')
0064         win=hamming(nwin);
0065     else
0066         win=ones(nwin,1);               % window function
0067     end
0068 end
0069 if ~numel(nlag)
0070     nlag=nwin;
0071 end
0072 nwl=nwin+nlag-1;
0073 nt=pow2(nextpow2(nwl));  % transform length
0074 nf=floor((nx-nwl+inc)/inc);  % number of frames
0075 i1=repmat((1:nwl)',1,nc)+repmat(0:nx:nx*nc-1,nwl,1);
0076 nb=min(nf,max(1,floor(memsize/(16*nwl*nc))));    % chunk size for calculating
0077 nl=ceil(nf/nb);                  % number of chunks
0078 jx0=nf-(nl-1)*nb;                % size of first chunk in frames
0079 wincg=(1:nwin)*win.^2/(win'*win);
0080 fwin=conj(fft(win,nt,1)); % fft of window
0081 y=zeros(nlag,nc,nf);
0082 % first do partial chunk
0083 jx=jx0;
0084 x2=zeros(nwl,nc*jx);
0085 x2(:)=x(repmat(i1(:),1,jx)+repmat((0:jx-1)*inc,nwl*nc,1));
0086 v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1));
0087 w=real(ifft(fwin(:,ones(1,nc*jx)).*fft(x2.*conj(x2),nt,1)));
0088 w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:));
0089 if isreal(x)
0090     y(:,:,1:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,jx);
0091 else
0092     y(:,:,1:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,jx);
0093 end
0094 % now do remaining chunks
0095 x2=zeros(nwl,nc*nb);
0096 for il=2:nl
0097     ix=jx+1;            % first frame in this chunk
0098     jx=jx+nb;           % last frame in this chunk
0099     x2(:)=x(repmat(i1(:),1,nb)+repmat((ix-1:jx-1)*inc,nwl*nc,1));
0100     v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1));
0101     w=real(ifft(fwin(:,ones(1,nc*nb)).*fft(x2.*conj(x2),nt,1)));
0102     w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:));
0103     if isreal(x)
0104         y(:,:,ix:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,nb);
0105     else
0106         y(:,:,ix:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,nb);
0107     end
0108 end
0109 ty=(0:nf-1)'*inc+wincg;       % calculate times of window centres
0110 if ~nargout
0111     imagesc(ty/fs,(1:nlag)/fs,squeeze(mean(y,2)));
0112     if nargin<6
0113         us='samp';
0114     else
0115         us='s';
0116     end
0117     xlabel(['Time (' xticksi us ')']);
0118     ylabel(['Lag (' yticksi us ')']);
0119     axis 'xy';
0120     title('Summary Correlogram');
0121 end
0122 
0123