1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | #if defined(HAVE_CONFIG_H1) |
30 | #include "config.h" |
31 | #endif |
32 | |
33 | #include "typedef.h" |
34 | #include "bv16cnst.h" |
35 | #include "bv16externs.h" |
36 | |
37 | int pitchtapquan( |
38 | Float *x, |
39 | int pp, |
40 | Float *b, |
41 | Float *re) |
42 | { |
43 | Float p[9], t, s0, s1, s2, cormax, cor; |
44 | Float t0, t1, t2; |
45 | Float *xt; |
46 | Float *fp0; |
47 | Float *fp1; |
48 | const Float *fp2; |
49 | int ppm2, qidx=0, i, j; |
50 | |
51 | ppm2 = pp - 2; |
52 | xt = x + XOFF(137 +1); |
53 | |
54 | for (i = 0; i < 3; i++) |
55 | { |
56 | fp0 = xt; |
57 | fp1 = x + XOFF(137 +1) - ppm2 - i - 1; |
58 | t = 0; |
59 | for (j = 0; j < FRSZ40; j++) |
60 | t += (*fp0++) * (*fp1++); |
61 | p[i] = t; |
62 | } |
63 | |
64 | fp0 = x+XOFF(137 +1)-ppm2-3; |
65 | s0 = *fp0++; |
66 | s1 = *fp0++; |
67 | s2 = *fp0--; |
68 | t0 = p[8] = s0*s0; |
69 | t1 = p[4] = s0*s1; |
70 | p[5] = s0*s2; |
71 | s0 = *fp0++; |
72 | s1 = *fp0++; |
73 | s2 = *fp0--; |
74 | t2 = s0*s0; |
75 | p[8] += t2; |
76 | p[4] += s0*s1; |
77 | p[5] += s0*s2; |
78 | for (i = 0; i < (FRSZ40 - 2); i++) |
79 | { |
80 | s0 = *fp0++; |
81 | s1 = *fp0++; |
82 | s2 = *fp0--; |
83 | p[8] += s0*s0; |
84 | p[4] += s0*s1; |
85 | p[5] += s0*s2; |
86 | } |
87 | s0 = *fp0++; |
88 | s1 = *fp0++; |
89 | s2 = *fp0--; |
| Value stored to 's2' is never read |
90 | p[7] = p[8] + (s0*s0) - t0; |
91 | p[3] = p[4] + (s0*s1) - t1; |
92 | p[6] = p[7] + (s1*s1) - t2; |
93 | |
94 | cormax = -1.0e30; |
95 | fp2 = bv16_pp9cb; |
96 | for (i = 0; i < PPCBSZ32; i++) |
97 | { |
98 | cor = 0.0; |
99 | fp1 = p; |
100 | for (j = 0; j < 9; j++) |
101 | cor += (*fp2++)*(*fp1++); |
102 | if (cor > cormax) |
103 | { |
104 | cormax = cor; |
105 | qidx = i; |
106 | } |
107 | } |
108 | |
109 | fp2 = bv16_pp9cb + qidx*9; |
110 | for (i = 0; i < 3; i++) |
111 | b[i] = fp2[i]*0.5; |
112 | |
113 | fp0 = x + XOFF(137 +1); |
114 | fp1 = x + XOFF(137 +1) - ppm2 - 1; |
115 | t = 0; |
116 | for (i = 0; i < FRSZ40; i++) |
117 | { |
118 | t0 = *fp0++ - b[0]*fp1[0] - b[1]*fp1[-1] - b[2]*fp1[-2]; |
119 | fp1++; |
120 | t += t0*t0; |
121 | } |
122 | *re = t; |
123 | return qidx; |
124 | } |
125 | |