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