Dawn@4
|
1 function CPP = func_GetCPP(y, Fs, F0, variables)
|
Dawn@4
|
2 % CPP = func_GetCPP(y, Fs, F0, variables)
|
Dawn@4
|
3 % Input: y, Fs - from wavread
|
Dawn@4
|
4 % F0 - vector of fundamental frequencies
|
Dawn@4
|
5 % variables - global settings
|
Dawn@4
|
6 % Output: CPP vector
|
Dawn@4
|
7 % Notes: Calculates the CPP according to the formula from Hillenbrand, 1994.
|
Dawn@4
|
8 % Currently does not use textgrid information. For very very long files this
|
Dawn@4
|
9 % may be a process-bottleneck.
|
Dawn@4
|
10 %
|
Dawn@4
|
11 % Author: Yen-Liang Shue, Speech Processing and Auditory Perception Laboratory, UCLA
|
Dawn@4
|
12 % Copyright UCLA SPAPL 2009
|
Dawn@4
|
13
|
Dawn@4
|
14 N_periods = variables.Nperiods_EC;
|
Dawn@4
|
15 sampleshift = (Fs / 1000 * variables.frameshift);
|
Dawn@4
|
16
|
Dawn@4
|
17 CPP = zeros(length(F0), 1) * NaN;
|
Dawn@4
|
18
|
Dawn@4
|
19 N_ms = round(Fs / 1000); % Quefrency below 1ms are ignored as per Hillenbrand
|
Dawn@4
|
20
|
Dawn@4
|
21 for k=1:length(F0)
|
Dawn@4
|
22 ks = round(k * sampleshift);
|
Dawn@4
|
23
|
Dawn@4
|
24 if (ks <= 0 || ks > length(y))
|
Dawn@4
|
25 continue;
|
Dawn@4
|
26 end
|
Dawn@4
|
27
|
Dawn@4
|
28 F0_curr = F0(k);
|
Dawn@4
|
29
|
Dawn@4
|
30 if (isnan(F0_curr))
|
Dawn@4
|
31 continue;
|
Dawn@4
|
32 end
|
Dawn@4
|
33
|
Dawn@4
|
34 N0_curr = Fs / F0_curr;
|
Dawn@4
|
35
|
Dawn@4
|
36 ystart = round(ks - N_periods/2*N0_curr);
|
Dawn@4
|
37 yend = round(ks + N_periods/2*N0_curr) - 1;
|
Dawn@4
|
38
|
Dawn@4
|
39 if (ystart <= 0)
|
Dawn@4
|
40 continue;
|
Dawn@4
|
41 end;
|
Dawn@4
|
42
|
Dawn@4
|
43 if (yend > length(y))
|
Dawn@4
|
44 continue;
|
Dawn@4
|
45 end;
|
Dawn@4
|
46
|
Dawn@4
|
47 yseg = y(ystart:yend);
|
Dawn@4
|
48 yseg = yseg .* hamming(length(yseg)); %windowing
|
Dawn@4
|
49 YSEG = fft(yseg);
|
Dawn@4
|
50 yseg_c = ifft(log(abs(YSEG)));
|
Dawn@4
|
51 yseg_c_db = 10*log10(yseg_c .^ 2);
|
Dawn@4
|
52
|
Dawn@4
|
53 yseg_c_db = yseg_c_db(1:floor(length(yseg)/2));
|
Dawn@4
|
54 [vals, inx] = func_pickpeaks(yseg_c_db(N_ms:end)', 2*N0_curr);
|
Dawn@4
|
55 [pos, pinx] = min(abs(inx - N0_curr));
|
Dawn@4
|
56
|
Dawn@4
|
57 if (~isempty(pinx))
|
Dawn@4
|
58 inx = inx(pinx(1));
|
Dawn@4
|
59 vals = yseg_c_db(inx+N_ms-1);
|
Dawn@4
|
60 p = polyfit([N_ms:length(yseg_c_db)], yseg_c_db(N_ms:end)', 1);
|
Dawn@4
|
61 base_val = p(1) * (inx+N_ms-1) + p(2);
|
Dawn@4
|
62 CPP(k) = vals - base_val;
|
Dawn@4
|
63
|
Dawn@4
|
64 % figure;
|
Dawn@4
|
65 % plot(yseg_c_db); hold on;
|
Dawn@4
|
66 % plot(inx+N_ms-1, vals, 'rx');
|
Dawn@4
|
67 % plot([N_ms:length(yseg_c_db)], polyval(p,[N_ms:length(yseg_c_db)]), 'g');
|
Dawn@4
|
68
|
Dawn@4
|
69 end
|
Dawn@4
|
70
|
Dawn@4
|
71 end
|
Dawn@4
|
72
|