1e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent/*
2e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *
4e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *  Use of this source code is governed by a BSD-style license
5e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *  that can be found in the LICENSE file in the root of the source
6e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *  tree. An additional intellectual property rights grant can be found
7e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *  in the file PATENTS.  All contributing project authors may
8e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent *  be found in the AUTHORS file in the root of the source tree.
9e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent */
10e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
11e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent/* Resamples a signal to an arbitrary rate. Used by the AEC to compensate for clock
12e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * skew by resampling the farend signal.
13e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent */
14e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
15c55a96383497a772a307b346368133960b02ad03Eric Laurent#include "aec_resampler.h"
16c55a96383497a772a307b346368133960b02ad03Eric Laurent
17e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent#include <assert.h>
18e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent#include <stdlib.h>
19e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent#include <string.h>
20e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent#include <math.h>
21e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
22e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent#include "aec_core.h"
23e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
24e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentenum { kEstimateLengthFrames = 400 };
25e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
26e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurenttypedef struct {
27c55a96383497a772a307b346368133960b02ad03Eric Laurent    short buffer[kResamplerBufferSize];
28e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float position;
29e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
30e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int deviceSampleRateHz;
31e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int skewData[kEstimateLengthFrames];
32e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int skewDataIndex;
33e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float skewEstimate;
34e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} resampler_t;
35e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
36e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentstatic int EstimateSkew(const int* rawSkew,
37e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                        int size,
38e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                        int absLimit,
39e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                        float *skewEst);
40e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
41e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentint WebRtcAec_CreateResampler(void **resampInst)
42e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{
43e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    resampler_t *obj = malloc(sizeof(resampler_t));
44e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    *resampInst = obj;
45e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    if (obj == NULL) {
46e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        return -1;
47e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
48e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
49e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    return 0;
50e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent}
51e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
52e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentint WebRtcAec_InitResampler(void *resampInst, int deviceSampleRateHz)
53e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{
54e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    resampler_t *obj = (resampler_t*) resampInst;
55e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    memset(obj->buffer, 0, sizeof(obj->buffer));
56e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    obj->position = 0.0;
57e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
58e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    obj->deviceSampleRateHz = deviceSampleRateHz;
59e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    memset(obj->skewData, 0, sizeof(obj->skewData));
60e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    obj->skewDataIndex = 0;
61e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    obj->skewEstimate = 0.0;
62e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
63e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    return 0;
64e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent}
65e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
66e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentint WebRtcAec_FreeResampler(void *resampInst)
67e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{
68e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    resampler_t *obj = (resampler_t*) resampInst;
69e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    free(obj);
70e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
71e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    return 0;
72e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent}
73e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
74e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentint WebRtcAec_ResampleLinear(void *resampInst,
75e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                             const short *inspeech,
76e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                             int size,
77e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                             float skew,
78e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                             short *outspeech)
79e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{
80e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    resampler_t *obj = (resampler_t*) resampInst;
81e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
82e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    short *y;
83e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float be, tnew, interp;
84e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int tn, outsize, mm;
85e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
86e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    if (size < 0 || size > 2 * FRAME_LEN) {
87e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        return -1;
88e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
89e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
90e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    // Add new frame data in lookahead
91e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    memcpy(&obj->buffer[FRAME_LEN + kResamplingDelay],
92e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent           inspeech,
93e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent           size * sizeof(short));
94e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
95e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    // Sample rate ratio
96e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    be = 1 + skew;
97e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
98e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    // Loop over input frame
99e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    mm = 0;
100e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    y = &obj->buffer[FRAME_LEN]; // Point at current frame
101e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
102e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    tnew = be * mm + obj->position;
103e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    tn = (int) tnew;
104e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
105e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    while (tn < size) {
106e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
107e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        // Interpolation
108e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        interp = y[tn] + (tnew - tn) * (y[tn+1] - y[tn]);
109e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
110e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        if (interp > 32767) {
111e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            interp = 32767;
112e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        }
113e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        else if (interp < -32768) {
114e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            interp = -32768;
115e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        }
116e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
117e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        outspeech[mm] = (short) interp;
118e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        mm++;
119e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
120e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        tnew = be * mm + obj->position;
121e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        tn = (int) tnew;
122e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
123e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
124e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    outsize = mm;
125e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    obj->position += outsize * be - size;
126e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
127e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    // Shift buffer
128e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    memmove(obj->buffer,
129e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            &obj->buffer[size],
130c55a96383497a772a307b346368133960b02ad03Eric Laurent            (kResamplerBufferSize - size) * sizeof(short));
131e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
132e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    return outsize;
133e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent}
134e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
135e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentint WebRtcAec_GetSkew(void *resampInst, int rawSkew, float *skewEst)
136e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{
137e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    resampler_t *obj = (resampler_t*)resampInst;
138e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int err = 0;
139e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
140e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    if (obj->skewDataIndex < kEstimateLengthFrames) {
141e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        obj->skewData[obj->skewDataIndex] = rawSkew;
142e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        obj->skewDataIndex++;
143e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
144e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    else if (obj->skewDataIndex == kEstimateLengthFrames) {
145e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        err = EstimateSkew(obj->skewData,
146e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                           kEstimateLengthFrames,
147e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                           obj->deviceSampleRateHz,
148e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                           skewEst);
149e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        obj->skewEstimate = *skewEst;
150e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        obj->skewDataIndex++;
151e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
152e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    else {
153e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        *skewEst = obj->skewEstimate;
154e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
155e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
156e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    return err;
157e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent}
158e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
159e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentint EstimateSkew(const int* rawSkew,
160c55a96383497a772a307b346368133960b02ad03Eric Laurent                 int size,
161c55a96383497a772a307b346368133960b02ad03Eric Laurent                 int deviceSampleRateHz,
162e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent                 float *skewEst)
163e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{
164e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    const int absLimitOuter = (int)(0.04f * deviceSampleRateHz);
165e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    const int absLimitInner = (int)(0.0025f * deviceSampleRateHz);
166e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int i = 0;
167e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int n = 0;
168e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float rawAvg = 0;
169e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float err = 0;
170e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float rawAbsDev = 0;
171e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int upperLimit = 0;
172e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    int lowerLimit = 0;
173e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float cumSum = 0;
174e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float x = 0;
175e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float x2 = 0;
176e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float y = 0;
177e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float xy = 0;
178e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float xAvg = 0;
179e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float denom = 0;
180e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    float skew = 0;
181e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
182e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    *skewEst = 0; // Set in case of error below.
183e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    for (i = 0; i < size; i++) {
184e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent      if ((rawSkew[i] < absLimitOuter && rawSkew[i] > -absLimitOuter)) {
185e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        n++;
186e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        rawAvg += rawSkew[i];
187e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent      }
188e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
189e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
190e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    if (n == 0) {
191e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent      return -1;
192e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
193e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    assert(n > 0);
194e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    rawAvg /= n;
195e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
196e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    for (i = 0; i < size; i++) {
197e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent      if ((rawSkew[i] < absLimitOuter && rawSkew[i] > -absLimitOuter)) {
198e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        err = rawSkew[i] - rawAvg;
199e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        rawAbsDev += err >= 0 ? err : -err;
200e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent      }
201e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
202e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    assert(n > 0);
203e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    rawAbsDev /= n;
204e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    upperLimit = (int)(rawAvg + 5 * rawAbsDev + 1); // +1 for ceiling.
205e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    lowerLimit = (int)(rawAvg - 5 * rawAbsDev - 1); // -1 for floor.
206e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
207e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    n = 0;
208e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    for (i = 0; i < size; i++) {
209e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        if ((rawSkew[i] < absLimitInner && rawSkew[i] > -absLimitInner) ||
210e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            (rawSkew[i] < upperLimit && rawSkew[i] > lowerLimit)) {
211e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            n++;
212e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            cumSum += rawSkew[i];
213e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            x += n;
214e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            x2 += n*n;
215e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            y += cumSum;
216e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent            xy += n * cumSum;
217e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        }
218e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
219e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
220e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    if (n == 0) {
221e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent      return -1;
222e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
223e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    assert(n > 0);
224e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    xAvg = x / n;
225e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    denom = x2 - xAvg*x;
226e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
227e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    if (denom != 0) {
228e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent        skew = (xy - xAvg*y) / denom;
229e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    }
230e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent
231e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    *skewEst = skew;
232e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent    return 0;
233e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent}
234