Commit f3ebf6b6 authored by David Henningsson's avatar David Henningsson
Browse files

lfe-filter: Import code from the Chrome OS audio server



The chrome OS audio server has some already existing code, which
has been made available under a BSD-style license, which should be
safe to import by us.
Signed-off-by: David Henningsson's avatarDavid Henningsson <david.henningsson@canonical.com>
parent ff329cda
......@@ -29,6 +29,9 @@ considered too small and stable to be considered as an external library) use the
more permissive MIT license. This include the device reservation DBus protocol
and realtime kit implementations.
A more permissive BSD-style license is used for LFE filters, see
src/pulsecore/filter/LICENSE.WEBKIT for details.
Additionally, a more permissive Sun license is used for code that performs
u-law, A-law and linear PCM conversions.
......
/*
* Copyright (C) 2010 Google Inc. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of
* its contributors may be used to endorse or promote products derived
* from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
* EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/* Copyright (c) 2013 The Chromium OS Authors. All rights reserved.
* Use of this source code is governed by a BSD-style license that can be
* found in the LICENSE file.
*/
/* Copyright (C) 2010 Google Inc. All rights reserved.
* Use of this source code is governed by a BSD-style license that can be
* found in the LICENSE.WEBKIT file.
*/
#include <math.h>
#include "biquad.h"
#ifndef max
#define max(a, b) ({ __typeof__(a) _a = (a); \
__typeof__(b) _b = (b); \
_a > _b ? _a : _b; })
#endif
#ifndef min
#define min(a, b) ({ __typeof__(a) _a = (a); \
__typeof__(b) _b = (b); \
_a < _b ? _a : _b; })
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
static void set_coefficient(struct biquad *bq, double b0, double b1, double b2,
double a0, double a1, double a2)
{
double a0_inv = 1 / a0;
bq->b0 = b0 * a0_inv;
bq->b1 = b1 * a0_inv;
bq->b2 = b2 * a0_inv;
bq->a1 = a1 * a0_inv;
bq->a2 = a2 * a0_inv;
}
static void biquad_lowpass(struct biquad *bq, double cutoff, double resonance)
{
/* Limit cutoff to 0 to 1. */
cutoff = max(0.0, min(cutoff, 1.0));
if (cutoff == 1) {
/* When cutoff is 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
} else if (cutoff > 0) {
/* Compute biquad coefficients for lowpass filter */
resonance = max(0.0, resonance); /* can't go negative */
double g = pow(10.0, 0.05 * resonance);
double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
double theta = M_PI * cutoff;
double sn = 0.5 * d * sin(theta);
double beta = 0.5 * (1 - sn) / (1 + sn);
double gamma = (0.5 + beta) * cos(theta);
double alpha = 0.25 * (0.5 + beta - gamma);
double b0 = 2 * alpha;
double b1 = 2 * 2 * alpha;
double b2 = 2 * alpha;
double a1 = 2 * -gamma;
double a2 = 2 * beta;
set_coefficient(bq, b0, b1, b2, 1, a1, a2);
} else {
/* When cutoff is zero, nothing gets through the filter, so set
* coefficients up correctly.
*/
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
}
}
static void biquad_highpass(struct biquad *bq, double cutoff, double resonance)
{
/* Limit cutoff to 0 to 1. */
cutoff = max(0.0, min(cutoff, 1.0));
if (cutoff == 1) {
/* The z-transform is 0. */
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
} else if (cutoff > 0) {
/* Compute biquad coefficients for highpass filter */
resonance = max(0.0, resonance); /* can't go negative */
double g = pow(10.0, 0.05 * resonance);
double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
double theta = M_PI * cutoff;
double sn = 0.5 * d * sin(theta);
double beta = 0.5 * (1 - sn) / (1 + sn);
double gamma = (0.5 + beta) * cos(theta);
double alpha = 0.25 * (0.5 + beta + gamma);
double b0 = 2 * alpha;
double b1 = 2 * -2 * alpha;
double b2 = 2 * alpha;
double a1 = 2 * -gamma;
double a2 = 2 * beta;
set_coefficient(bq, b0, b1, b2, 1, a1, a2);
} else {
/* When cutoff is zero, we need to be careful because the above
* gives a quadratic divided by the same quadratic, with poles
* and zeros on the unit circle in the same place. When cutoff
* is zero, the z-transform is 1.
*/
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
}
}
static void biquad_bandpass(struct biquad *bq, double frequency, double Q)
{
/* No negative frequencies allowed. */
frequency = max(0.0, frequency);
/* Don't let Q go negative, which causes an unstable filter. */
Q = max(0.0, Q);
if (frequency > 0 && frequency < 1) {
double w0 = M_PI * frequency;
if (Q > 0) {
double alpha = sin(w0) / (2 * Q);
double k = cos(w0);
double b0 = alpha;
double b1 = 0;
double b2 = -alpha;
double a0 = 1 + alpha;
double a1 = -2 * k;
double a2 = 1 - alpha;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
} else {
/* When Q = 0, the above formulas have problems. If we
* look at the z-transform, we can see that the limit
* as Q->0 is 1, so set the filter that way.
*/
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
}
} else {
/* When the cutoff is zero, the z-transform approaches 0, if Q
* > 0. When both Q and cutoff are zero, the z-transform is
* pretty much undefined. What should we do in this case?
* For now, just make the filter 0. When the cutoff is 1, the
* z-transform also approaches 0.
*/
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
}
}
static void biquad_lowshelf(struct biquad *bq, double frequency, double db_gain)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = max(0.0, min(frequency, 1.0));
double A = pow(10.0, db_gain / 40);
if (frequency == 1) {
/* The z-transform is a constant gain. */
set_coefficient(bq, A * A, 0, 0, 1, 0, 0);
} else if (frequency > 0) {
double w0 = M_PI * frequency;
double S = 1; /* filter slope (1 is max value) */
double alpha = 0.5 * sin(w0) *
sqrt((A + 1 / A) * (1 / S - 1) + 2);
double k = cos(w0);
double k2 = 2 * sqrt(A) * alpha;
double a_plus_one = A + 1;
double a_minus_one = A - 1;
double b0 = A * (a_plus_one - a_minus_one * k + k2);
double b1 = 2 * A * (a_minus_one - a_plus_one * k);
double b2 = A * (a_plus_one - a_minus_one * k - k2);
double a0 = a_plus_one + a_minus_one * k + k2;
double a1 = -2 * (a_minus_one + a_plus_one * k);
double a2 = a_plus_one + a_minus_one * k - k2;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
} else {
/* When frequency is 0, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
}
}
static void biquad_highshelf(struct biquad *bq, double frequency,
double db_gain)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = max(0.0, min(frequency, 1.0));
double A = pow(10.0, db_gain / 40);
if (frequency == 1) {
/* The z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
} else if (frequency > 0) {
double w0 = M_PI * frequency;
double S = 1; /* filter slope (1 is max value) */
double alpha = 0.5 * sin(w0) *
sqrt((A + 1 / A) * (1 / S - 1) + 2);
double k = cos(w0);
double k2 = 2 * sqrt(A) * alpha;
double a_plus_one = A + 1;
double a_minus_one = A - 1;
double b0 = A * (a_plus_one + a_minus_one * k + k2);
double b1 = -2 * A * (a_minus_one + a_plus_one * k);
double b2 = A * (a_plus_one + a_minus_one * k - k2);
double a0 = a_plus_one - a_minus_one * k + k2;
double a1 = 2 * (a_minus_one - a_plus_one * k);
double a2 = a_plus_one - a_minus_one * k - k2;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
} else {
/* When frequency = 0, the filter is just a gain, A^2. */
set_coefficient(bq, A * A, 0, 0, 1, 0, 0);
}
}
static void biquad_peaking(struct biquad *bq, double frequency, double Q,
double db_gain)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = max(0.0, min(frequency, 1.0));
/* Don't let Q go negative, which causes an unstable filter. */
Q = max(0.0, Q);
double A = pow(10.0, db_gain / 40);
if (frequency > 0 && frequency < 1) {
if (Q > 0) {
double w0 = M_PI * frequency;
double alpha = sin(w0) / (2 * Q);
double k = cos(w0);
double b0 = 1 + alpha * A;
double b1 = -2 * k;
double b2 = 1 - alpha * A;
double a0 = 1 + alpha / A;
double a1 = -2 * k;
double a2 = 1 - alpha / A;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
} else {
/* When Q = 0, the above formulas have problems. If we
* look at the z-transform, we can see that the limit
* as Q->0 is A^2, so set the filter that way.
*/
set_coefficient(bq, A * A, 0, 0, 1, 0, 0);
}
} else {
/* When frequency is 0 or 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
}
}
static void biquad_notch(struct biquad *bq, double frequency, double Q)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = max(0.0, min(frequency, 1.0));
/* Don't let Q go negative, which causes an unstable filter. */
Q = max(0.0, Q);
if (frequency > 0 && frequency < 1) {
if (Q > 0) {
double w0 = M_PI * frequency;
double alpha = sin(w0) / (2 * Q);
double k = cos(w0);
double b0 = 1;
double b1 = -2 * k;
double b2 = 1;
double a0 = 1 + alpha;
double a1 = -2 * k;
double a2 = 1 - alpha;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
} else {
/* When Q = 0, the above formulas have problems. If we
* look at the z-transform, we can see that the limit
* as Q->0 is 0, so set the filter that way.
*/
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
}
} else {
/* When frequency is 0 or 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
}
}
static void biquad_allpass(struct biquad *bq, double frequency, double Q)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = max(0.0, min(frequency, 1.0));
/* Don't let Q go negative, which causes an unstable filter. */
Q = max(0.0, Q);
if (frequency > 0 && frequency < 1) {
if (Q > 0) {
double w0 = M_PI * frequency;
double alpha = sin(w0) / (2 * Q);
double k = cos(w0);
double b0 = 1 - alpha;
double b1 = -2 * k;
double b2 = 1 + alpha;
double a0 = 1 + alpha;
double a1 = -2 * k;
double a2 = 1 - alpha;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
} else {
/* When Q = 0, the above formulas have problems. If we
* look at the z-transform, we can see that the limit
* as Q->0 is -1, so set the filter that way.
*/
set_coefficient(bq, -1, 0, 0, 1, 0, 0);
}
} else {
/* When frequency is 0 or 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
}
}
void biquad_set(struct biquad *bq, enum biquad_type type, double freq, double Q,
double gain)
{
/* Default is an identity filter. Also clear history values. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
bq->x1 = 0;
bq->x2 = 0;
bq->y1 = 0;
bq->y2 = 0;
switch (type) {
case BQ_LOWPASS:
biquad_lowpass(bq, freq, Q);
break;
case BQ_HIGHPASS:
biquad_highpass(bq, freq, Q);
break;
case BQ_BANDPASS:
biquad_bandpass(bq, freq, Q);
break;
case BQ_LOWSHELF:
biquad_lowshelf(bq, freq, gain);
break;
case BQ_HIGHSHELF:
biquad_highshelf(bq, freq, gain);
break;
case BQ_PEAKING:
biquad_peaking(bq, freq, Q, gain);
break;
case BQ_NOTCH:
biquad_notch(bq, freq, Q);
break;
case BQ_ALLPASS:
biquad_allpass(bq, freq, Q);
break;
case BQ_NONE:
break;
}
}
/* Copyright (c) 2013 The Chromium OS Authors. All rights reserved.
* Use of this source code is governed by a BSD-style license that can be
* found in the LICENSE file.
*/
#ifndef BIQUAD_H_
#define BIQUAD_H_
#ifdef __cplusplus
extern "C" {
#endif
/* The biquad filter parameters. The transfer function H(z) is (b0 + b1 * z^(-1)
* + b2 * z^(-2)) / (1 + a1 * z^(-1) + a2 * z^(-2)). The previous two inputs
* are stored in x1 and x2, and the previous two outputs are stored in y1 and
* y2.
*
* We use double during the coefficients calculation for better accurary, but
* float is used during the actual filtering for faster computation.
*/
struct biquad {
float b0, b1, b2;
float a1, a2;
float x1, x2;
float y1, y2;
};
/* The type of the biquad filters */
enum biquad_type {
BQ_NONE,
BQ_LOWPASS,
BQ_HIGHPASS,
BQ_BANDPASS,
BQ_LOWSHELF,
BQ_HIGHSHELF,
BQ_PEAKING,
BQ_NOTCH,
BQ_ALLPASS
};
/* Initialize a biquad filter parameters from its type and parameters.
* Args:
* bq - The biquad filter we want to set.
* type - The type of the biquad filter.
* frequency - The value should be in the range [0, 1]. It is relative to
* half of the sampling rate.
* Q - Quality factor. See Web Audio API for details.
* gain - The value is in dB. See Web Audio API for details.
*/
void biquad_set(struct biquad *bq, enum biquad_type type, double freq, double Q,
double gain);
#ifdef __cplusplus
} /* extern "C" */
#endif
#endif /* BIQUAD_H_ */
/* Copyright (c) 2013 The Chromium OS Authors. All rights reserved.
* Use of this source code is governed by a BSD-style license that can be
* found in the LICENSE file.
*/
#include "crossover.h"
#include "biquad.h"
static void lr4_set(struct lr4 *lr4, enum biquad_type type, float freq)
{
struct biquad q;
biquad_set(&q, type, freq, 0, 0);
lr4->b0 = q.b0;
lr4->b1 = q.b1;
lr4->b2 = q.b2;
lr4->a1 = q.a1;
lr4->a2 = q.a2;
lr4->x1 = 0;
lr4->x2 = 0;
lr4->y1 = 0;
lr4->y2 = 0;
lr4->z1 = 0;
lr4->z2 = 0;
}
/* Split input data using two LR4 filters, put the result into the input array
* and another array.
*
* data0 --+-- lp --> data0
* |
* \-- hp --> data1
*/
static void lr4_split(struct lr4 *lp, struct lr4 *hp, int count, float *data0,
float *data1)
{
float lx1 = lp->x1;
float lx2 = lp->x2;
float ly1 = lp->y1;
float ly2 = lp->y2;
float lz1 = lp->z1;
float lz2 = lp->z2;
float lb0 = lp->b0;
float lb1 = lp->b1;
float lb2 = lp->b2;
float la1 = lp->a1;
float la2 = lp->a2;
float hx1 = hp->x1;
float hx2 = hp->x2;
float hy1 = hp->y1;
float hy2 = hp->y2;
float hz1 = hp->z1;
float hz2 = hp->z2;
float hb0 = hp->b0;
float hb1 = hp->b1;
float hb2 = hp->b2;
float ha1 = hp->a1;
float ha2 = hp->a2;
int i;
for (i = 0; i < count; i++) {
float x, y, z;
x = data0[i];
y = lb0*x + lb1*lx1 + lb2*lx2 - la1*ly1 - la2*ly2;
z = lb0*y + lb1*ly1 + lb2*ly2 - la1*lz1 - la2*lz2;
lx2 = lx1;
lx1 = x;
ly2 = ly1;
ly1 = y;
lz2 = lz1;
lz1 = z;
data0[i] = z;
y = hb0*x + hb1*hx1 + hb2*hx2 - ha1*hy1 - ha2*hy2;
z = hb0*y + hb1*hy1 + hb2*hy2 - ha1*hz1 - ha2*hz2;
hx2 = hx1;
hx1 = x;
hy2 = hy1;
hy1 = y;
hz2 = hz1;
hz1 = z;
data1[i] = z;
}
lp->x1 = lx1;
lp->x2 = lx2;
lp->y1 = ly1;
lp->y2 = ly2;
lp->z1 = lz1;
lp->z2 = lz2;
hp->x1 = hx1;
hp->x2 = hx2;
hp->y1 = hy1;
hp->y2 = hy2;
hp->z1 = hz1;
hp->z2 = hz2;
}
/* Split input data using two LR4 filters and sum them back to the original
* data array.
*
* data --+-- lp --+--> data
* | |
* \-- hp --/
*/
static void lr4_merge(struct lr4 *lp, struct lr4 *hp, int count, float *data)
{
float lx1 = lp->x1;
float lx2 = lp->x2;
float ly1 = lp->y1;
float ly2 = lp->y2;
float lz1 = lp->z1;
float lz2 = lp->z2;
float lb0 = lp->b0;
float lb1 = lp->b1;
float lb2 = lp->b2;
float la1 = lp->a1;
float la2 = lp->a2;
float hx1 = hp->x1;
float hx2 = hp->x2;
float hy1 = hp->y1;
float hy2 = hp->y2;
float hz1 = hp->z1;
float hz2 = hp->z2;
float hb0 = hp->b0;
float hb1 = hp->b1;
float hb2 = hp->b2;
float ha1 = hp->a1;
float ha2 = hp->a2;
int i;
for (i = 0; i < count; i++) {
float x, y, z;
x = data[i];
y = lb0*x + lb1*lx1 + lb2*lx2 - la1*ly1 - la2*ly2;
z = lb0*y + lb1*ly1 + lb2*ly2 - la1*lz1 - la2*lz2;
lx2 = lx1;