pa_dither.c (7128B)
1 /* 2 * $Id$ 3 * Portable Audio I/O Library triangular dither generator 4 * 5 * Based on the Open Source API proposed by Ross Bencina 6 * Copyright (c) 1999-2002 Phil Burk, Ross Bencina 7 * 8 * Permission is hereby granted, free of charge, to any person obtaining 9 * a copy of this software and associated documentation files 10 * (the "Software"), to deal in the Software without restriction, 11 * including without limitation the rights to use, copy, modify, merge, 12 * publish, distribute, sublicense, and/or sell copies of the Software, 13 * and to permit persons to whom the Software is furnished to do so, 14 * subject to the following conditions: 15 * 16 * The above copyright notice and this permission notice shall be 17 * included in all copies or substantial portions of the Software. 18 * 19 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 20 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 21 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 22 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR 23 * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF 24 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 25 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 26 */ 27 28 /* 29 * The text above constitutes the entire PortAudio license; however, 30 * the PortAudio community also makes the following non-binding requests: 31 * 32 * Any person wishing to distribute modifications to the Software is 33 * requested to send the modifications to the original developer so that 34 * they can be incorporated into the canonical version. It is also 35 * requested that these non-binding requests be included along with the 36 * license above. 37 */ 38 39 /** @file 40 @ingroup common_src 41 42 @brief Functions for generating dither noise 43 */ 44 45 #include "pa_types.h" 46 #include "pa_dither.h" 47 48 49 /* Note that the linear congruential algorithm requires 32 bit integers 50 * because it uses arithmetic overflow. So use PaUint32 instead of 51 * unsigned long so it will work on 64 bit systems. 52 */ 53 54 #define PA_DITHER_BITS_ (15) 55 56 57 void PaUtil_InitializeTriangularDitherState( PaUtilTriangularDitherGenerator *state ) 58 { 59 state->previous = 0; 60 state->randSeed1 = 22222; 61 state->randSeed2 = 5555555; 62 } 63 64 65 PaInt32 PaUtil_Generate16BitTriangularDither( PaUtilTriangularDitherGenerator *state ) 66 { 67 PaInt32 current, highPass; 68 69 /* Generate two random numbers. */ 70 state->randSeed1 = (state->randSeed1 * 196314165) + 907633515; 71 state->randSeed2 = (state->randSeed2 * 196314165) + 907633515; 72 73 /* Generate triangular distribution about 0. 74 * Shift before adding to prevent overflow which would skew the distribution. 75 * Also shift an extra bit for the high pass filter. 76 */ 77 #define DITHER_SHIFT_ ((sizeof(PaInt32)*8 - PA_DITHER_BITS_) + 1) 78 79 current = (((PaInt32)state->randSeed1)>>DITHER_SHIFT_) + 80 (((PaInt32)state->randSeed2)>>DITHER_SHIFT_); 81 82 /* High pass filter to reduce audibility. */ 83 highPass = current - state->previous; 84 state->previous = current; 85 return highPass; 86 } 87 88 89 /* Multiply by PA_FLOAT_DITHER_SCALE_ to get a float between -2.0 and +1.99999 */ 90 #define PA_FLOAT_DITHER_SCALE_ (1.0f / ((1<<PA_DITHER_BITS_)-1)) 91 static const float const_float_dither_scale_ = PA_FLOAT_DITHER_SCALE_; 92 93 float PaUtil_GenerateFloatTriangularDither( PaUtilTriangularDitherGenerator *state ) 94 { 95 PaInt32 current, highPass; 96 97 /* Generate two random numbers. */ 98 state->randSeed1 = (state->randSeed1 * 196314165) + 907633515; 99 state->randSeed2 = (state->randSeed2 * 196314165) + 907633515; 100 101 /* Generate triangular distribution about 0. 102 * Shift before adding to prevent overflow which would skew the distribution. 103 * Also shift an extra bit for the high pass filter. 104 */ 105 current = (((PaInt32)state->randSeed1)>>DITHER_SHIFT_) + 106 (((PaInt32)state->randSeed2)>>DITHER_SHIFT_); 107 108 /* High pass filter to reduce audibility. */ 109 highPass = current - state->previous; 110 state->previous = current; 111 return ((float)highPass) * const_float_dither_scale_; 112 } 113 114 115 /* 116 The following alternate dither algorithms (from musicdsp.org) could be 117 considered 118 */ 119 120 /*Noise shaped dither (March 2000) 121 ------------------- 122 123 This is a simple implementation of highpass triangular-PDF dither with 124 2nd-order noise shaping, for use when truncating floating point audio 125 data to fixed point. 126 127 The noise shaping lowers the noise floor by 11dB below 5kHz (@ 44100Hz 128 sample rate) compared to triangular-PDF dither. The code below assumes 129 input data is in the range +1 to -1 and doesn't check for overloads! 130 131 To save time when generating dither for multiple channels you can do 132 things like this: r3=(r1 & 0x7F)<<8; instead of calling rand() again. 133 134 135 136 int r1, r2; //rectangular-PDF random numbers 137 float s1, s2; //error feedback buffers 138 float s = 0.5f; //set to 0.0f for no noise shaping 139 float w = pow(2.0,bits-1); //word length (usually bits=16) 140 float wi= 1.0f/w; 141 float d = wi / RAND_MAX; //dither amplitude (2 lsb) 142 float o = wi * 0.5f; //remove dc offset 143 float in, tmp; 144 int out; 145 146 147 //for each sample... 148 149 r2=r1; //can make HP-TRI dither by 150 r1=rand(); //subtracting previous rand() 151 152 in += s * (s1 + s1 - s2); //error feedback 153 tmp = in + o + d * (float)(r1 - r2); //dc offset and dither 154 155 out = (int)(w * tmp); //truncate downwards 156 if(tmp<0.0f) out--; //this is faster than floor() 157 158 s2 = s1; 159 s1 = in - wi * (float)out; //error 160 161 162 163 -- 164 paul.kellett@maxim.abel.co.uk 165 http://www.maxim.abel.co.uk 166 */ 167 168 169 /* 170 16-to-8-bit first-order dither 171 172 Type : First order error feedforward dithering code 173 References : Posted by Jon Watte 174 175 Notes : 176 This is about as simple a dithering algorithm as you can implement, but it's 177 likely to sound better than just truncating to N bits. 178 179 Note that you might not want to carry forward the full difference for infinity. 180 It's probably likely that the worst performance hit comes from the saturation 181 conditionals, which can be avoided with appropriate instructions on many DSPs 182 and integer SIMD type instructions, or CMOV. 183 184 Last, if sound quality is paramount (such as when going from > 16 bits to 16 185 bits) you probably want to use a higher-order dither function found elsewhere 186 on this site. 187 188 189 Code : 190 // This code will down-convert and dither a 16-bit signed short 191 // mono signal into an 8-bit unsigned char signal, using a first 192 // order forward-feeding error term dither. 193 194 #define uchar unsigned char 195 196 void dither_one_channel_16_to_8( short * input, uchar * output, int count, int * memory ) 197 { 198 int m = *memory; 199 while( count-- > 0 ) { 200 int i = *input++; 201 i += m; 202 int j = i + 32768 - 128; 203 uchar o; 204 if( j < 0 ) { 205 o = 0; 206 } 207 else if( j > 65535 ) { 208 o = 255; 209 } 210 else { 211 o = (uchar)((j>>8)&0xff); 212 } 213 m = ((j-32768+128)-i); 214 *output++ = o; 215 } 216 *memory = m; 217 } 218 */