Commit bf77bf36 authored by Monty Montgomery's avatar Monty Montgomery

Ongoing filter work; add declipper (not inline yet, but avoid losing work)


git-svn-id: https://svn.xiph.org/trunk/postfish@5707 0101bb08-14d6-0310-b084-bc0e0c8e3800
parent c39f81d6
......@@ -5,8 +5,10 @@
CC=gcc
LD=gcc
SRC = main.c mainpanel.c multibar.c readout.c input.c output.c clippanel.c
OBJ = main.o mainpanel.o multibar.o readout.o input.o output.o clippanel.o
SRC = main.c mainpanel.c multibar.c readout.c input.c output.c clippanel.c declip.c \
reconstruct.c smallft.c
OBJ = main.o mainpanel.o multibar.o readout.o input.o output.o clippanel.o declip.o \
reconstruct.o smallft.o
GCF = `pkg-config --cflags gtk+-2.0` -DG_DISABLE_DEPRECATED -DGDK_DISABLE_DEPRECATED -DGTK_DISABLE_DEPRECATED -DGDK_PIXBUF_DISABLE_DEPRECATED
all:
......
/*
*
* postfish
*
* Copyright (C) 2002-2004 Monty and Xiph.Org
*
* 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 <math.h>
#include <sys/types.h>
#include "smallft.h"
#include "reconstruct.h"
#include <stdio.h>
extern int input_rate;
extern int input_ch;
extern int input_size;
void _analysis(char *base,int i,double *v,int n,int bark,int dB){
int j;
FILE *of;
char buffer[80];
sprintf(buffer,"%s_%d.m",base,i);
of=fopen(buffer,"w");
if(!of)perror("failed to open data dump file");
for(j=0;j<n;j++){
if(bark){
float b=toBark((4000.f*j/n)+.25);
fprintf(of,"%f ",b);
}else
fprintf(of,"%f ",(double)j);
if(dB){
if(j==0||j==n-1)
fprintf(of,"%f\n",todB(v[j]));
else{
fprintf(of,"%f\n",todB(hypot(v[j],v[j+1])));
j++;
}
}else{
fprintf(of,"%f\n",v[j]);
}
}
fclose(of);
}
/* accessed only in playback thread/setup */
static drft_lookup fft;
static int blocksize=0;
static int lopad=0,hipad=0;
static u_int32_t *widthlookup=0;
static double *window=0;
static double width=.5;
static double **lap=0;
static double **cache;
static int cache_samples;
static int fillstate=0; /* 0: uninitialized
1: normal
2: eof processed */
static time_linkage out;
/* accessed across threads */
static sig_atomic_t declip_active=0;
static sig_atomic_t declip_converge=2; /* 0=over, 1=full, 2=half, 3=partial, 4=approx */
static sig_atomic_t *chtrigger=0;
static sig_atomic_t *chactive=0;
static sig_atomic_t pending_blocksize=0;
/* called only by initial setup */
int declip_load(void){
int i;
chtrigger=malloc(input_ch*sizeof(*chtrigger));
chactive=malloc(input_ch*sizeof(*chactive));
for(i=0;i<input_ch;i++){
chtrigger[i]=0x80000000UL;
chactive[i]=1;
}
out.size=input_size;
out.channels=input_ch;
out.rate=input_rate;
out.data=malloc(input_ch*sizeof(*out.data));
for(i=0;i<input_ch;i++)
out.data[i]=malloc(input_size*sizeof(**out.data));
fillstate=0;
cache=malloc(input_ch*sizeof(*cache));
for(i=0;i<input_ch;i++)
cache[i]=malloc(input_size*sizeof(**cache));
lap=malloc(input_ch*sizeof(*lap));
for(i=0;i<input_ch;i++)
lap[i]=malloc(input_size*sizeof(**lap));
return(0);
}
int declip_setblock(int n){
if(n<32)return -1;
if(n>input_size*2)return -1;
pending_blocksize=n;
return 0;
}
int declip_setactive(int activep,int ch){
if(ch<0 || ch>=input_ch)return -1;
chactive[ch]=activep;
return 0;
}
int declip_settrigger(double trigger,int ch){
if(ch<0 || ch>=input_ch)return -1;
chtrigger[ch]=rint(trigger*(1.*0x80000000UL));
return 0;
}
/* called only in playback thread */
int declip_reset(void){
/* reset cached pipe state */
fillstate=0;
return 0;
}
static void sliding_bark_average(double *f,double *w, int n,double width){
int i=0,j;
double acc=0.,del=0.;
double sec[hipad+1];
memset(sec,0,sizeof(sec));
{
double bark=toBark(0.);
int hi=widthlookup[0]>>16;
int lo=widthlookup[0]&(0xffff);
double del=fabs(f[0])/(lo-hi);
double hidel=del/(-hi+lopad);
double lodel=del/(lo-lopad);
sec[hi]+=hidel;
sec[lopad]-=hidel;
sec[lopad]-=lodel;
sec[lo]+=lodel;
}
for(i=1;i<n/2;i++){
double bark=toBark(44100.*i/n);
int hi=widthlookup[i]>>16;
int lo=widthlookup[i]&(0xffff);
double del=hypot(f[(i<<1)-1],f[i<<1])/(lo-hi);
double hidel=del/((i-hi+lopad));
double lodel=del/((lo-i-lopad));
sec[hi]+=hidel;
sec[i+lopad]-=hidel;
sec[i+lopad]-=lodel;
sec[lo]+=lodel;
}
for(i=0;i<lopad;i++){
del+=sec[i];
acc+=del;
}
w[0]=1./(acc*acc);
del+=sec[lopad];
acc+=del;
for(i=1;i<n/2;i++){
w[(i<<1)-1]=w[i<<1]=1./(acc*acc);
del+=sec[i+lopad];
acc+=del;
}
w[n-1]=w[n-2];
}
static void declip(double *data,double *lap,double *out,
int blocksize,unsigned long trigger){
double freq[blocksize];
int flag[blocksize];
double triggerlevel=trigger*(1./0x80000000UL);
double epsilon=1e-12;
int iterbound=blocksize,i,j,count=0;
for(i=0;i<blocksize/8;i++)flag[i]=0;
for(;i<blocksize*7/8;i++){
flag[i]=0;
if(data[i]>=trigger || data[i]<=-trigger){
flag[i]=1;
count++;
}
}
for(;i<blocksize;i++)flag[i]=0;
for(i=0;i<blocksize;i++)data[i]*=window[i];
memcpy(freq,data,sizeof(freq));
drft_forward(&fft,freq);
sliding_bark_average(freq,freq,blocksize,width);
switch(declip_converge){
case 0:
epsilon=1e-12;
iterbound=blocksize*2;
break;
case 1:
epsilon=1e-8;
iterbound=count;
break;
case 2:
epsilon=1e-6;
iterbound=count/2;
break;
case 3:
epsilon=1e-5;
iterbound=count/4;
break;
case 4:
epsilon=1e-3;
iterbound=count/8;
break;
}
if(iterbound<20)iterbound=20;
reconstruct(&fft,data,freq,flag,epsilon*count,iterbound,blocksize);
if(out)
for(i=0;i<blocksize/2;i++)
out[i]=lap[i]+data[i]*window[i];
for(i=blocksize/2,j=0;i<blocksize;i++)
lap[j]=data[i]*window[i];
}
/* called only by playback thread */
time_linkage *declip_read(time_linkage *in){
int i;
double work[blocksize];
if(pending_blocksize!=blocksize){
if(blocksize){
free(widthlookup);
free(window);
drft_clear(&fft);
}
blocksize=pending_blocksize;
widthlookup=malloc((blocksize>>1)*sizeof(*widthlookup));
for(i=0;i<blocksize/2;i++){
double bark=toBark(input_rate*i/blocksize);
int hi=rint(fromBark(bark-width)*blocksize/input_rate)-1+lopad;
int lo=rint(fromBark(bark+width)*blocksize/input_rate)+1+lopad;
if(hi<0 || lo<0 || hi>65535 || lo<65535) return 0;
widthlookup[i]=(hi<<16)+lo;
}
lopad=1-rint(fromBark(toBark(0.)-width)*blocksize/input_rate);
hipad=rint(fromBark(toBark(input_rate*.5)+width)*blocksize/input_rate)+lopad;
drft_init(&fft,blocksize);
window=malloc(blocksize*sizeof(*window));
for(i=0;i<blocksize/8;i++) window[i]=0.;
for(;i<blocksize*3/8;i++) window[i]=sin( (double)(i-blocksize/8)/blocksize*M_PIl );
for(;i<blocksize*5/8;i++) window[i]=1.;
for(;i<blocksize*7/8;i++) window[i]=sin( (double)(blocksize*7/8-i)/blocksize*M_PIl );
for(;i<blocksize;i++) window[i]=0.;
for(i=0;i<blocksize;i++) window[i]*=window[i];
for(i=0;i<blocksize;i++) window[i]=sin(window[i]*M_PIl);
}
switch(fillstate){
case 0: /* prime the lapping and cache */
for(i=0;i<input_ch;i++){
int j;
double *temp=in->data[i];
if(chactive[i] && declip_active){
memset(work,0,sizeof(*work)*blocksize/2);
memcpy(work+blocksize/2,temp,sizeof(*work)*blocksize/2);
declip(work,lap[i],0,blocksize,chtrigger[i]);
}else{
for(j=0;j<blocksize/2;j++)
lap[i][j]=window[j+blocksize/2]*temp[j];
}
memset(cache[i],0,sizeof(**cache)*input_size);
in->data[i]=cache[i];
cache[i]=temp;
}
cache_samples=in->samples;
fillstate=1;
if(in->samples==in->size)return 0;
in->samples=0;
/* fall through */
case 1: /* nominal processing */
for(i=0;i<input_ch;i++){
double *temp=cache[i];
int j;
if(chactive[i] && declip_active){
for(j=0;j+blocksize<out.size;j+=blocksize/2){
memcpy(work,temp+j,sizeof(*work)*blocksize);
declip(work,lap[i],out.data[i]+j,blocksize,chtrigger[i]);
}
memcpy(work,temp+j,sizeof(*work)*blocksize/2);
memcpy(work+blocksize/2,in->data[i],sizeof(*work)*blocksize/2);
declip(work,lap[i],out.data[i]+j,blocksize,chtrigger[i]);
}else{
memcpy(out.data[i],temp,out.size*sizeof(*temp));
for(j=0;j<blocksize/2;j++)
lap[i][j]=window[j+blocksize/2]*temp[j];
}
in->data[i]=cache[i];
cache[i]=temp;
}
out.samples=cache_samples;
cache_samples=in->samples;
if(out.samples<out.size)fillstate=2;
break;
case 2: /* we've pushed out EOF already */
return 0;
}
return &out;
}
/*
*
* postfish
*
* Copyright (C) 2002-2004 Monty and Xiph.Org
*
* 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.
*
*
*/
extern int declip_load(void);
extern int declip_setblock(int n);
extern int declip_setactive(int activep,int ch);
extern int declip_settrigger(double trigger,int ch);
extern int declip_reset(void);
extern time_linkage *declip_read(time_linkage *in);
......@@ -31,10 +31,11 @@ static off_t cursor=0;
sig_atomic_t loop_active;
int seekable;
static int input_rate;
int input_rate;
int input_ch;
static int inbytes;
static int signp;
int input_size;
typedef struct {
FILE *f;
......@@ -308,7 +309,7 @@ int input_load(int n,char *list[]){
}
}
out.size=2048;
input_size=out.size=2048;
input_ch=out.channels=ch;
input_rate=out.rate=rate;
out.data=malloc(sizeof(*out.data)*ch);
......
......@@ -47,6 +47,7 @@ Carsten Bormann
*********************************************************************/
#define _ISOC99_SOURCE
#define _GNU_SOURCE
#include <stdlib.h>
#include <string.h>
......@@ -66,23 +67,6 @@ static void autocorr(double *data,double *aut,int n,int m){
}
}
static void autocorr_incomplete(double *data,double *aut,int n,int m){
int i,j;
/* autocorrelation, p+1 lag coefficients */
j=m+1;
while(j--){
int count=0;
double d=0.; /* double needed for accumulator depth */
for(i=j;i<n;i++){
if(!isnan(data[i]) && !isnan(data[i-j])){
d+=data[i]*data[i-j];
count++;
}
}
aut[j]=d;
}
}
/* Autocorrelation LPC coeff generation algorithm invented by
N. Levinson in 1947, modified by J. Durbin in 1959. */
......@@ -91,14 +75,14 @@ static double levinson_durbin(double *aut,double *lpc,int m){
double error=aut[0];
int i,j;
if(error==0){
memset(lpc,0,m*sizeof(*lpc));
return 0.;
}
for(i=0;i<m;i++){
double r= -aut[i+1];
if(error==0){
memset(lpc,0,m*sizeof(*lpc));
return 0.;
}
/* Sum up this iteration's reflection coefficient; note that in
Vorbis we don't save it. If anyone wants to recycle this code
and needs reflection coefficients, save the results of 'r' from
......@@ -135,12 +119,6 @@ double lpc_from_data(double *data,double *lpc,int n,int m){
return levinson_durbin(aut,lpc,m);
}
double lpc_from_incomplete_data(double *data,double *lpc,int n,int m){
double *aut=alloca(sizeof(*aut)*(m+1));
autocorr_incomplete(data,aut,n,m);
return levinson_durbin(aut,lpc,m);
}
void lpc_extrapolate(double *lpc,double *prime,int m,
double *data,int n){
......@@ -165,15 +143,14 @@ void lpc_extrapolate(double *lpc,double *prime,int m,
p=m;
for(j=0;j<m;j++)
y-=work[o++]*lpc[--p];
data[i]=work[o]=y+data[i];
data[i]=work[o]=y;
}
}
void lpc_subtract(double *lpc,double *prime,int m,
double *data,int n){
long i,j,o,p;
void lpc_extrapolateB(double *lpc,double *prime,int m,
double *data,int n){
long i,j,o,p;
double y;
double *work=alloca(sizeof(*work)*(m+n));
......@@ -182,7 +159,7 @@ void lpc_subtract(double *lpc,double *prime,int m,
work[i]=0.;
else
for(i=0;i<m;i++)
work[i]=prime[i];
work[i]=prime[m-i-1];
for(i=0;i<n;i++){
y=0;
......@@ -190,148 +167,11 @@ void lpc_subtract(double *lpc,double *prime,int m,
p=m;
for(j=0;j<m;j++)
y-=work[o++]*lpc[--p];
work[o]=data[i];
data[i]-=y;
}
}
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#define todB(x) ((x)==0?-140.f:log((x)*(x))*4.34294480f)
#define toBARK(n) (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n))
void _analysis(char *base,int i,double *v,int n,int bark,int dB){
int j;
FILE *of;
char buffer[80];
sprintf(buffer,"%s_%d.m",base,i);
of=fopen(buffer,"w");
if(!of)perror("failed to open data dump file");
for(j=0;j<n;j++){
if(bark){
float b=toBARK((4000.f*j/n)+.25);
fprintf(of,"%f ",b);
}else
fprintf(of,"%f ",(double)j);
if(dB){
fprintf(of,"%f\n",todB(v[j]));
}else{
fprintf(of,"%f\n",v[j]);
}
data[n-i-1]=work[o]=y;
}
fclose(of);
}
void lpc_to_curve(double *curve,double *lpc,int n,int m,double amp){
int i;
drft_lookup fft;
double *work=alloca(sizeof(*work)*n*2);
memset(work,0,sizeof(double)*n*2);
for(i=0;i<m;i++){
work[i*2+1]=lpc[i]/(4*amp);
work[i*2+2]=-lpc[i]/(4*amp);
}
drft_init(&fft,n*2);
drft_backward(&fft,work); /* reappropriated ;-) */
drft_clear(&fft);
{
int n2=n*2;
double unit=1./amp;
curve[0]=(1./(work[0]*2+unit));
for(i=1;i<n;i++){
double real=(work[i]+work[n2-i]);
double imag=(work[i]-work[n2-i]);
double a = real + unit;
curve[i] = 1.0 / hypot(a, imag);
}
}
}
#define N 1024
#define M 32
#define READ (N+M)
signed char readbuffer[READ*4+44];
main(){
static int seq;
int eos=0,ret;
int i, founddata;
readbuffer[0] = '\0';
for (i=0, founddata=0; i<30 && ! feof(stdin) && ! ferror(stdin); i++)
{
fread(readbuffer,1,2,stdin);
if ( ! strncmp((char*)readbuffer, "da", 2) )
{
founddata = 1;
fread(readbuffer,1,6,stdin);
break;
}
}
while(!eos){
double samples[READ];
double curve[N];
double lpc[M];
long i;
long bytes=fread(readbuffer,1,READ*4,stdin); /* stereo hardwired here */
if(bytes<READ*4)break;
/* uninterleave samples */
for(i=0;i<bytes/4;i++)
samples[i]=
((readbuffer[i*4+1]<<8)|
(0x00ff&(int)readbuffer[i*4]))/32768. +
((readbuffer[i*4+3]<<8)|
(0x00ff&(int)readbuffer[i*4+2]))/32768.f;
_analysis("pcm",seq,samples,READ,0,0);
lpc_from_data(samples,lpc,READ,M);
_analysis("lpc",seq,lpc,M,0,0);
lpc_to_curve(curve,lpc,N,M,1.0);
_analysis("H",seq,curve,N,1,1);
/* kill off samples */
{
double max=0;
for(i=0;i<READ;i++)
if(samples[i]>max)max=samples[i];
for(i=0;i<READ;i++){
if(samples[i]>max*.50)samples[i]=max*.50;
if(samples[i]<-max*.50)samples[i]=-max*.50;
}
}
_analysis("pcmi",seq,samples,READ,0,0);
lpc_from_incomplete_data(samples,lpc,READ,M);
_analysis("lpci",seq,lpc,M,0,0);
lpc_to_curve(curve,lpc,N,M,1.0);
_analysis("Hi",seq,curve,N,1,1);
lpc_subtract(lpc,samples,M,samples+M,N);
_analysis("resi",seq,samples,READ,0,0);
seq++;
}
return 0;
}
......@@ -900,9 +900,9 @@ void mainpanel_go(int argc,char *argv[], int ch){
int i;
memset(&p,0,sizeof(p));
gtk_rc_add_default_file("/etc/postfish/postfishrc");
gtk_rc_add_default_file("/etc/postfish/postfish-gtkrc");
if(homedir){
char *rcfile="/.postfishrc";
char *rcfile="/.postfish-gtkrc";
char *homerc=calloc(1,strlen(homedir)+strlen(rcfile)+1);
strcat(homerc,homedir);
strcat(homerc,rcfile);
......
......@@ -20,123 +20,270 @@
*
*
*/
#define _GNU_SOURCE
#include <stdio.h>
#include "smallft.h"
/* Implment a linear system least-squares minimizer. Typical
implementations are inefficient due to our need for a weighting
vector that changes from frame to frame, eliminating our ability to
presolve the first- and second-order differential matricies. This
alternate to Newton's Method allows fast worst-case O(N^2) direct
solution of the least-squares fit. */
/* if no 'x' vector, assume zeroes. If no slope vector, assume unit slope */