________________________________ Interesting numbers ... I just don't understand the 32-bit pentium 3 numbers. Where is the reference from which the factors displayed in this table for this CPU are derived ? J.

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? On Thu, Aug 2, 2012 at 7:39 AM, Martin Homuth-Rosemann <linuxaudio(a)cryptomys.de> wrote: If we remove the "static" scope of this function, could we make a header and re-use this code efficiently? 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? Thanks for the code--any further explanation would be enlightening. Chuck

Am 02.08.2012 17:13, schrieb Charles Henry: You're right - but the compiler didn't complain. 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; 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; }

[David Robillard] Not to worry about the aliased pointer ... I think it's almost always a better idea to add an inaudible DC offset or a square wave at the block interval or at Nyquist, whatever, all less problematic than this (adorable but potentially slow) solution or relying on specific FTZ/DAZ-capable hardware to run your code. Tim

On 08/02/2012 10:56 PM, Harry van Haaren wrote: 1e-20 ( ~= signal at -400db ) should be inaudible enough. Denormals start < 1e-31 (~ -620db). You should to add it to each sample on the input, before the signal enters the effect. ciao, robin

On Fri, 2012-08-03 at 15:50 +0100, Harry van Haaren wrote: The problem is unique for Intel, in which case you can use: -msse -mfpmath=sse -ffast-math

[Harry van Haaren] Some useful numeric suggestions have already been made, and the why amounts to this: http://en.wikipedia.org/wiki/Single-precision_floating-point_format : """The minimum positive normal value [of a 32-bit float] is 2−126 ≈ 1.18 × 10−38""" 1.0f corresponding to 0dB makes 1.18e-38 equivalent to around -758dB. Roughly speaking, 16 bit samples cover around 90 dB (1./(1<<15)), 24 bits give you 138 dB (1./(1<<23)). For all intents and purposes involving current digital audio hardware and human beings at their current point of hearing evolution, anything past that is surely perfectly inaudible, so a sensible approach is to choose any amplitude between this and the above mentioned smallest normal float for your injected "normalising" signal, keeping in mind subsequent amplification factors in your signal path, and that amplitudes can decrease quickly in IIR filters, especially when you stack them (one feeding the next). Tim

On Thu, 2012-08-02 at 11:28 -0400, Paul Davis wrote: I mean set the flush-to-zero flag (and/or any other processor state stuff required). I know Ardour 'injects' non-zero data to fix denormals, but last I tinkered with this, setting that flag made the problems go away. This was many processor generations ago, though. -dr

On Thu, Aug 2, 2012 at 9:40 AM, Paul Davis &lt;paul(a)linuxaudiosystems.com&gt;wrote;wrote: Yes, please. I like this idea. Make it a JACK option. That way the black magic can be centralized instead of everybody having to copy code out of Ardour.