Line data Source code
1 : /*
2 : * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 : *
4 : * Use of this source code is governed by a BSD-style license
5 : * that can be found in the LICENSE file in the root of the source
6 : * tree. An additional intellectual property rights grant can be found
7 : * in the file PATENTS. All contributing project authors may
8 : * be found in the AUTHORS file in the root of the source tree.
9 : */
10 :
11 : #include <memory.h>
12 : #include <string.h>
13 : #ifdef WEBRTC_ANDROID
14 : #include <stdlib.h>
15 : #endif
16 : #include "pitch_estimator.h"
17 : #include "lpc_analysis.h"
18 : #include "codec.h"
19 :
20 :
21 :
22 0 : void WebRtcIsac_AllPoleFilter(double* InOut,
23 : double* Coef,
24 : size_t lengthInOut,
25 : int orderCoef) {
26 : /* the state of filter is assumed to be in InOut[-1] to InOut[-orderCoef] */
27 : double scal;
28 : double sum;
29 : size_t n;
30 : int k;
31 :
32 : //if (fabs(Coef[0]-1.0)<0.001) {
33 0 : if ( (Coef[0] > 0.9999) && (Coef[0] < 1.0001) )
34 : {
35 0 : for(n = 0; n < lengthInOut; n++)
36 : {
37 0 : sum = Coef[1] * InOut[-1];
38 0 : for(k = 2; k <= orderCoef; k++){
39 0 : sum += Coef[k] * InOut[-k];
40 : }
41 0 : *InOut++ -= sum;
42 : }
43 : }
44 : else
45 : {
46 0 : scal = 1.0 / Coef[0];
47 0 : for(n=0;n<lengthInOut;n++)
48 : {
49 0 : *InOut *= scal;
50 0 : for(k=1;k<=orderCoef;k++){
51 0 : *InOut -= scal*Coef[k]*InOut[-k];
52 : }
53 0 : InOut++;
54 : }
55 : }
56 0 : }
57 :
58 :
59 0 : void WebRtcIsac_AllZeroFilter(double* In,
60 : double* Coef,
61 : size_t lengthInOut,
62 : int orderCoef,
63 : double* Out) {
64 : /* the state of filter is assumed to be in In[-1] to In[-orderCoef] */
65 :
66 : size_t n;
67 : int k;
68 : double tmp;
69 :
70 0 : for(n = 0; n < lengthInOut; n++)
71 : {
72 0 : tmp = In[0] * Coef[0];
73 :
74 0 : for(k = 1; k <= orderCoef; k++){
75 0 : tmp += Coef[k] * In[-k];
76 : }
77 :
78 0 : *Out++ = tmp;
79 0 : In++;
80 : }
81 0 : }
82 :
83 :
84 0 : void WebRtcIsac_ZeroPoleFilter(double* In,
85 : double* ZeroCoef,
86 : double* PoleCoef,
87 : size_t lengthInOut,
88 : int orderCoef,
89 : double* Out) {
90 : /* the state of the zero section is assumed to be in In[-1] to In[-orderCoef] */
91 : /* the state of the pole section is assumed to be in Out[-1] to Out[-orderCoef] */
92 :
93 0 : WebRtcIsac_AllZeroFilter(In,ZeroCoef,lengthInOut,orderCoef,Out);
94 0 : WebRtcIsac_AllPoleFilter(Out,PoleCoef,lengthInOut,orderCoef);
95 0 : }
96 :
97 :
98 0 : void WebRtcIsac_AutoCorr(double* r, const double* x, size_t N, size_t order) {
99 : size_t lag, n;
100 : double sum, prod;
101 : const double *x_lag;
102 :
103 0 : for (lag = 0; lag <= order; lag++)
104 : {
105 0 : sum = 0.0f;
106 0 : x_lag = &x[lag];
107 0 : prod = x[0] * x_lag[0];
108 0 : for (n = 1; n < N - lag; n++) {
109 0 : sum += prod;
110 0 : prod = x[n] * x_lag[n];
111 : }
112 0 : sum += prod;
113 0 : r[lag] = sum;
114 : }
115 :
116 0 : }
117 :
118 :
119 0 : void WebRtcIsac_BwExpand(double* out, double* in, double coef, size_t length) {
120 : size_t i;
121 : double chirp;
122 :
123 0 : chirp = coef;
124 :
125 0 : out[0] = in[0];
126 0 : for (i = 1; i < length; i++) {
127 0 : out[i] = chirp * in[i];
128 0 : chirp *= coef;
129 : }
130 0 : }
131 :
132 0 : void WebRtcIsac_WeightingFilter(const double* in,
133 : double* weiout,
134 : double* whiout,
135 : WeightFiltstr* wfdata) {
136 : double tmpbuffer[PITCH_FRAME_LEN + PITCH_WLPCBUFLEN];
137 : double corr[PITCH_WLPCORDER+1], rc[PITCH_WLPCORDER+1];
138 : double apol[PITCH_WLPCORDER+1], apolr[PITCH_WLPCORDER+1];
139 0 : double rho=0.9, *inp, *dp, *dp2;
140 : double whoutbuf[PITCH_WLPCBUFLEN + PITCH_WLPCORDER];
141 : double weoutbuf[PITCH_WLPCBUFLEN + PITCH_WLPCORDER];
142 : double *weo, *who, opol[PITCH_WLPCORDER+1], ext[PITCH_WLPCWINLEN];
143 : int k, n, endpos, start;
144 :
145 : /* Set up buffer and states */
146 0 : memcpy(tmpbuffer, wfdata->buffer, sizeof(double) * PITCH_WLPCBUFLEN);
147 0 : memcpy(tmpbuffer+PITCH_WLPCBUFLEN, in, sizeof(double) * PITCH_FRAME_LEN);
148 0 : memcpy(wfdata->buffer, tmpbuffer+PITCH_FRAME_LEN, sizeof(double) * PITCH_WLPCBUFLEN);
149 :
150 0 : dp=weoutbuf;
151 0 : dp2=whoutbuf;
152 0 : for (k=0;k<PITCH_WLPCORDER;k++) {
153 0 : *dp++ = wfdata->weostate[k];
154 0 : *dp2++ = wfdata->whostate[k];
155 0 : opol[k]=0.0;
156 : }
157 0 : opol[0]=1.0;
158 0 : opol[PITCH_WLPCORDER]=0.0;
159 0 : weo=dp;
160 0 : who=dp2;
161 :
162 0 : endpos=PITCH_WLPCBUFLEN + PITCH_SUBFRAME_LEN;
163 0 : inp=tmpbuffer + PITCH_WLPCBUFLEN;
164 :
165 0 : for (n=0; n<PITCH_SUBFRAMES; n++) {
166 : /* Windowing */
167 0 : start=endpos-PITCH_WLPCWINLEN;
168 0 : for (k=0; k<PITCH_WLPCWINLEN; k++) {
169 0 : ext[k]=wfdata->window[k]*tmpbuffer[start+k];
170 : }
171 :
172 : /* Get LPC polynomial */
173 0 : WebRtcIsac_AutoCorr(corr, ext, PITCH_WLPCWINLEN, PITCH_WLPCORDER);
174 0 : corr[0]=1.01*corr[0]+1.0; /* White noise correction */
175 0 : WebRtcIsac_LevDurb(apol, rc, corr, PITCH_WLPCORDER);
176 0 : WebRtcIsac_BwExpand(apolr, apol, rho, PITCH_WLPCORDER+1);
177 :
178 : /* Filtering */
179 0 : WebRtcIsac_ZeroPoleFilter(inp, apol, apolr, PITCH_SUBFRAME_LEN, PITCH_WLPCORDER, weo);
180 0 : WebRtcIsac_ZeroPoleFilter(inp, apolr, opol, PITCH_SUBFRAME_LEN, PITCH_WLPCORDER, who);
181 :
182 0 : inp+=PITCH_SUBFRAME_LEN;
183 0 : endpos+=PITCH_SUBFRAME_LEN;
184 0 : weo+=PITCH_SUBFRAME_LEN;
185 0 : who+=PITCH_SUBFRAME_LEN;
186 : }
187 :
188 : /* Export filter states */
189 0 : for (k=0;k<PITCH_WLPCORDER;k++) {
190 0 : wfdata->weostate[k]=weoutbuf[PITCH_FRAME_LEN+k];
191 0 : wfdata->whostate[k]=whoutbuf[PITCH_FRAME_LEN+k];
192 : }
193 :
194 : /* Export output data */
195 0 : memcpy(weiout, weoutbuf+PITCH_WLPCORDER, sizeof(double) * PITCH_FRAME_LEN);
196 0 : memcpy(whiout, whoutbuf+PITCH_WLPCORDER, sizeof(double) * PITCH_FRAME_LEN);
197 0 : }
198 :
199 :
200 : static const double APupper[ALLPASSSECTIONS] = {0.0347, 0.3826};
201 : static const double APlower[ALLPASSSECTIONS] = {0.1544, 0.744};
202 :
203 :
204 0 : void WebRtcIsac_AllpassFilterForDec(double* InOut,
205 : const double* APSectionFactors,
206 : size_t lengthInOut,
207 : double* FilterState) {
208 : //This performs all-pass filtering--a series of first order all-pass sections are used
209 : //to filter the input in a cascade manner.
210 : size_t n,j;
211 : double temp;
212 0 : for (j=0; j<ALLPASSSECTIONS; j++){
213 0 : for (n=0;n<lengthInOut;n+=2){
214 0 : temp = InOut[n]; //store input
215 0 : InOut[n] = FilterState[j] + APSectionFactors[j]*temp;
216 0 : FilterState[j] = -APSectionFactors[j]*InOut[n] + temp;
217 : }
218 : }
219 0 : }
220 :
221 0 : void WebRtcIsac_DecimateAllpass(const double* in,
222 : double* state_in,
223 : size_t N,
224 : double* out) {
225 : size_t n;
226 : double data_vec[PITCH_FRAME_LEN];
227 :
228 : /* copy input */
229 0 : memcpy(data_vec+1, in, sizeof(double) * (N-1));
230 :
231 0 : data_vec[0] = state_in[2*ALLPASSSECTIONS]; //the z^(-1) state
232 0 : state_in[2*ALLPASSSECTIONS] = in[N-1];
233 :
234 0 : WebRtcIsac_AllpassFilterForDec(data_vec+1, APupper, N, state_in);
235 0 : WebRtcIsac_AllpassFilterForDec(data_vec, APlower, N, state_in+ALLPASSSECTIONS);
236 :
237 0 : for (n=0;n<N/2;n++)
238 0 : out[n] = data_vec[2*n] + data_vec[2*n+1];
239 :
240 0 : }
241 :
242 :
243 : /* create high-pass filter ocefficients
244 : * z = 0.998 * exp(j*2*pi*35/8000);
245 : * p = 0.94 * exp(j*2*pi*140/8000);
246 : * HP_b = [1, -2*real(z), abs(z)^2];
247 : * HP_a = [1, -2*real(p), abs(p)^2]; */
248 : static const double a_coef[2] = { 1.86864659625574, -0.88360000000000};
249 : static const double b_coef[2] = {-1.99524591718270, 0.99600400000000};
250 :
251 : /* second order high-pass filter */
252 0 : void WebRtcIsac_Highpass(const double* in,
253 : double* out,
254 : double* state,
255 : size_t N) {
256 : size_t k;
257 :
258 0 : for (k=0; k<N; k++) {
259 0 : *out = *in + state[1];
260 0 : state[1] = state[0] + b_coef[0] * *in + a_coef[0] * *out;
261 0 : state[0] = b_coef[1] * *in++ + a_coef[1] * *out++;
262 : }
263 0 : }
|
