PitchTimeAuditory

computes the f_0 via a simple "auditory" approach
> called by ::ComputePitch
>
> @param x: audio signal
> @param iBlockLength: block length in samples
> @param iHopLength: hop length in samples
> @param f_s: sample rate of audio data 
>
> @retval f_0 fundamental frequency estimate (in Hz)
> @retval t time stamp of f_0 estimate (in s)
 ======================================================================

0001 %computes the f_0 via a simple "auditory" approach
0002 %> called by ::ComputePitch
0003 %>
0004 %> @param x: audio signal
0005 %> @param iBlockLength: block length in samples
0006 %> @param iHopLength: hop length in samples
0007 %> @param f_s: sample rate of audio data
0008 %>
0009 %> @retval f_0 fundamental frequency estimate (in Hz)
0010 %> @retval t time stamp of f_0 estimate (in s)
0011 % ======================================================================
0012 function [f_0, t] = PitchTimeAuditory(x, iBlockLength, iHopLength, f_s)
0013 
0014     % number of results
0015     iNumOfBlocks = ceil (length(x)/iHopLength);
0016     
0017     % compute time stamps
0018     t = ((0:iNumOfBlocks-1) * iHopLength + (iBlockLength/2)) / f_s;
0019     
0020     % allocate memory
0021     f_0 = zeros(1, iNumOfBlocks);
0022 
0023     %initialization
0024     f_max = 2000;
0025     eta_min = round(f_s/f_max);
0026     iNumBands = 20;
0027     fLengthLpInS = 0.001;
0028     iLengthLp = ceil(fLengthLpInS*f_s);
0029 
0030     % apply filterbank
0031     X = ToolGammatoneFb(x, f_s, iNumBands);
0032     
0033     % half wave rectification
0034     X(X < 0) = 0;
0035 
0036     % smooth the results
0037     b = ones(iLengthLp,1)/iLengthLp;
0038     X = filtfilt (b,1,X')';
0039         
0040     for n = 1:iNumOfBlocks
0041         afSumCorr = zeros(1, iBlockLength-1);
0042         i_start = (n-1)*iHopLength + 1;
0043         i_stop = min(length(x), i_start + iBlockLength - 1);
0044  
0045         % compute ACF per band and summarize
0046         for k = 1: iNumBands
0047             if sum(X(k, i_start:i_stop)) < 1e-20
0048                 continue;
0049             end
0050             afCorr = xcorr( [X(k, i_start:i_stop), ...
0051                             zeros(1, iBlockLength-(i_stop-i_start)-1)], ...
0052                             'coeff');
0053             afSumCorr = afSumCorr + afCorr((ceil((length(afCorr)/2))+1):end);
0054         end
0055         
0056         % find local maxima
0057         [afPeaks, eta_peak] = findpeaks(afSumCorr);
0058 
0059         i_tmp = find(eta_peak > eta_min,1);
0060         
0061         % set all entries below eta_min and the first local maximum to 0
0062         afSumCorr(1:eta_peak(i_tmp)-1) = 0;
0063 
0064         % calculate the acf maximum
0065         [fDummy, f_0(n)] = max(afSumCorr);
0066     end
0067  
0068     % convert to Hz
0069     ind = f_0 < eta_min;
0070     f_0 = f_s ./ f_0;
0071     f_0(ind==true) = 0;
0072 end