Commit 0b49d510 authored by Monty Montgomery's avatar Monty Montgomery

Farily large bugfixes to the Bessel followers and peak tracking



git-svn-id: https://svn.xiph.org/trunk/postfish@7913 0101bb08-14d6-0310-b084-bc0e0c8e3800
parent 51022785
......@@ -183,138 +183,104 @@ double mkbessel(double raw_alpha,int order,double *ycoeff){
return hypot(dc_gain.re,dc_gain.im);
}
/* assymetrical attack/decay filter computation */
/* this one is designed for fast attack, slow decay */
void compute_iir_fast_attack2(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],d_c1=decay->c[1];
double a_g=attack->g, d_g=decay->g;
/* applies a 2nd order filter (attack) that is decay-limited by a
first-order freefall filter (decay) */
void compute_iir_symmetric_limited(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit){
double a_c0=attack->c[0],l_c0=limit->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];
int state=is->state;
int i=0;
if(zerome(y0) && zerome(y1)){
y0=y1=0.;
if(zerome(y0) && zerome(y1)) y0=y1=0.;
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;
x1=x0;x0=x[i];
y1=y0;x[i]=y0=ya;
i++;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
}
/* applies a 2nd order filter (decay) to decay from peak value only,
decay limited by a first-order freefall filter (limit) with the
same alpha as decay */
void compute_iir_decay_limited(float *x, int n, iir_state *is,
iir_filter *decay, iir_filter *limit){
double d_c0=decay->c[0],l_c0=limit->c[0];
double d_c1=decay->c[1];
double d_g=1./decay->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.;
if(x[0]>y0)state=0;
while(i<n){
double yd= (x[i]+x0*2.+x1)*d_g + y0*d_c0+y1*d_c1;
double yl= y0*l_c0;
if(yd<yl)yd=yl;
if(yd<x[i])y1=y0=yd=x[i];
if(state==0){
/* attack case */
while(i<n){
double ya= (x[i]+x0*2.+x1)/a_g + y0*a_c0+y1*a_c1;
if(x[i]<y0 && ya<y0){
state=1;
break;
}
x1=x0;x0=x[i];
y1=y0;x[i]=y0=ya;
i++;
}
}
if(state==1){
/* decay case */
if(y1<y0){
/* decay fixup needed because we're in discontinuous time */
y1=y0;
}
while(1){
double yd = (x[i]+x0*2.+x1)/d_g + y0*d_c0+y1*d_c1;
x1=x0;x0=x[i];
y1=y0;x[i]=y0=yd;
i++;
if(i>=n)break;
if(x[i]>y0){
state=0;
break;
}
}
}
x1=x0;x0=x[i];
y1=y0;x[i]=y0=yd;
i++;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
is->state=state;
}
/* this one is designed for fast decay, slow attack */
void compute_iir_fast_decay2(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],d_c1=decay->c[1];
double a_g=attack->g, d_g=decay->g;
/* applies a 2nd order filter (attack) that is decay-limited by a
first-order filter (decay) */
void compute_iir_freefall_limited(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit){
double a_c0=attack->c[0],l_c0=limit->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];
int state=is->state;
int i=0;
if(zerome(y0) && zerome(y1)){
y0=y1=0.;
}
if(zerome(y0) && zerome(y1)) y0=y1=0.;
if(x[0]<y0)state=1;
while(i<n){
if(state==1){
/* decay case */
while(i<n){
double yd= (x[i]+x0*2.+x1)/d_g + y0*d_c0+y1*d_c1;
if(x[i]>y0 && yd>y0){
state=0;
break;
}
x1=x0;x0=x[i];
y1=y0;x[i]=y0=yd;
i++;
}
}
double ya= (x0*2.+x1)*a_g + y0*a_c0+y1*a_c1;
double yl= y0*l_c0;
if(state==0){
/* attack case */
if(y1>y0){
/* attack fixup needed because we're in discontinuous time */
y1=y0;
}
while(1){
double ya = (x[i]+x0*2.+x1)/a_g + y0*a_c0+y1*a_c1;
x1=x0;x0=x[i];
y1=y0;x[i]=y0=ya;
i++;
if(i>=n)break;
if(x[i]<y0){
state=1;
break;
}
}
if(x[i]<ya){
x1=x0;x0=0;
}else{
ya+= x[i]*a_g;
x1=x0;x0=x[i];
}
if(ya<y0 && ya<yl)ya=yl;
y1=y0;x[i]=y0=ya;
i++;
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
is->state=state;
}
/* allow decay to proceed in freefall */
void compute_iir_only_freefall1(float *x, int n, iir_state *is,
iir_filter *decay){
void compute_iir_freefallonly1(float *x, int n, iir_state *is,
iir_filter *decay){
double d_c0=decay->c[0];
double x0=is->x[0];
......@@ -347,123 +313,40 @@ void compute_iir_only_freefall1(float *x, int n, iir_state *is,
}
void compute_iir_decayonly2(float *x, int n, iir_state *is,
iir_filter *decay){
double d_c0=decay->c[0];
double d_c1=decay->c[1];
double d_g=decay->g;
/* 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];
double x0=is->x[0];
double x1=is->x[1];
double y0=is->y[0];
double y1=is->y[1];
int i=0;
if(zerome(y0) && zerome(y1)){
y0=y1=0.;
}
if(zerome(y0) && zerome(y1)) y0=y1=0.;
if(zerome(yd)) yd=0.;
while(i<n){
double yd;
if(y1<y0)y1=y0; // slope fixup
yd = (x[i]+x0*2.+x1)/d_g + y0*d_c0+y1*d_c1;
if(x[i]>yd)yd=x[i];
yd*=d_c0;
if(yd<x[i])yd=x[i];
x1=x0;x0=x[i];
y1=y0;x[i]=y0=yd;
i++;
}
is->x[0]=x0;
is->x[1]=x1;
is->y[0]=y0;
is->y[1]=y1;
}
/* symmetric filter computation */
void compute_iir_symmetric2(float *x, int n, iir_state *is,
iir_filter *filter){
double c0=filter->c[0];
double c1=filter->c[1];
double g=filter->g;
double x0=is->x[0];
double x1=is->x[1];
double y0=is->y[0];
double y1=is->y[1];
int i=0;
if(zerome(y0) && zerome(y1)){
y0=y1=0.;
}
while(i<n){
double yd= (x[i]+x0*2.+x1)/g + y0*c0+y1*c1;
x1=x0;x0=x[i];
y1=y0;x[i]=y0=yd;
i++;
}
is->x[0]=x0;
is->x[1]=x1;
is->y[0]=y0;
is->y[1]=y1;
}
void compute_iir_symmetric_freefall2(float *x, int n, iir_state *is,
iir_filter *filter){
double c0=filter->c[0];
double c1=filter->c[1];
double g=filter->g;
double x0=is->x[0];
double x1=is->x[1];
double y0=is->y[0];
double y1=is->y[1];
int i=0;
{
double ya= (yd+x0*2.+x1)*a_g + y0*a_c0+y1*a_c1;
if(zerome(y0) && zerome(y1)){
y0=y1=0.;
}
while(i<n){
double yd= (x0*2.+x1)/g + y0*c0+y1*c1;
if(x[i]<yd){
x1=x0;x0=0;
}else{
yd+= x[i]/g;
x1=x0;x0=x[i];
x1=x0;x0=yd;
y1=y0;x[i]=y0=ya;
}
y1=y0;x[i]=y0=yd;
i++;
}
is->x[0]=x0;
is->x[1]=x1;
is->y[0]=y0;
is->y[1]=y1;
}
/* filter decision wrapper */
void compute_iir2(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *decay){
if (attack->alpha > decay->alpha){
/* fast attack, slow decay */
compute_iir_fast_attack2(x, n, is, attack, decay);
}else{
compute_iir_symmetric2(x, n, is, attack);
}
is->x[0]=x0;is->x[1]=x1;
is->y[0]=y0;is->y[1]=y1;
is->y[2]=yd;
}
......@@ -24,7 +24,7 @@
#include "postfish.h"
extern int input_rate;
#define MAXORDER 4
#define MAXORDER 2
typedef struct {
double c[MAXORDER];
......@@ -66,25 +66,26 @@ static inline float impulse_freq4(long ahead){
typedef struct {
double x[MAXORDER];
double y[MAXORDER];
int state;
double y[MAXORDER+1];
} iir_state;
extern double mkbessel(double raw_alpha,int order,double *ycoeff);
extern void compute_iir_fast_attack2(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *decay);
extern void compute_iir_fast_decay2(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *decay);
extern void compute_iir_symmetric2(float *x, int n, iir_state *is,
iir_filter *filter);
extern void compute_iir2(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *decay);
extern void compute_iir_only_freefall1(float *x, int n, iir_state *is,
iir_filter *decay);
extern void compute_iir_decayonly2(float *x, int n, iir_state *is,
iir_filter *decay);
extern void compute_iir_symmetric_freefall2(float *x, int n, iir_state *is,
iir_filter *filter);
extern void compute_iir_symmetric_limited(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit);
extern void compute_iir_decay_limited(float *x, int n, iir_state *is,
iir_filter *decay, iir_filter *limit);
extern void compute_iir_freefall_limited(float *x, int n, iir_state *is,
iir_filter *attack, iir_filter *limit);
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);
......@@ -49,6 +49,7 @@ typedef struct {
subband_state ss;
iir_filter smooth;
iir_filter smoothlimit;
iir_filter release;
iir_state *iirS[deverb_freqs];
......@@ -113,28 +114,6 @@ int deverb_load(void){
return 0;
}
#if 0
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);
}
#endif
static void deverb_work_helper(void *vs, deverb_settings *sset){
deverb_state *sss=(deverb_state *)vs;
subband_state *ss=&sss->ss;
......@@ -142,10 +121,14 @@ static void deverb_work_helper(void *vs, deverb_settings *sset){
float smoothms=sset->smooth*.1;
float releasems=sset->release*.1;
iir_filter *smooth=&sss->smooth;
iir_filter *smoothlimit=&sss->smoothlimit;
iir_filter *release=&sss->release;
int ahead;
if(smoothms!=smooth->ms)filter_set(ss,smoothms,smooth,1,2);
if(smoothms!=smooth->ms){
filter_set(ss,smoothms,smooth,1,2);
filter_set(ss,smoothms,smoothlimit,0,1);
}
if(releasems!=release->ms)filter_set(ss,releasems,release,0,1);
ahead=impulse_ahead2(smooth->alpha);
......@@ -186,17 +169,14 @@ static void deverb_work_helper(void *vs, deverb_settings *sset){
if(sset->linkp==0 || firstlink==1){
//memcpy(slow,fast,sizeof(slow));
compute_iir_freefall_limited(fast, input_size, &sss->iirS[i][j],
smooth,smoothlimit);
compute_iir_symmetric_freefall2(fast, input_size, &sss->iirS[i][j],
smooth);
memcpy(slow,fast,sizeof(slow));
compute_iir_only_freefall1(slow, input_size, &sss->iirR[i][j],
release);
compute_iir_freefallonly1(slow, input_size, &sss->iirR[i][j],
release);
//_analysis("fast",i,fast,input_size,1,offset);
//_analysis("slow",i,slow,input_size,1,offset);
//_analysis("fast3",i,fast,input_size,1,offset);
if(multiplier==sss->prevratio[i]){
......@@ -214,7 +194,7 @@ static void deverb_work_helper(void *vs, deverb_settings *sset){
}
//_analysis("adj",i,fast,input_size,1,offset);
//_analysis("adj3",i,fast,input_size,1,offset);
if(sset->linkp && firstlink==1){
......@@ -245,6 +225,7 @@ static void deverb_work_helper(void *vs, deverb_settings *sset){
sss->prevratio[i]=multiplier;
}
//offset+=input_size;
}
static void deverb_work_channel(void *vs){
......
......@@ -219,8 +219,8 @@ static deverb_panel_state *deverbpanel_create_helper(postfish_mainpanel *mp,
multibar_callback(MULTIBAR(slider),timing_change,&ps->timing);
multibar_thumb_set(MULTIBAR(slider),0,100);
multibar_thumb_set(MULTIBAR(slider),1,400);
multibar_thumb_set(MULTIBAR(slider),100,0);
multibar_thumb_set(MULTIBAR(slider),400,1);
gtk_table_attach(GTK_TABLE(table),slider,1,2,1,2,
GTK_FILL|GTK_EXPAND,GTK_EXPAND,5,0);
......
......@@ -27,6 +27,7 @@
#include "feedback.h"
#include "freq.h"
#include "lpc.h"
#include "window.h"
extern int input_rate;
extern int input_size;
......@@ -76,11 +77,7 @@ int freq_class_load(freq_class_setup *f,const float *frequencies, int bands){
f->bands=bands;
/* fill in time window */
f->window=malloc(blocksize*sizeof(*f->window));
for(i=0;i<blocksize;i++)f->window[i]=sin(M_PIl*i/blocksize);
for(i=0;i<blocksize;i++)f->window[i]*=f->window[i];
for(i=0;i<blocksize;i++)f->window[i]=sin(f->window[i]*M_PIl*.5);
for(i=0;i<blocksize;i++)f->window[i]*=f->window[i];
f->window=window_get(3,f->qblocksize);
f->fftwf_buffer = fftwf_malloc(sizeof(*f->fftwf_buffer) *
(f->qblocksize*4+2));
......@@ -290,11 +287,9 @@ static void fill_freq_buffer_helper(float *buffer,float *window,
memset(buffer+qblocksize*3,0,sizeof(*buffer)*qblocksize);
/* window (if nonzero) */
if(!muted0 || !mutedC){
buffer+=qblocksize;
for(i=0;i<qblocksize*2;i++)
buffer[i] *= window[i]*scale;
}
if(!muted0 || !mutedC)
window_apply(buffer+qblocksize,window,scale,qblocksize);
}
static void freq_work(freq_class_setup *fc,
......@@ -307,12 +302,14 @@ static void freq_work(freq_class_setup *fc,
int i,j,ch=f->out.channels;
int have_feedback=0;
u_int32_t outactive=0;
f->cache_samples+=in->samples;
for(i=0;i<ch;i++){
int mutedC=mute_channel_muted(in->active,i);
int muted0=mute_channel_muted(f->mutemask0,i);
int muted1=mute_channel_muted(f->mutemask1,i);
int activeC=active[i] && !(muted0 && mutedC);
int active0=f->active0[i];
......@@ -386,9 +383,11 @@ static void freq_work(freq_class_setup *fc,
float *temp=out->data[i];
out->data[i]=f->cache1[i];
f->cache1[i]=temp;
outactive|= (f->mutemask1 & (1<<i));
}
}else{
float *w2=fc->window+input_size;
float *w2=fc->window;
float *l0=f->lap0[i];
float *l1=f->lap1[i];
float *lC=f->lapC[i];
......@@ -400,8 +399,10 @@ static void freq_work(freq_class_setup *fc,
float *tC=t0+input_size;
float *tN=tC+input_size;
outactive|= (1<<i);
if(!trans_active){
/* lap the cache into the trasform vector */
/* lap the cache into the transform vector */
fill_freq_buffer_helper(fc->fftwf_buffer,
fc->window,
f->cache0[i],in->data[i],
......@@ -411,10 +412,15 @@ static void freq_work(freq_class_setup *fc,
if(!activeP && !active1 && !active0){
/* previously in a bypassed/inactive state; the lapping cache
will need to be re-prepared */
memcpy(l1,c1,sizeof(*l1)*input_size);
for(j=0;j<input_size;j++)
l0[j]= c0[j]*w2[j];
if(muted1)
memset(l1,0,sizeof(*l1)*input_size);
else
memcpy(l1,c1,sizeof(*l1)*input_size);
if(muted0)
memset(l0,0,sizeof(*l0)*input_size);
else
for(j=0;j<input_size;j++)
l0[j]= c0[j]*w2[input_size-j];
memset(lC,0,sizeof(*lC)*input_size);
}
......@@ -478,7 +484,7 @@ static void freq_work(freq_class_setup *fc,
/* complete output linkage */
if(out){
out->active=f->mutemask1;
out->active=outactive;
f->out.samples=(f->cache_samples>input_size?input_size:f->cache_samples);
f->cache_samples-=f->out.samples;
}
......
......@@ -25,6 +25,7 @@
#include "feedback.h"
#include "bessel.h"
#include "limit.h"
#include "window.h"
extern int input_size;
extern int input_rate;
......@@ -38,6 +39,7 @@ typedef struct{
iir_state *iir;
iir_filter decay;
iir_filter limit;
int prev_active;
int initted;
......@@ -48,7 +50,7 @@ typedef struct{
limit_settings limitset;
limit_state limitstate;
float *window;
static float *window;
/* feedback! */
typedef struct limit_feedback{
......@@ -87,22 +89,19 @@ int limit_load(int ch){
for(i=0;i<ch;i++)
limitstate.out.data[i]=malloc(input_size*sizeof(**limitstate.out.data));
window=malloc(input_size*sizeof(*window));
for(i=0;i<input_size;i++){
window[i]=sin((i+.5)/input_size*M_PI*.5);
window[i]*=window[i];
}
window=window_get(1,input_size);
return(0);
}
static void filter_set(float msec,
int order,
iir_filter *filter){
float alpha;
float corner_freq= 500./msec;
alpha=corner_freq/input_rate;
filter->g=mkbessel(alpha,2,filter->c);
filter->g=mkbessel(alpha,order,filter->c);
filter->alpha=alpha;
filter->Hz=alpha*input_rate;
filter->ms=msec;
......@@ -141,7 +140,10 @@ time_linkage *limit_read(time_linkage *in){
float localatt;
float decayms=limitset.decay*.1;
if(decayms!=limitstate.decay.ms)filter_set(decayms,&limitstate.decay);
if(decayms!=limitstate.decay.ms){
filter_set(decayms,2,&limitstate.decay);
filter_set(decayms,1,&limitstate.limit);
}
if(in->samples==0){
limitstate.out.samples=0;
......@@ -191,8 +193,8 @@ time_linkage *limit_read(time_linkage *in){
prev_thresh+=thresh_add;
}
compute_iir_decayonly2(x,input_size,limitstate.iir+i,&limitstate.decay);
compute_iir_decay_limited(x,input_size,limitstate.iir+i,
&limitstate.decay,&limitstate.limit);
for(k=0;k<in->samples;k++)
x[k]=inx[k]*fromdB(-x[k]);
......
......@@ -24,6 +24,7 @@
#include "postfish.h"
#include "feedback.h"
#include "mix.h"
#include "window.h"
extern int input_ch;
extern int input_size;
......@@ -58,7 +59,7 @@ typedef struct{
mix_state ms;
/* this should be moved somewhere obvious/generic */
float *frame_window;
static float *window;
/* feedback! */
typedef struct limit_feedback{
......@@ -131,11 +132,7 @@ int mix_load(int outch){
for(i=0;i<outch;i++)
ms.out.data[i]=malloc(input_size*sizeof(**ms.out.data));
frame_window=malloc(input_size*sizeof(*frame_window));
for(i=0;i<input_size;i++){
frame_window[i]= sin( (i+.5)/input_size*M_PI*.5 );
frame_window[i]*=frame_window[i];
}
window=window_get(1,input_size);
return 0;
}
......@@ -182,20 +179,20 @@ static void mixwork(float *data,float *cacheP,float *cachePP,
/* current settings */
i=input_size*2+del;
while(i<input_size && offset<input_size){
out[offset]+=cachePP[i++]*att*frame_window[offset];
out[offset]+=cachePP[i++]*att*window[offset];
offset++;
}
i=input_size+del;
if(i<0)i=0;
while(i<input_size && offset<input_size){
out[offset]+=cacheP[i++]*att*frame_window[offset];
out[offset]+=cacheP[i++]*att*window[offset];
offset++;
}
i