Commit 3da534ec authored by Monty Montgomery's avatar Monty Montgomery

A second major upgrade to the math behind the Bessel followers used by

the companders.  

a) Soft knees are now computed as a freebie during filter following.
No explicit soft-knee means no sqrt() means considerable optimization

b) streamlined the computation pathway of the followers, partially as
an optimization, partially to prepare for adding the scond tier
over/under filters before release

c) found a bug in the behavior of the new peak-follower math (that was
instigated to repair a more serious bug in the original version);
decay was still using a second order bessel which could still ring
about zero.  Restored the peak-hold behavior to near-original
behavior, but now use the same first-order-decay-limited filters as
the RMS follower.

Depending on platform and compander usage, these changes are good for
getting back 10-25% of the CPU.



git-svn-id: https://svn.xiph.org/trunk/postfish@8033 0101bb08-14d6-0310-b084-bc0e0c8e3800
parent 1436f573
......@@ -30,13 +30,13 @@ SRC = main.c mainpanel.c multibar.c readout.c input.c output.c clippanel.c \
feedback.c freq.c eq.c eqpanel.c compandpanel.c subband.c lpc.c \
bessel.c deverbpanel.c deverb.c singlecomp.c singlepanel.c \
limit.c limitpanel.c mute.c mixpanel.c mix.c reverb.c reverbpanel.c \
outpanel.c config.c window.c
outpanel.c config.c window.c follower.c
OBJ = main.o mainpanel.o multibar.o readout.o input.o output.o clippanel.o \
declip.o reconstruct.o multicompand.o windowbutton.o subpanel.o \
feedback.o freq.o eq.o eqpanel.o compandpanel.o subband.o lpc.o \
bessel.o deverbpanel.o deverb.o singlecomp.o singlepanel.o \
limit.o limitpanel.o mute.o mixpanel.o mix.o reverb.o reverbpanel.o \
outpanel.o config.o window.o
outpanel.o config.o window.o follower.o
GCF = -DETCDIR=\\\"$(ETCDIR)\\\" `pkg-config --cflags gtk+-2.0` -DG_DISABLE_DEPRECATED -DGDK_DISABLE_DEPRECATED -DGTK_DISABLE_DEPRECATED -DGDK_PIXBUF_DISABLE_DEPRECATED
all:
......
......@@ -201,7 +201,7 @@ void compute_iir_symmetric_limited(float *x, int n, iir_state *is,
while(i<n){
double ya= (x[i]+x0*2.+x1)*a_g + y0*a_c0+y1*a_c1;
double yl= y0*l_c0;
if(ya<y0 && ya<yl)ya=yl;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;x[i]=y0=ya;
......@@ -269,7 +269,7 @@ void compute_iir_freefall_limited(float *x, int n, iir_state *is,
x1=x0;x0=x[i];
}
if(ya<y0 && ya<yl)ya=yl;
if(ya<yl)ya=yl;
y1=y0;x[i]=y0=ya;
i++;
......@@ -313,40 +313,225 @@ void compute_iir_freefallonly1(float *x, int n, iir_state *is,
}
/* applies a 1st order freefall filter, followed by a 2nd order attack
filter */
void compute_iir_freefall1_then_symmetric2(float *x, int n, iir_state *is,
iir_filter *attack,
iir_filter *decay){
double a_c0=attack->c[0],d_c0=decay->c[0];
double a_c1=attack->c[1];
double a_g=1./attack->g;
double x0=is->x[0],x1=is->x[1];
double y0=is->y[0],y1=is->y[1];
double yd=is->y[2];
int i=0;
/* a new experiment; bessel followers constructed specifically for
compand functions. Hard and soft knee are implemented as part of
the follower. */
#define prologue double a_c0=attack->c[0],l_c0=limit->c[0]; \
double a_c1=attack->c[1]; \
double a_g=1./(knee * attack->g);\
double x0=is->x[0],x1=is->x[1]; \
double y0=is->y[0],y1=is->y[1]; \
int i=0; \
if(zerome(y0) && zerome(y1)) y0=y1=0.; \
while(i<n)
/* Three delicious filters in one: The 'attack' filter is a fast
second-order Bessel used two ways; as a direct RMS follower and as
a filter pegged to free-fall to the knee value. The 'limit' filter
is a first-order free-fall (to 0) Bessel used to limit the 'attack'
filter to a linear-dB decay. Output is normalized to place the knee
at 1.0 (0dB) */
void compute_iir_over_soft(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj){
double xag=4./attack->g;
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
if(ya<xag)ya=xag;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
adj[i++]-= todB_a((float)ya)*mult;
}
if(zerome(y0) && zerome(y1)) y0=y1=0.;
if(zerome(yd)) yd=0.;
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
while(i<n){
yd*=d_c0;
if(yd<x[i])yd=x[i];
{
double ya= (yd+x0*2.+x1)*a_g + y0*a_c0+y1*a_c1;
x1=x0;x0=yd;
y1=y0;x[i]=y0=ya;
}
void compute_iir_under_soft(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj){
double xag=4./attack->g;
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
if(ya>xag)ya=xag;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
adj[i++]-= todB_a((float)ya)*mult;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
void compute_iir_over_hard(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj){
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
if(ya>1.)adj[i]-= todB_a((float)ya)*mult;
i++;
}
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
is->y[2]=yd;
}
void compute_iir_under_hard(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj){
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
if(ya<1.)adj[i]-= todB_a((float)ya)*mult;
i++;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
/* One more take on the above; we need to be able to gently slew
multiplier changes from one block to the next. Although it seems
absurd to make eight total compander follower variants, doing this
inline avoids a copy later on, and these are enough of a CPU sink
that the silliness is actually worth it. */
void compute_iir_over_soft_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj){
float multadd=(mult2-mult)/n;
double xag=4./attack->g;
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
if(ya<xag)ya=xag;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
adj[i++]-= todB_a((float)ya)*mult;
mult+=multadd;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
void compute_iir_under_soft_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj){
float multadd=(mult2-mult)/n;
double xag=4./attack->g;
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
if(ya>xag)ya=xag;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
adj[i++]-= todB_a((float)ya)*mult;
mult+=multadd;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
void compute_iir_over_hard_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj){
float multadd=(mult2-mult)/n;
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
if(ya>1.)adj[i]-= todB_a((float)ya)*mult;
mult+=multadd;
i++;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
void compute_iir_under_hard_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj){
float multadd=(mult2-mult)/n;
prologue {
double ya = (x[i]+x0*2.+x1)*a_g;
ya += y0*a_c0;
ya += y1*a_c1;
double yl = y0*l_c0;
if(ya<yl)ya=yl;
x1=x0;x0=x[i];
y1=y0;y0=ya;
if(ya<1.)adj[i]-= todB_a((float)ya)*mult;
mult+=multadd;
i++;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
void reset_iir(iir_state *is,float value){
int i;
for(i=0;i<MAXORDER;i++){
is->x[i]=(double)value;
is->y[i]=(double)value;
}
}
......@@ -66,7 +66,7 @@ static inline float impulse_freq4(long ahead){
typedef struct {
double x[MAXORDER];
double y[MAXORDER+1];
double y[MAXORDER];
} iir_state;
extern double mkbessel(double raw_alpha,int order,double *ycoeff);
......@@ -85,7 +85,42 @@ extern void compute_iir_freefall_limited(float *x, int n, iir_state *is,
extern void compute_iir_freefallonly1(float *x, int n, iir_state *is,
iir_filter *decay);
extern void compute_iir_freefall1_then_symmetric2(float *x, int n,
iir_state *is,
iir_filter *attack,
iir_filter *decay);
extern void compute_iir_over_soft(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj);
extern void compute_iir_under_soft(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj);
extern void compute_iir_over_hard(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj);
extern void compute_iir_under_hard(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float *adj);
extern void compute_iir_over_soft_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj);
extern void compute_iir_under_soft_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj);
extern void compute_iir_over_hard_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj);
extern void compute_iir_under_hard_del(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit,
float knee, float mult, float mult2,
float *adj);
extern void reset_iir(iir_state *is,float value);
......@@ -317,8 +317,6 @@ static int determine_average(multi_panel_state *mp,int x){
static void average_change(GtkWidget *w,gpointer in){
multi_panel_state *mp=(multi_panel_state *)in;
if(!mp->updating_av_slider){
multicompand_settings *ms=mp->ms;
banked_multicompand_settings *bc=ms->bc;
cbar *b=mp->bars+multicomp_freqs_max;
cbar *bars=mp->bars;
int bank_active=mp->bank_active;
......@@ -370,8 +368,6 @@ static void slider_change(GtkWidget *w,gpointer in){
int bank_active=mp->bank_active;
int o,u;
int i;
int adj;
u=multibar_get_value(MULTIBAR(b->slider),0);
sprintf(buffer,"%+4ddB",u);
......
......@@ -213,14 +213,14 @@ static void deverb_work_helper(void *vs, deverb_settings *sset){
if(multiplier==sss->prevratio[i]){
for(k=0;k<input_size;k++)
fast[k]=fromdB_a((todB_a(slow+k)-todB_a(fast+k))*.5*multiplier);
fast[k]=fromdB_a((todB_a(slow[k])-todB_a(fast[k]))*.5*multiplier);
}else{
float multiplier_add=(multiplier-sss->prevratio[i])/input_size;
multiplier=sss->prevratio[i];
for(k=0;k<input_size;k++){
fast[k]=fromdB_a((todB_a(slow+k)-todB_a(fast+k))*.5*multiplier);
fast[k]=fromdB_a((todB_a(slow[k])-todB_a(fast[k]))*.5*multiplier);
multiplier+=multiplier_add;
}
......
/*
*
* postfish
*
* Copyright (C) 2002-2004 Monty
*
* Postfish is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2, or (at your option)
* any later version.
*
* Postfish 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Postfish; see the file COPYING. If not, write to the
* Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
*
*
*/
#include "postfish.h"
#include "follower.h"
static void _analysis(char *base,int seq, float *data, int n,int dB,
off_t offset){
FILE *f;
char buf[80];
sprintf(buf,"%s_%d.m",base,seq);
f=fopen(buf,"a");
if(f){
int i;
for(i=0;i<n;i++)
if(dB)
fprintf(f,"%d %f\n",(int)(i+offset),todB(data[i]));
else
fprintf(f,"%d %f\n",(int)(i+offset),(data[i]));
}
fclose(f);
}
off_t offset=0;
int offch;
/* Common follower code */
static void prepare_peak(float *peak, float *x, int n, int ahead,
int hold, peak_state *ps){
int ii,jj;
int loc=ps->loc;
float val=ps->val;
/* Although we have two input_size blocks of zeroes after a
reset, we may still need to look ahead explicitly after a
reset if the lookahead is exceptionally long */
if(loc==0 && val==0){
for(ii=0;ii<ahead;ii++)
if((x[ii]*x[ii])>val){
val=(x[ii]*x[ii]);
loc=ii+hold;
}
}
if(val>peak[0])peak[0]=val;
for(ii=1;ii<n;ii++){
if((x[ii+ahead]*x[ii+ahead])>val){
val=(x[ii+ahead]*x[ii+ahead]);
loc=ii+ahead+hold;
}
if(ii>=loc){
/* backfill */
val=0;
for(jj=ii+ahead-1;jj>=ii;jj--){
if((x[jj]*x[jj])>val)val=(x[jj]*x[jj]);
if(jj<n && val>peak[jj])peak[jj]=val;
}
val=(x[ii+ahead]*x[ii+ahead]);
loc=ii+ahead+hold;
}
if(val>peak[ii])peak[ii]=val;
}
ps->loc=loc-input_size;
ps->val=val;
}
static void fill_work(float *A, float *B, float *work,
int ahead, int hold, int mode, peak_state *ps){
int k;
if(mode){
/* peak mode */
memset(work,0,sizeof(*work)*input_size);
if(B){
float bigcache[input_size*2];
memcpy(bigcache,A,sizeof(*work)*input_size);
memcpy(bigcache+input_size,B,sizeof(*work)*input_size);
prepare_peak(work, bigcache, input_size, ahead, hold, ps);
}else{
prepare_peak(work, A, input_size, ahead, hold, ps);
}
}else{
/* rms mode */
float *cachea=A+ahead;
if(B){
float *worka=work+input_size-ahead;
for(k=0;k<input_size-ahead;k++)
work[k]=cachea[k]*cachea[k];
for(k=0;k<ahead;k++)
worka[k]=B[k]*B[k];
}else{
for(k=0;k<input_size;k++)
work[k]=cachea[k]*cachea[k];
}
}
}
void bi_compand(float *A,float *B,float *adj,
float corner,
float multiplier,
float currmultiplier,
float lookahead,
int mode,int softknee,
iir_filter *attack, iir_filter *decay,
iir_state *iir, peak_state *ps,
int active,
int over){
float work[input_size];
float kneelevel=fromdB(corner*2);
int hold,ahead=(mode?step_ahead(attack->alpha):impulse_ahead2(attack->alpha));
if(ahead>input_size)ahead=input_size;
hold=ahead*(1.-lookahead);
ahead*=lookahead;
fill_work(A,B,work,ahead,hold,mode,ps);
multiplier*=.5;
currmultiplier*=.5;
if(!active || !adj)adj=work;
if(over){
if(multiplier!=currmultiplier){
if(softknee){
compute_iir_over_soft_del(work, input_size, iir, attack, decay,
kneelevel, multiplier, currmultiplier, adj);
}else{
compute_iir_over_hard_del(work, input_size, iir, attack, decay,
kneelevel, multiplier, currmultiplier, adj);
}
}else{
if(softknee){
compute_iir_over_soft(work, input_size, iir, attack, decay,
kneelevel, multiplier, adj);
}else{
compute_iir_over_hard(work, input_size, iir, attack, decay,
kneelevel, multiplier, adj);
}
}
}else{
if(multiplier!=currmultiplier){
if(softknee){
compute_iir_under_soft_del(work, input_size, iir, attack, decay,
kneelevel, multiplier, currmultiplier, adj);
}else{
compute_iir_under_hard_del(work, input_size, iir, attack, decay,
kneelevel, multiplier, currmultiplier, adj);
}
}else{
if(softknee){
compute_iir_under_soft(work, input_size, iir, attack, decay,
kneelevel, multiplier, adj);
}else{
compute_iir_under_hard(work, input_size, iir, attack, decay,
kneelevel, multiplier, adj);
}
}
}
}
void full_compand(float *A,float *B,float *adj,
float multiplier,float currmultiplier,
int mode,
iir_filter *attack, iir_filter *decay,
iir_state *iir, peak_state *ps,
int active){
int k;
float work[input_size];
int ahead=(mode?step_ahead(attack->alpha):impulse_ahead2(attack->alpha));
fill_work(A,B,work,ahead,0,mode,ps);
multiplier*=.5;
currmultiplier*=.5;
compute_iir_symmetric_limited(work, input_size, iir, attack, decay);
if(active){
if(multiplier!=currmultiplier){
float multiplier_add=(currmultiplier-multiplier)/input_size;
for(k=0;k<input_size;k++){
adj[k]-=(todB_a2(work[k])+adj[k])*multiplier;
multiplier+=multiplier_add;
}
}else{
for(k=0;k<input_size;k++)
adj[k]-=(todB_a2(work[k])+adj[k])*multiplier;
}
}
}
/*
*
* postfish
*
* Copyright (C) 2002-2004 Monty
*
* Postfish is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2, or (at your option)
* any later version.
*
* Postfish 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Postfish; see the file COPYING. If not, write to the
* Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
*
*
*/
#include "bessel.h"
typedef struct {
int loc;
float val;
} peak_state;
extern void bi_compand(float *A,float *B,float *adj,
float corner,
float multiplier,
float currmultiplier,
float lookahead,
int mode,int softknee,
iir_filter *attack, iir_filter *decay,
iir_state *iir, peak_state *ps,
int active,
int over);
extern void full_compand(float *A,float *B,float *adj,
float multiplier,float currmultiplier,
int mode,
iir_filter *attack, iir_filter *decay,
iir_state *iir, peak_state *ps,
int active);
......@@ -26,7 +26,10 @@
#include "multicompand.h"
#include <fftw3.h>
#include "subband.h"
#include "bessel.h"
#include "follower.h"
extern off_t offset;
extern int offch;
/* feedback! */
typedef struct multicompand_feedback{
......@@ -37,11 +40,6 @@ typedef struct multicompand_feedback{
int bypass;
} multicompand_feedback;
typedef struct {
int loc;
float val;
} peak_state;
typedef struct {
sig_atomic_t static_u[multicomp_freqs_max];
sig_atomic_t under_ratio;
......@@ -134,9 +132,12 @@ int pull_multicompand_feedback_channel(float **peak,float **rms,int *b){
static void reset_filters_onech(multicompand_state *ms,int ch){
int i;
ms->base_delay[ch]=2;
ms->over_delay[ch]=2;
ms->under_delay[ch]=2;
/* delays are only used to softstart individual filters that went
inactive at unity; the key is that we know they're starting from
mult of zero, which is not necessarily true at a reset */
ms->base_delay[ch]=0;
ms->over_delay[ch]=0;
ms->under_delay[ch]=0;
for(i=0;i<multicomp_freqs_max;i++){
memset(&ms->over_peak[i][ch],0,sizeof(peak_state));
memset(&ms->under_peak[i][ch],0,sizeof(peak_state));
......@@ -249,218 +250,6 @@ static void filterbank_set2(multicompand_state *ms,
}
static void prepare_rms(float *rms, float *xx, int n, int ahead){
int i;
float *x=xx+ahead;
for(i=0;i<n;i++)
rms[i]+=x[i]*x[i];
}
static void prepare_peak(float *peak, float *x, int n, int ahead,
peak_state *ps){
int ii,jj;
int loc=ps->loc;
float val=ps->val;
/* Although we have two input_size blocks of zeroes after a
reset, we may still need to look ahead explicitly after a
reset if the lookahead is exceptionally long */
if(loc==0 && val==0){
for(ii=0;ii<ahead;ii++)
if((x[ii]*x[ii])>val){
val=(x[ii]*x[ii]);
loc=ii;
}
}