2 Aug
2012
2 Aug
'12
8:31 p.m.
Am 02.08.2012 17:13, schrieb Charles Henry:
As I use this function quite often in reverb it's coded as inline to
avoid the call/return overhead:
/*****************************************************************************
*
* reverb.cpp
*
* Simulation of an electronic organ like Vox Continental
* with JACK MIDI input and JACK audio output
*
* Copyright (C) 2009,2010 Martin Homuth-Rosemann
*
* This program 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; version 2 of the License
*
* This program 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 this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
02110-1301 USA.
*
******************************************************************************/
#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <linux/types.h>
#include "reverb.h"
// modified JCRev (with output feedback)
// -> ccrma.stanford.edu/~jos/pasp/Schroeder_Reverberator_called_JCRev.html
// JCRev uses 3 all pass filtes in series
// and 4 parallel feed forward comb filters
// to replace the FFCF with IIR filters uncomment next line
// #define IIR
// three all pass filters
// ======================
// length of delay line
#define NA1 1051
#define NA2 337
#define NA3 113
// gain
#define GA1 0.707
#define GA2 0.707
#define GA3 0.707
// delay lines
static float ap1[NA1];
static float ap2[NA2];
static float ap3[NA3];
// position in delay line
static int ia1;
static int ia2;
static int ia3;
// four comb filters
// =================
#ifndef IIR
// length of delay line
#define NC1 4799
#define NC2 4999
#define NC3 5399
#define NC4 5801
#else
// length of delay line
#define NC1 479
#define NC2 499
#define NC3 539
#define NC4 581
#endif
// delay line
static float cf1[NC1];
static float cf2[NC2];
static float cf3[NC3];
static float cf4[NC4];
// position in delay line
static int ic1;
static int ic2;
static int ic3;
static int ic4;
//
// DENORMALS ARE EVIL
//
// 32 bit float
// SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMM
// E = 0, M != 0 -> denormal
// processing denormals uses lot of cpu.
// problem: an IIR feeds back 0.7*y.
// a value > 0 will decay until the smallest float is reached:
// 00000000000000000000000000000001
// multiplying with 0.7 and rounding (to nearest, default) gives again:
// 00000000000000000000000000000001
// this value circulates forever and consumes lot of cpu cycles :(
// even with "round to zero" - set in main() -
// it takes about 5 seconds until the denorm fades to zero...
//
// solution:
// "it's better to burn out than to fade away"
//
// denormals are zero
static inline float daz( float f )
{
// define an aliasing type to perform a "reinterpret cast"
typedef __u32 __attribute__ (( __may_alias__ )) u32bit;
if ( *( (u32bit*)&f ) & 0x7F000000 ) // E > 1 : normal.
return f;
else // E <= 1 : zero or _almost_ denormal
// (may become denormal with next operation)
return 0.0;
}
//
// reverb for one sample
//
float reverb( float xin )
{
static float yout = 0.0;
static float xv0, xv1, yv0, yv1;
float x, y;
// additional feedback
x = daz( xin/8 + yout/64 );
// three all pass filters
y = ap1[ia1];
ap1[ia1] = daz( GA1 * (x + y) );
x = y - x;
if ( ++ia1 >= NA1 )
ia1 = 0;
y = ap2[ia2];
ap2[ia2] = daz( GA2 * (x + y) );
x = y - x;
if ( ++ia2 >= NA2 )
ia2 = 0;
y = ap3[ia3];
ap3[ia3] = daz( GA3 * (x + y) );
x = y - x;
if ( ++ia3 >= NA3 )
ia3 = 0;
#ifndef IIR
// four feed forward comb filters
// gain
#define GC1 0.742
#define GC2 0.733
#define GC3 0.715
#define GC4 0.697
yout = 0;
yout += x + GC1 * cf1[ic1];
cf1[ic1] = x;
if ( ++ic1 >= NC1 )
ic1 = 0;
yout += x + GC2 * cf2[ic2];
cf2[ic2] = x;
if ( ++ic2 >= NC2 )
ic2 = 0;
yout += x + GC3 * cf3[ic3];
cf3[ic3] = x;
if ( ++ic3 >= NC3 )
ic3 = 0;
yout += x + GC4 * cf4[ic4];
cf4[ic4] = x;
if ( ++ic4 >= NC4 )
ic4 = 0;
#else
// four recursive comb filters
// gain
#define GC1 0.7
#define GC2 0.7
#define GC3 0.7
#define GC4 0.7
yout = 0.0;
y = cf1[ic1];
cf1[ic1] = daz( x + GC1 * y );
if ( ++ic1 >= NC1 )
ic1 = 0;
yout += y;
y = cf2[ic2];
cf2[ic2] = daz( x + GC2 * y );
if ( ++ic2 >= NC2 )
ic2 = 0;
yout += y;
y = cf3[ic3];
cf3[ic3] = daz( x + GC3 * y );
if ( ++ic3 >= NC3 )
ic3 = 0;
yout += y;
y = cf4[ic4];
cf4[ic4] = daz( x + GC4 * y );
if ( ++ic4 >= NC4 )
ic4 = 0;
yout += y;
#endif
// IIR LP filter 3000 Hz
xv0 = xv1;
xv1 = yout/6;
yv0 = yv1;
yout = yv1 = daz( xv0 + xv1 + 0.668 * yv0 );
return yout;
}
Hi Martin,
Can I pick your brain on how this works?
My biggest question is why to use the typedef__u32__attribute__ line
inside an inline function. Don't you only have to do this once? If
so, wouldn't you place that line outside the function?
You're right - but
the compiler didn't complain.
On Thu, Aug 2, 2012 at 7:39 AM, Martin Homuth-Rosemann
<linuxaudio(a)cryptomys.de> wrote:
No, it checks if the
exponent is >1 because we have
SEEE.EEEE.EMMM.MMMM.MMMM.MMMM.MMMM.MMMM
It's a special function to prevent denormals in recursive filters of a
reverb in my organ project connie. As my coeficients are .5 < coef < 1
the operation
out = in * coef;
may create denormals (E=0) if E is 1. So the code above changes to
out = daz( in ) * coef;
static inline float daz( float f )
If we
remove the "static" scope of this function, could we make a
header and re-use this code efficiently?
{
// define an aliasing type to perform a "reinterpret cast"
typedef __u32 __attribute__ (( __may_alias__ )) u32bit;
if ( *( (u32bit*)&f )& 0x7F000000 ) // E> 1 : normal.
Is there a
64-bit version of the same? This appears to work with a
32-bit unsigned int as a pointer, meaning that it's applicable to
32-bits. How does a "reinterpret cast" work?
What this does is check the exponent is non-zero, right? Doesn't this
also depend on the mantissa? for example, if the mantissa is 2^24-1
and the exponent is zero, is that number still a denormal? return
f;
else // E<= 1 : zero or _almost_ denormal
// (may become denormal with next operation)
return 0.0;
}
Ciao, Martin
Thanks for the code--any further explanation would be enlightening.
Chuck