Dawn@4: function CPP = func_GetCPP(y, Fs, F0, variables)
Dawn@4: % CPP = func_GetCPP(y, Fs, F0, variables)
Dawn@4: % Input:  y, Fs - from wavread
Dawn@4: %         F0 - vector of fundamental frequencies
Dawn@4: %         variables - global settings
Dawn@4: % Output: CPP vector
Dawn@4: % Notes:  Calculates the CPP according to the formula from Hillenbrand, 1994.
Dawn@4: % Currently does not use textgrid information. For very very long files this 
Dawn@4: % may be a process-bottleneck.
Dawn@4: %
Dawn@4: % Author: Yen-Liang Shue, Speech Processing and Auditory Perception Laboratory, UCLA
Dawn@4: % Copyright UCLA SPAPL 2009
Dawn@4: 
Dawn@4: N_periods = variables.Nperiods_EC;
Dawn@4: sampleshift = (Fs / 1000 * variables.frameshift);
Dawn@4: 
Dawn@4: CPP = zeros(length(F0), 1) * NaN;
Dawn@4: 
Dawn@4: N_ms = round(Fs / 1000); % Quefrency below 1ms are ignored as per Hillenbrand
Dawn@4: 
Dawn@4: for k=1:length(F0)
Dawn@4:     ks = round(k * sampleshift);
Dawn@4:     
Dawn@4:     if (ks <= 0 || ks > length(y))
Dawn@4:         continue;
Dawn@4:     end
Dawn@4:     
Dawn@4:     F0_curr = F0(k);
Dawn@4:     
Dawn@4:     if (isnan(F0_curr))
Dawn@4:         continue;
Dawn@4:     end
Dawn@4:     
Dawn@4:     N0_curr = Fs / F0_curr;
Dawn@4:   
Dawn@4:     ystart = round(ks - N_periods/2*N0_curr);
Dawn@4:     yend = round(ks + N_periods/2*N0_curr) - 1;
Dawn@4:     
Dawn@4:     if (ystart <= 0)
Dawn@4:         continue;
Dawn@4:     end;
Dawn@4:     
Dawn@4:     if (yend > length(y)) 
Dawn@4:         continue; 
Dawn@4:     end;
Dawn@4:   
Dawn@4:     yseg = y(ystart:yend);  
Dawn@4:     yseg = yseg .* hamming(length(yseg));  %windowing
Dawn@4:     YSEG = fft(yseg);
Dawn@4:     yseg_c = ifft(log(abs(YSEG)));
Dawn@4:     yseg_c_db = 10*log10(yseg_c .^ 2);
Dawn@4:   
Dawn@4:     yseg_c_db = yseg_c_db(1:floor(length(yseg)/2));
Dawn@4:     [vals, inx] = func_pickpeaks(yseg_c_db(N_ms:end)', 2*N0_curr);
Dawn@4:     [pos, pinx] = min(abs(inx - N0_curr));
Dawn@4:     
Dawn@4:     if (~isempty(pinx))
Dawn@4:         inx = inx(pinx(1));
Dawn@4:         vals = yseg_c_db(inx+N_ms-1);
Dawn@4:         p = polyfit([N_ms:length(yseg_c_db)], yseg_c_db(N_ms:end)', 1);        
Dawn@4:         base_val = p(1) * (inx+N_ms-1) + p(2);
Dawn@4:         CPP(k) = vals - base_val;
Dawn@4:         
Dawn@4: %     figure;
Dawn@4: %     plot(yseg_c_db); hold on;
Dawn@4: %     plot(inx+N_ms-1, vals, 'rx');
Dawn@4: %     plot([N_ms:length(yseg_c_db)], polyval(p,[N_ms:length(yseg_c_db)]), 'g');
Dawn@4:     
Dawn@4:     end
Dawn@4:     
Dawn@4: end
Dawn@4: