We are currently experiencing downtime impacting viewing & cloning the Mesa repo, and some GitLab pages returning 503. Please see #freedesktop on IRC for more updates.

adding resample lib

Original commit message from CVS:
adding resample lib
parent 9161ba1c
libdir = $(libdir)/gst
lib_LTLIBRARIES = libresample.la
if HAVE_CPU_I386
ARCHCFLAGS = -march=i486
else
if HAVE_CPU_PPC
ARCHCFLAGS = -Wa,-m7400
else
ARCHCFLAGS =
endif
endif
libresample_la_SOURCES = dtos.c functable.c resample.c resample.h
libresample_la_LIBADD = $(GST_LIBS)
libresample_la_CFLAGS = $(GST_CFLAGS) -ffast-math $(ARCHCFLAGS)
noinst_HEADERS = resample.h
#check_PROGRAMS = test
#test_SOURCES = test.c
#test_LDADD = libresample.la
This is a snapshot of my current work developing an audio
resampling library. While working on this library, I started
writing lots of general purpose functions that should really
be part of a larger library. Rather than have a constantly
changing library, and since the current code is capable, I
decided to freeze this codebase for use with gstreamer, and
move active development of the code elsewhere.
The algorithm used is based on Shannon's theorem, which says
that you can recreate an input signal from equidistant samples
using a sin(x)/x filter; thus, you can create new samples from
the regenerated input signal. Since sin(x)/x is expensive to
evaluate, an interpolated lookup table is used. Also, a
windowing function (1-x^2)^2 is used, which aids the convergence
of sin(x)/x for lower frequencies at the expense of higher.
There is one tunable parameter, which is the filter length.
Longer filter lengths are obviously slower, but more accurate.
There's not much reason to use a filter length longer than 64,
since other approximations start to dominate. Filter lengths
as short as 8 are audially acceptable, but should not be
considered for serious work.
Performance: A PowerPC G4 at 400 Mhz can resample 2 audio
channels at almost 10x speed with a filter length of 64, without
using Altivec extensions. (My goal was 10x speed, which I almost
reached. Maybe later.)
Limitations: Currently only supports streams in the form of
interleaved signed 16-bit samples.
The test.c program is a simple regression test. It creates a
test input pattern (1 sec at 48 khz) that is a frequency ramp
from 0 to 24000 hz, and then converts it to 44100 hz using a
filter length of 64. It then compares the result to the same
pattern generated at 44100 hz, and outputs the result to the
file "out".
A graph of the correct output should have field 2 and field 4
almost equal (plus/minus 1) up to about sample 40000 (which
corresponds to 20 khz), and then field 2 should be close to 0
above that. Running the test program will print to stdout
something like the following:
time 0.112526
average error 10k=0.4105 22k=639.34
The average error is RMS error over the range [0-10khz] and
[0-22khz], and is expressed in sample values, for an input
amplitude of 16000. Note that RMS errors below 1.0 can't
really be compared, but basically this shows that below
10 khz, the resampler is nearly perfect. Most of the error
is concentrated above 20 khz.
If the average error is significantly larger after modifying
the code, it's probably not good.
dave...
/* Resampling library
* Copyright (C) <2001> David A. Schleef <ds@schleef.org>
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
//#include <ml.h>
#include <resample.h>
#define short_to_double_table
//#define short_to_double_altivec
#define short_to_double_unroll
#ifdef short_to_double_table
static float ints_high[256];
static float ints_low[256];
void conv_double_short_table(double *dest, short *src, int n)
{
static int init = 0;
int i;
unsigned int idx;
if(!init){
for(i=0;i<256;i++){
ints_high[i]=256.0*((i<128)?i:i-256);
ints_low[i]=i;
}
init = 1;
}
if(n&1){
idx = (unsigned short)*src++;
*dest++ = ints_high[(idx>>8)] + ints_low[(idx&0xff)];
n-=1;
}
for(i=0;i<n;i+=2){
idx = (unsigned short)*src++;
*dest++ = ints_high[(idx>>8)] + ints_low[(idx&0xff)];
idx = (unsigned short)*src++;
*dest++ = ints_high[(idx>>8)] + ints_low[(idx&0xff)];
}
}
#endif
#ifdef short_to_double_unroll
void conv_double_short_unroll(double *dest, short *src, int n)
{
if(n&1){
*dest++ = *src++;
n--;
}
if(n&2){
*dest++ = *src++;
*dest++ = *src++;
n-=2;
}
while(n>0){
*dest++ = *src++;
*dest++ = *src++;
*dest++ = *src++;
*dest++ = *src++;
n-=4;
}
}
#endif
void conv_double_short_ref(double *dest, short *src, int n)
{
int i;
for(i=0;i<n;i++){
dest[i]=src[i];
}
}
#ifdef HAVE_CPU_PPC
#if 0
static union { int i[4]; float f[4]; } av_tmp __attribute__ ((__aligned__ (16)));
void conv_double_short_altivec(double *dest, short *src, int n)
{
int i;
for(i=0;i<n;i+=4){
av_tmp.i[0] = src[0];
av_tmp.i[1] = src[1];
av_tmp.i[2] = src[2];
av_tmp.i[3] = src[3];
asm(
" lvx 0,0,%0\n"
" vcfsx 1,0,0\n"
" stvx 1,0,%0\n"
: : "r" (&av_tmp)
);
dest[0]=av_tmp.f[0];
dest[1]=av_tmp.f[1];
dest[2]=av_tmp.f[2];
dest[3]=av_tmp.f[3];
src += 4;
dest += 4;
}
}
#endif
#endif
/* double to short */
void conv_short_double_ref(short *dest, double *src, int n)
{
int i;
double x;
for(i=0;i<n;i++){
x = *src++;
if(x<-32768.0)x=-32768.0;
if(x>32767.0)x=32767.0;
*dest++ = rint(x);
}
}
#ifdef HAVE_CPU_PPC
void conv_short_double_ppcasm(short *dest, double *src, int n)
{
int tmp[2];
double min = -32768.0;
double max = 32767.0;
double ftmp0, ftmp1;
asm __volatile__(
"\taddic. %3,%3,-8\n"
"\taddic. %6,%6,-2\n"
"loop:\n"
"\tlfdu %0,8(%3)\n"
"\tfsub %1,%0,%4\n"
"\tfsel %0,%1,%0,%4\n"
"\tfsub %1,%0,%5\n"
"\tfsel %0,%1,%5,%0\n"
"\tfctiw %1,%0\n"
"\taddic. 5,5,-1\n"
"\tstfd %1,0(%2)\n"
"\tlhz 9,6(%2)\n"
"\tsthu 9,2(%6)\n"
"\tbne loop\n"
: "=&f" (ftmp0), "=&f" (ftmp1)
: "b" (tmp), "r" (src), "f" (min), "f" (max), "r" (dest)
: "r9", "r5" );
}
#endif
void conv_double_short_dstr(double *dest, short *src, int n, int dstr)
{
int i;
void *d = dest;
for(i=0;i<n;i++){
(*(double *)d)=*src++;
d += dstr;
}
}
void conv_short_double_sstr(short *dest, double *src, int n, int sstr)
{
int i;
double x;
void *s = src;
for(i=0;i<n;i++){
x = *(double *)s;
if(x<-32768.0)x=-32768.0;
if(x>32767.0)x=32767.0;
*dest++ = rint(x);
s += sstr;
}
}
/* Resampling library
* Copyright (C) <2001> David A. Schleef <ds@schleef.org>
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <resample.h>
double functable_sinc(void *p,double x)
{
if(x==0)return 1;
return sin(x)/x;
}
double functable_dsinc(void *p,double x)
{
if(x==0)return 0;
return cos(x)/x - sin(x)/(x*x);
}
double functable_window_boxcar(void *p,double x)
{
if(x<-1 || x>1)return 0;
return 1;
}
double functable_window_dboxcar(void *p,double x)
{
return 0;
}
double functable_window_std(void *p,double x)
{
if(x<-1 || x>1)return 0;
return (1-x*x)*(1-x*x);
}
double functable_window_dstd(void *p,double x)
{
if(x<-1 || x>1)return 0;
return -4*x*(1-x*x);
}
void functable_init(functable_t *t)
{
int i;
double x;
t->fx = malloc(sizeof(double)*(t->len+1));
t->fdx = malloc(sizeof(double)*(t->len+1));
t->invoffset = 1.0 / t->offset;
for(i=0;i<t->len+1;i++){
x = t->start + t->offset * i;
x *= t->scale;
t->fx[i] = t->func_x(t->priv,x);
t->fdx[i] = t->scale * t->func_dx(t->priv,x);
}
if(t->func2_x){
double f1x,f1dx;
double f2x,f2dx;
for(i=0;i<t->len+1;i++){
x = t->start + t->offset * i;
x *= t->scale2;
f2x = t->func2_x(t->priv,x);
f2dx = t->scale2 * t->func2_dx(t->priv,x);
f1x = t->fx[i];
f1dx = t->fdx[i];
t->fx[i] = f1x * f2x;
t->fdx[i] = f1x * f2dx + f1dx * f2x;
}
}
}
double functable_eval(functable_t *t,double x)
{
int i;
double f0, f1, w0, w1;
double x2, x3;
double w;
if(x<t->start || x>(t->start+(t->len+1)*t->offset)){
printf("x out of range %g\n",x);
}
x -= t->start;
x /= t->offset;
i = floor(x);
x -= i;
x2 = x * x;
x3 = x2 * x;
f1 = 3 * x2 - 2 * x3;
f0 = 1 - f1;
w0 = (x - 2 * x2 + x3) * t->offset;
w1 = (-x2 + x3) * t->offset;
//printf("i=%d x=%g f0=%g f1=%g w0=%g w1=%g\n",i,x,f0,f1,w0,w1);
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
//w = t->fx[i] * (1-x) + t->fx[i+1] * x;
return w;
}
double functable_fir(functable_t *t, double x, int n, double *data, int len)
{
int i,j;
double f0, f1, w0, w1;
double x2, x3;
double w;
double sum;
x -= t->start;
x /= t->offset;
i = floor(x);
x -= i;
x2 = x * x;
x3 = x2 * x;
f1 = 3 * x2 - 2 * x3;
f0 = 1 - f1;
w0 = (x - 2 * x2 + x3) * t->offset;
w1 = (-x2 + x3) * t->offset;
sum = 0;
for(j=0;j<len;j++){
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
sum += data[j*2] * w;
i += n;
}
return sum;
}
void functable_fir2(functable_t *t, double *r0, double *r1, double x,
int n, double *data, int len)
{
int i,j;
double f0, f1, w0, w1;
double x2, x3;
double w;
double sum0, sum1;
double floor_x;
x -= t->start;
x *= t->invoffset;
floor_x = floor(x);
i = floor_x;
x -= floor_x;
x2 = x * x;
x3 = x2 * x;
f1 = 3 * x2 - 2 * x3;
f0 = 1 - f1;
w0 = (x - 2 * x2 + x3) * t->offset;
w1 = (-x2 + x3) * t->offset;
sum0 = 0;
sum1 = 0;
for(j=0;j<len;j++){
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
sum0 += data[j*2] * w;
sum1 += data[j*2+1] * w;
i += n;
#define unroll2
#define unroll3
#define unroll4
#ifdef unroll2
j++;
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
sum0 += data[j*2] * w;
sum1 += data[j*2+1] * w;
i += n;
#endif
#ifdef unroll3
j++;
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
sum0 += data[j*2] * w;
sum1 += data[j*2+1] * w;
i += n;
#endif
#ifdef unroll4
j++;
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
sum0 += data[j*2] * w;
sum1 += data[j*2+1] * w;
i += n;
#endif
}
*r0 = sum0;
*r1 = sum1;
}
#ifdef unused
void functable_fir2_altivec(functable_t *t, float *r0, float *r1,
double x, int n, float *data, int len)
{
int i,j;
double f0, f1, w0, w1;
double x2, x3;
double w;
double sum0, sum1;
double floor_x;
x -= t->start;
x *= t->invoffset;
floor_x = floor(x);
i = floor_x;
x -= floor_x;
x2 = x * x;
x3 = x2 * x;
f1 = 3 * x2 - 2 * x3;
f0 = 1 - f1;
w0 = (x - 2 * x2 + x3) * t->offset;
w1 = (-x2 + x3) * t->offset;
sum0 = 0;
sum1 = 0;
for(j=0;j<len;j++){
// t->fx, t->fdx needs to be multiplexed by n
// we need 5 consecutive floats, which fit into 2 vecs
// load v0, t->fx[i]
// load v1, t->fx[i+n]
// v2 = v0 (not correct)
// v3 = (v0>>32) || (v1<<3*32) (not correct)
//
// load v4, t->dfx[i]
// load v5, t->dfx[i+n]
// v6 = v4 (not correct)
// v7 = (v4>>32) || (v5<<3*32) (not correct)
//
// v8 = splat(f0)
// v9 = splat(f1)
// v10 = splat(w0)
// v11 = splat(w1)
//
// v12 = v2 * v8
// v12 += v3 * v9
// v12 += v6 * v10
// v12 += v7 * v11
w = t->fx[i] * f0 + t->fx[i + 1] * f1 +
t->fdx[i] * w0 + t->fdx[i + 1] * w1;
// v13 = data[j*2]
// v14 = data[j*2+4]
// v15 = deinterlace_high(v13,v14)
// v16 = deinterlace_low(v13,v14)
// (sum0) v17 += multsum(v13,v15)
// (sum1) v18 += multsum(v14,v16)
sum0 += data[j*2] * w;
sum1 += data[j*2+1] * w;
i += n;
}
*r0 = sum0;
*r1 = sum1;
}
#endif
This diff is collapsed.
/* Resampling library
* Copyright (C) <2001> David Schleef <ds@schleef.org>
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef __RESAMPLE_H__
#define __RESAMPLE_H__
#include <config.h>
typedef struct resample_s resample_t;
struct resample_s {
/* parameters */
int method;
int channels;
int verbose;
int filter_length;
double i_rate;