correlogram.zip_correlogram_teager

function [y,ty]=correlogram(x,inc,nw,nlag,m,fs) % 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 memsize=voicebox('memsize'); % set memory size to use [nx,nc]=size(x); % number of sampes and channels if nargin<6 fs=1; if nargin<5 m='h'; if nargin<4 nlag=[]; if nargin<3 nw=[]; if nargin<2 inc=[]; end end end end end if ~numel(inc) inc=128; end if ~numel(nw) nw=inc; end nwin=length(nw); if nwin>1 win=nw(:); else nwin=nw; if any(m=='h') win=hamming(nwin); else win=ones(nwin,1); % window function end end if ~numel(nlag) nlag=nwin; end nwl=nwin+nlag-1; nt=pow2(nextpow2(nwl)); % transform length nf=floor((nx-nwl+inc)/inc); % number of frames i1=repmat((1:nwl)',1,nc)+repmat(0:nx:nx*nc-1,nwl,1); nb=min(nf,max(1,floor(memsize/(16*nwl*nc)))); % chunk size for calculating nl=ceil(nf/nb); % number of chunks jx0=nf-(nl-1)*nb; % size of first chunk in frames wincg=(1:nwin)*win.^2/(win'*win); fwin=conj(fft(win,nt,1)); % fft of window y=zeros(nlag,nc,nf); % first do partial chunk jx=jx0; x2=zeros(nwl,nc*jx); x2(:)=x(repmat(i1(:),1,jx)+repmat((0:jx-1)*inc,nwl*nc,1)); v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1)); w=real(ifft(fwin(:,ones(1,nc*jx)).*fft(x2.*conj(x2),nt,1))); w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:)); if isreal(x) y(:,:,1:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,jx); else y(:,:,1:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,jx); end % now do remaining chunks x2=zeros(nwl,nc*nb); for il=2:nl ix=jx+1; % first frame in this chunk jx=jx+nb; % last frame in this chunk x2(:)=x(repmat(i1(:),1,nb)+repmat((ix-1:jx-1)*inc,nwl*nc,1)); v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1)); w=real(ifft(fwin(:,ones(1,nc*nb)).*fft(x2.*conj(x2),nt,1))); w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:)); if isreal(x) y(:,:,ix:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,nb); else y(:,:,ix:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,nb); end end ty=(0:nf-1)'*inc+wincg; % calculate times of window centres if ~nargout imagesc(ty/fs,(1:nlag)/fs,squeeze(mean(y,2))); if nargin<6 us='samp'; else us='s'; end xlabel(['Time (' xticksi us ')']); ylabel(['Lag (' yticksi us ')']); axis 'xy'; title('Summary Correlogram'); end