From 1861b5d7bdf9fb669633b9ee471f6c57a5633b83 Mon Sep 17 00:00:00 2001 From: James Courtier-Dutton Date: Sat, 15 May 2004 18:22:26 +0000 Subject: Port subwoofer filter from mplayer. CVS patchset: 6543 CVS date: 2004/05/15 18:22:26 --- src/post/audio/Makefile.am | 2 +- src/post/audio/dsp.h | 22 +++ src/post/audio/filter.c | 434 +++++++++++++++++++++++++++++++++++++++++++++ src/post/audio/filter.h | 71 ++++++++ src/post/audio/upmix.c | 138 +++++++------- src/post/audio/window.c | 203 +++++++++++++++++++++ src/post/audio/window.h | 33 ++++ 7 files changed, 832 insertions(+), 71 deletions(-) create mode 100644 src/post/audio/dsp.h create mode 100644 src/post/audio/filter.c create mode 100644 src/post/audio/filter.h create mode 100644 src/post/audio/window.c create mode 100644 src/post/audio/window.h (limited to 'src') diff --git a/src/post/audio/Makefile.am b/src/post/audio/Makefile.am index 9c246ea86..3b2475d19 100644 --- a/src/post/audio/Makefile.am +++ b/src/post/audio/Makefile.am @@ -5,7 +5,7 @@ libdir = $(XINE_PLUGINDIR)/post lib_LTLIBRARIES = xineplug_post_audio_filter_upmix.la xineplug_post_audio_filter_upmix_la_SOURCES = \ - upmix.c + upmix.c filter.c window.c xineplug_post_audio_filter_upmix_la_LIBADD = $(XINE_LIB) xineplug_post_audio_filter_upmix_la_LDFLAGS = -avoid-version -module @XINE_PLUGIN_MIN_SYMS@ -lm diff --git a/src/post/audio/dsp.h b/src/post/audio/dsp.h new file mode 100644 index 000000000..8f5d5eb51 --- /dev/null +++ b/src/post/audio/dsp.h @@ -0,0 +1,22 @@ +/*============================================================================= +// +// This software has been released under the terms of the GNU Public +// license. See http://www.gnu.org/copyleft/gpl.html for details. +// +// Copyright 2002 Anders Johansson ajh@atri.curtin.edu.au +// +//============================================================================= +*/ + +#ifndef _DSP_H +#define _DSP_H 1 + +/* Implementation of routines used for DSP */ + +/* Size of floating point type used in routines */ +#define _ftype_t float + +#include +#include + +#endif diff --git a/src/post/audio/filter.c b/src/post/audio/filter.c new file mode 100644 index 000000000..10f1c975d --- /dev/null +++ b/src/post/audio/filter.c @@ -0,0 +1,434 @@ +/*============================================================================= +// +// This software has been released under the terms of the GNU Public +// license. See http://www.gnu.org/copyleft/gpl.html for details. +// +// Copyright 2001 Anders Johansson ajh@atri.curtin.edu.au +// +//============================================================================= +*/ + +/* Design and implementation of different types of digital filters + +*/ +#include +#include +#include "dsp.h" + +/****************************************************************************** +* FIR filter implementations +******************************************************************************/ + +/* C implementation of FIR filter y=w*x + + n number of filter taps, where mod(n,4)==0 + w filter taps + x input signal must be a circular buffer which is indexed backwards +*/ +inline _ftype_t fir(register unsigned int n, _ftype_t* w, _ftype_t* x) +{ + register _ftype_t y; // Output + y = 0.0; + do{ + n--; + y+=w[n]*x[n]; + }while(n != 0); + return y; +} + +/* C implementation of parallel FIR filter y(k)=w(k) * x(k) (where * denotes convolution) + + n number of filter taps, where mod(n,4)==0 + d number of filters + xi current index in xq + w filter taps k by n big + x input signal must be a circular buffers which are indexed backwards + y output buffer + s output buffer stride +*/ +inline _ftype_t* pfir(unsigned int n, unsigned int d, unsigned int xi, _ftype_t** w, _ftype_t** x, _ftype_t* y, unsigned int s) +{ + register _ftype_t* xt = *x + xi; + register _ftype_t* wt = *w; + register int nt = 2*n; + while(d-- > 0){ + *y = fir(n,wt,xt); + wt+=n; + xt+=nt; + y+=s; + } + return y; +} + +/* Add new data to circular queue designed to be used with a parallel + FIR filter, with d filters. xq is the circular queue, in pointing + at the new samples, xi current index in xq and n the length of the + filter. xq must be n*2 by k big, s is the index for in. +*/ +inline int updatepq(unsigned int n, unsigned int d, unsigned int xi, _ftype_t** xq, _ftype_t* in, unsigned int s) +{ + register _ftype_t* txq = *xq + xi; + register int nt = n*2; + + while(d-- >0){ + *txq= *(txq+n) = *in; + txq+=nt; + in+=s; + } + return (++xi)&(n-1); +} + +/****************************************************************************** +* FIR filter design +******************************************************************************/ + +/* Design FIR filter using the Window method + + n filter length must be odd for HP and BS filters + w buffer for the filter taps (must be n long) + fc cutoff frequencies (1 for LP and HP, 2 for BP and BS) + 0 < fc < 1 where 1 <=> Fs/2 + flags window and filter type as defined in filter.h + variables are ored together: i.e. LP|HAMMING will give a + low pass filter designed using a hamming window + opt beta constant used only when designing using kaiser windows + + returns 0 if OK, -1 if fail +*/ +int design_fir(unsigned int n, _ftype_t* w, _ftype_t* fc, unsigned int flags, _ftype_t opt) +{ + unsigned int o = n & 1; // Indicator for odd filter length + unsigned int end = ((n + 1) >> 1) - o; // Loop end + unsigned int i; // Loop index + + _ftype_t k1 = 2 * M_PI; // 2*pi*fc1 + _ftype_t k2 = 0.5 * (_ftype_t)(1 - o);// Constant used if the filter has even length + _ftype_t k3; // 2*pi*fc2 Constant used in BP and BS design + _ftype_t g = 0.0; // Gain + _ftype_t t1,t2,t3; // Temporary variables + _ftype_t fc1,fc2; // Cutoff frequencies + + // Sanity check + if(!w || (n == 0)) return -1; + + // Get window coefficients + switch(flags & WINDOW_MASK){ + case(BOXCAR): + boxcar(n,w); break; + case(TRIANG): + triang(n,w); break; + case(HAMMING): + hamming(n,w); break; + case(HANNING): + hanning(n,w); break; + case(BLACKMAN): + blackman(n,w); break; + case(FLATTOP): + flattop(n,w); break; + case(KAISER): + kaiser(n,w,opt); break; + default: + return -1; + } + + if(flags & (LP | HP)){ + fc1=*fc; + // Cutoff frequency must be < 0.5 where 0.5 <=> Fs/2 + fc1 = ((fc1 <= 1.0) && (fc1 > 0.0)) ? fc1/2 : 0.25; + k1 *= fc1; + + if(flags & LP){ // Low pass filter + + // If the filter length is odd, there is one point which is exactly + // in the middle. The value at this point is 2*fCutoff*sin(x)/x, + // where x is zero. To make sure nothing strange happens, we set this + // value separately. + if (o){ + w[end] = fc1 * w[end] * 2.0; + g=w[end]; + } + + // Create filter + for (i=0 ; i Fs/2 + fc1 = ((fc1 <= 1.0) && (fc1 > 0.0)) ? fc1/2 : 0.25; + fc2 = ((fc2 <= 1.0) && (fc2 > 0.0)) ? fc2/2 : 0.25; + k3 = k1 * fc2; // 2*pi*fc2 + k1 *= fc1; // 2*pi*fc1 + + if(flags & BP){ // Band pass + // Calculate center tap + if (o){ + g=w[end]*(fc1+fc2); + w[end] = (fc2 - fc1) * w[end] * 2.0; + } + + // Create filter + for (i=0 ; i HP filter + + returns 0 if OK, -1 if fail +*/ +int design_pfir(unsigned int n, unsigned int k, _ftype_t* w, _ftype_t** pw, _ftype_t g, unsigned int flags) +{ + int l = (int)n/k; // Length of individual FIR filters + int i; // Counters + int j; + _ftype_t t; // g * w[i] + + // Sanity check + if(l<1 || k<1 || !w || !pw) + return -1; + + // Do the stuff + if(flags&REW){ + for(j=l-1;j>-1;j--){//Columns + for(i=0;i<(int)k;i++){//Rows + t=g * *w++; + pw[i][j]=t * ((flags & ODD) ? ((j & 1) ? -1 : 1) : 1); + } + } + } + else{ + for(j=0;j1000.0 || Q< 1.0)) + return -1; + + memcpy(at,a,3*sizeof(_ftype_t)); + memcpy(bt,b,3*sizeof(_ftype_t)); + + bt[1]/=Q; + + /* Calculate a and b and overwrite the original values */ + prewarp(at, fc, fs); + prewarp(bt, fc, fs); + /* Execute bilinear transform */ + bilinear(at, bt, k, fs, coef); + + return 0; +} + diff --git a/src/post/audio/filter.h b/src/post/audio/filter.h new file mode 100644 index 000000000..b64eeedbd --- /dev/null +++ b/src/post/audio/filter.h @@ -0,0 +1,71 @@ +/*============================================================================= +// +// This software has been released under the terms of the GNU Public +// license. See http://www.gnu.org/copyleft/gpl.html for details. +// +// Copyright 2001 Anders Johansson ajh@atri.curtin.edu.au +// +//============================================================================= +*/ + +#if !defined _DSP_H +# error "Never use directly; include instead" +#endif + +#ifndef _FILTER_H +#define _FILTER_H 1 + + +// Design and implementation of different types of digital filters + + +// Flags used for filter design + +// Filter characteristics +#define LP 0x00010000 // Low pass +#define HP 0x00020000 // High pass +#define BP 0x00040000 // Band pass +#define BS 0x00080000 // Band stop +#define TYPE_MASK 0x000F0000 + +// Window types +#define BOXCAR 0x00000001 +#define TRIANG 0x00000002 +#define HAMMING 0x00000004 +#define HANNING 0x00000008 +#define BLACKMAN 0x00000010 +#define FLATTOP 0x00000011 +#define KAISER 0x00000012 +#define WINDOW_MASK 0x0000001F + +// Parallel filter design +#define FWD 0x00000001 // Forward indexing of polyphase filter +#define REW 0x00000002 // Reverse indexing of polyphase filter +#define ODD 0x00000010 // Make filter HP + +// Exported functions +extern _ftype_t fir(unsigned int n, _ftype_t* w, _ftype_t* x); + +extern _ftype_t* pfir(unsigned int n, unsigned int k, unsigned int xi, _ftype_t** w, _ftype_t** x, _ftype_t* y, unsigned int s); + +extern int updateq(unsigned int n, unsigned int xi, _ftype_t* xq, _ftype_t* in); +extern int updatepq(unsigned int n, unsigned int k, unsigned int xi, _ftype_t** xq, _ftype_t* in, unsigned int s); + +extern int design_fir(unsigned int n, _ftype_t* w, _ftype_t* fc, unsigned int flags, _ftype_t opt); + +extern int design_pfir(unsigned int n, unsigned int k, _ftype_t* w, _ftype_t** pw, _ftype_t g, unsigned int flags); +extern void prewarp(_ftype_t* a, _ftype_t fc, _ftype_t fs); +void bilinear(_ftype_t* a, _ftype_t* b, _ftype_t* k, _ftype_t fs, _ftype_t *coef); + +extern int szxform(_ftype_t* a, _ftype_t* b, _ftype_t Q, _ftype_t fc, _ftype_t fs, _ftype_t *k, _ftype_t *coef); + +/* Add new data to circular queue designed to be used with a FIR + filter. xq is the circular queue, in pointing at the new sample, xi + current index for xq and n the length of the filter. xq must be n*2 + long. +*/ +#define updateq(n,xi,xq,in)\ + xq[xi]=(xq)[(xi)+(n)]=*(in);\ + xi=(++(xi))&((n)-1); + +#endif diff --git a/src/post/audio/upmix.c b/src/post/audio/upmix.c index 777f552de..58c0ec55e 100644 --- a/src/post/audio/upmix.c +++ b/src/post/audio/upmix.c @@ -23,7 +23,7 @@ * process. It simply paints the screen a solid color and rotates through * colors on each iteration. * - * $Id: upmix.c,v 1.1 2004/05/15 15:32:47 jcdutton Exp $ + * $Id: upmix.c,v 1.2 2004/05/15 18:22:26 jcdutton Exp $ * */ @@ -32,6 +32,7 @@ #include "xine_internal.h" #include "xineutils.h" #include "post.h" +#include "dsp.h" #define FPS 20 @@ -49,6 +50,39 @@ struct post_class_upmix_s { xine_t *xine; }; + +/* Q value for low-pass filter */ +#define Q 1.0 + +/* Analog domain biquad section */ +typedef struct{ + float a[3]; // Numerator coefficients + float b[3]; // Denominator coefficients +} biquad_t; + +/* S-parameters for designing 4th order Butterworth filter */ +static biquad_t s_param[2] = {{{1.0,0.0,0.0},{1.0,0.765367,1.0}}, + {{1.0,0.0,0.0},{1.0,1.847759,1.0}}}; + +/* Data for specific instances of this filter */ +typedef struct af_sub_s +{ + float w[2][4]; /* Filter taps for low-pass filter */ + float q[2][2]; /* Circular queues */ + float fc; /* Cutoff frequency [Hz] for low-pass filter */ + float k; /* Filter gain */ +}af_sub_t; + +#ifndef IIR +#define IIR(in,w,q,out) { \ + float h0 = (q)[0]; \ + float h1 = (q)[1]; \ + float hn = (in) - h0 * (w)[0] - h1 * (w)[1]; \ + out = hn + h0 * (w)[2] + h1 * (w)[3]; \ + (q)[1] = h0; \ + (q)[0] = hn; \ +} +#endif struct post_plugin_upmix_s { post_plugin_t post; @@ -65,6 +99,7 @@ struct post_plugin_upmix_s { int data_idx; short data [2][NUMSAMPLES]; audio_buffer_t *buf; /* dummy buffer just to hold a copy of audio data */ + af_sub_t *sub; int channels; int channels_out; @@ -129,6 +164,22 @@ static int upmix_port_open(xine_audio_port_t *port_gen, xine_stream_t *stream, } else { this->channels_out=2; } + + this->sub = xine_xmalloc(sizeof(af_sub_t)); + if (!this->sub) { + return 0; + } + this->sub->fc = 60; /* LFE Cutoff frequency 60Hz */ + this->sub->k = 1.0; + if((-1 == szxform(s_param[0].a, s_param[0].b, Q, this->sub->fc, + (float)rate, &this->sub->k, this->sub->w[0])) || + (-1 == szxform(s_param[1].a, s_param[1].b, Q, this->sub->fc, + (float)rate, &this->sub->k, this->sub->w[1]))) { + free(this->sub); + this->sub=NULL; + return 0; + } + this->samples_per_frame = rate / FPS; this->data_idx = 0; @@ -153,7 +204,7 @@ static void upmix_port_close(xine_audio_port_t *port_gen, xine_stream_t *stream _x_post_dec_usage(port); } -int upmix_frames_2to51_16bit(uint8_t *dst8, uint8_t *src8, int num_frames) { +static int upmix_frames_2to51_16bit(uint8_t *dst8, uint8_t *src8, int num_frames, af_sub_t *sub) { int16_t *dst=(int16_t *)dst8; int16_t *src=(int16_t *)src8; @@ -163,16 +214,26 @@ int upmix_frames_2to51_16bit(uint8_t *dst8, uint8_t *src8, int num_frames) { int dst_num_channels=6; int src_frame; int dst_frame; + float sample; + int32_t sum; for (frame=0;frame < num_frames; frame++) { dst_frame=frame*dst_num_channels*bytes_per_sample; src_frame=frame*src_num_channels*bytes_per_sample; dst[dst_frame] = src[src_frame]; dst[dst_frame+(1*bytes_per_sample)] = src[src_frame+(1*bytes_per_sample)]; - dst[dst_frame+(2*bytes_per_sample)] = (src[src_frame] - src[src_frame+(1*bytes_per_sample)]) / 2; /* try a bit of dolby */ + /* try a bit of dolby */ + /* FIXME: Dobly surround is a bit more complicated than this, but this is a start. */ + dst[dst_frame+(2*bytes_per_sample)] = (src[src_frame] - src[src_frame+(1*bytes_per_sample)]) / 2; dst[dst_frame+(3*bytes_per_sample)] = (src[src_frame] - src[src_frame+(1*bytes_per_sample)]) / 2; - dst[dst_frame+(4*bytes_per_sample)] = (src[src_frame] + src[src_frame+(1*bytes_per_sample)]) / 2; - dst[dst_frame+(5*bytes_per_sample)] = 0; + sum = ((int32_t)src[src_frame] + (int32_t)src[src_frame+(1*bytes_per_sample)]) / 2; + dst[dst_frame+(4*bytes_per_sample)] = (int16_t)sum; + /* Create the LFE channel using a low pass filter */ + /* filter feature ported from mplayer */ + sample = (1.0/SHRT_MAX)*((float)sum); + IIR(sample * sub->k, sub->w[0], sub->q[0], sample); + IIR(sample , sub->w[1], sub->q[1], sample); + dst[dst_frame+(5*bytes_per_sample)] = (int16_t)(sample * SHRT_MAX); } return frame; } @@ -213,8 +274,6 @@ static void upmix_port_put_buffer (xine_audio_port_t *port_gen, this->buf->format.rate = port->rate; this->buf->format.mode = AO_CAP_MODE_5_1CHANNEL; _x_extra_info_merge( this->buf->extra_info, buf->extra_info); - // xine_stream_t *stream; /* stream that send that buffer */ - /* FIXME: This does 2 to 5.1 channel upmix */ step_channel = this->buf->format.bits>>3; dst_step_frame = this->channels_out*step_channel; src_step_frame = this->channels*step_channel; @@ -226,78 +285,18 @@ static void upmix_port_put_buffer (xine_audio_port_t *port_gen, data8src=(int8_t*)buf->mem; data8src+=num_frames_processed*src_step_frame; data8dst=(int8_t*)this->buf->mem; - num_frames_done = upmix_frames_2to51_16bit(data8dst, data8src, num_frames); + num_frames_done = upmix_frames_2to51_16bit(data8dst, data8src, num_frames, this->sub); this->buf->num_frames = num_frames_done; num_frames_processed+= num_frames_done; /* pass data to original port */ port->original_port->put_buffer(port->original_port, this->buf, stream ); } } - //printf("num_frames_done=%d, num_frames=%d\n",num_frames_done, num_frames); /* free data from origial buffer */ buf->num_frames=0; /* UNDOCUMENTED, but hey, it works! Force old audio_out buffer free. */ port->original_port->put_buffer(port->original_port, buf, stream ); - /* pass data to original port */ - /* we must not use original data anymore, it should have already being moved - * to the fifo of free audio buffers. just use our private copy instead. - */ - return; - - buf = this->buf; - - this->sample_counter += buf->num_frames; - - j = (this->channels >= 2) ? 1 : 0; - - do { - - if( port->bits == 8 ) { - data8 = (int8_t *)buf->mem; - data8 += samples_used * this->channels; - - /* scale 8 bit data to 16 bits and convert to signed as well */ - for( i = 0; i < buf->num_frames && this->data_idx < NUMSAMPLES; - i++, this->data_idx++, data8 += this->channels ) { - this->data[0][this->data_idx] = ((int16_t)data8[0] << 8) - 0x8000; - this->data[1][this->data_idx] = ((int16_t)data8[j] << 8) - 0x8000; - } - } else { - data = buf->mem; - data += samples_used * this->channels; - - for( i = 0; i < buf->num_frames && this->data_idx < NUMSAMPLES; - i++, this->data_idx++, data += this->channels ) { - this->data[0][this->data_idx] = data[0]; - this->data[1][this->data_idx] = data[j]; - } - } - - if( this->sample_counter >= this->samples_per_frame && - this->data_idx == NUMSAMPLES ) { - - this->data_idx = 0; - samples_used += this->samples_per_frame; - - frame = this->vo_port->get_frame (this->vo_port, FOO_WIDTH, FOO_HEIGHT, - this->ratio, XINE_IMGFMT_YUY2, - VO_BOTH_FIELDS); - frame->extra_info->invalid = 1; - frame->bad_frame = 0; - frame->duration = 90000 * this->samples_per_frame / port->rate; - frame->pts = pts; - this->metronom->got_video_frame(this->metronom, frame); - - this->sample_counter -= this->samples_per_frame; - - memset(frame->base[0], this->current_yuv_byte, FOO_WIDTH * FOO_HEIGHT * 2); - this->current_yuv_byte += 3; - - frame->draw(frame, NULL); - frame->free(frame); - } - } while( this->sample_counter >= this->samples_per_frame ); } static void upmix_dispose(post_plugin_t *this_gen) @@ -305,9 +304,8 @@ static void upmix_dispose(post_plugin_t *this_gen) post_plugin_upmix_t *this = (post_plugin_upmix_t *)this_gen; if (_x_post_dispose(this_gen)) { - this->metronom->exit(this->metronom); - + if (this->sub) free(this->sub); free(this); } } diff --git a/src/post/audio/window.c b/src/post/audio/window.c new file mode 100644 index 000000000..25b92bc31 --- /dev/null +++ b/src/post/audio/window.c @@ -0,0 +1,203 @@ +/*============================================================================= +// +// This software has been released under the terms of the GNU Public +// license. See http://www.gnu.org/copyleft/gpl.html for details. +// +// Copyright 2001 Anders Johansson ajh@atri.curtin.edu.au +// +//============================================================================= +*/ + +/* Calculates a number of window functions. The following window + functions are currently implemented: Boxcar, Triang, Hanning, + Hamming, Blackman, Flattop and Kaiser. In the function call n is + the number of filter taps and w the buffer in which the filter + coefficients will be stored. +*/ + +#include +#include "dsp.h" + +/* +// Boxcar +// +// n window length +// w buffer for the window parameters +*/ +void boxcar(int n, _ftype_t* w) +{ + int i; + // Calculate window coefficients + for (i=0 ; i> 1; + int i; + + // Calculate window coefficients + for (i=0 ; i= BIZ_EPSILON * sum); + return(sum); +} + +/* +// Kaiser +// +// n window length +// w buffer for the window parameters +// b beta parameter of Kaiser window, Beta >= 1 +// +// Beta trades the rejection of the low pass filter against the +// transition width from passband to stop band. Larger Beta means a +// slower transition and greater stop band rejection. See Rabiner and +// Gold (Theory and Application of DSP) under Kaiser windows for more +// about Beta. The following table from Rabiner and Gold gives some +// feel for the effect of Beta: +// +// All ripples in dB, width of transition band = D*N where N = window +// length +// +// BETA D PB RIP SB RIP +// 2.120 1.50 +-0.27 -30 +// 3.384 2.23 0.0864 -40 +// 4.538 2.93 0.0274 -50 +// 5.658 3.62 0.00868 -60 +// 6.764 4.32 0.00275 -70 +// 7.865 5.0 0.000868 -80 +// 8.960 5.7 0.000275 -90 +// 10.056 6.4 0.000087 -100 +*/ +void kaiser(int n, _ftype_t* w, _ftype_t b) +{ + _ftype_t tmp; + _ftype_t k1 = 1.0/besselizero(b); + int k2 = 1 - (n & 1); + int end = (n + 1) >> 1; + int i; + + // Calculate window coefficients + for (i=0 ; i directly; include instead" +#endif + +#ifndef _WINDOW_H +#define _WINDOW_H 1 + +extern void boxcar(int n, _ftype_t* w); +extern void triang(int n, _ftype_t* w); +extern void hanning(int n, _ftype_t* w); +extern void hamming(int n,_ftype_t* w); +extern void blackman(int n,_ftype_t* w); +extern void flattop(int n,_ftype_t* w); +extern void kaiser(int n, _ftype_t* w,_ftype_t b); + +#endif -- cgit v1.2.3