src/ST2095_PinkNoise.cpp

   1 /*
   2 Copyright (c) 2015, John Hurst
   3 All rights reserved.
   4
   5 Redistribution and use in source and binary forms, with or without
   6 modification, are permitted provided that the following conditions
   7 are met:
   8 1. Redistributions of source code must retain the above copyright
   9    notice, this list of conditions and the following disclaimer.
  10 2. Redistributions in binary form must reproduce the above copyright
  11    notice, this list of conditions and the following disclaimer in the
  12    documentation and/or other materials provided with the distribution.
  13 3. The name of the author may not be used to endorse or promote products
  14    derived from this software without specific prior written permission.
  15
  16 THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  17 IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  18 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  19 IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  20 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  21 NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  22 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  23 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  24 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  25 THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26 */
  27 /*! \file    ST2095_PinkNoise.cpp
  28     \version $Id$
  29     \brief   Pink Noise filter and LCG generator
  30 */
  31
  32 #include "ST2095_PinkNoise.h"
  33
  34 //
  35 // This file is full of magic numbers.  Details behind the
  36 // selection of these values can be found in SMPTE ST 2095-1:2015.
  37 //
  38
  39 static ui32_t const C_rand_step = 52737;
  40 static float const C_max_ampl_32 = pow(2.0, 31) - 1.0;
  41 static float const C_max_peak = -9.5;     // Clipping Threshold in dB FS (+/-1.0 = 0 dB)
  42 static float const C_max_amp = pow(10.0, C_max_peak / 20.0);
  43
  44
  45 //
  46 ASDCP::LinearCongruentialGenerator::LinearCongruentialGenerator(const ui32_t sample_rate) : m_Seed(0)
  47 {
  48   ui32_t samples_per_period = 524288;
  49
  50   if ( sample_rate > 48000 )
  51     {
  52       samples_per_period = 1048576;
  53     }
  54
  55   m_RandMax = samples_per_period - 1;
  56   m_ScaleFactor = 2.0 / float(m_RandMax);
  57 }
  58
  59 //
  60 float
  61 ASDCP::LinearCongruentialGenerator::GetNextSample()
  62 {
  63   m_Seed = (1664525 * m_Seed + C_rand_step) & m_RandMax;
  64   float out = float(m_Seed) * m_ScaleFactor - 1.0;
  65   return out;
  66 }
  67
  68 //
  69 ASDCP::PinkFilter::PinkFilter(const i32_t sample_rate, float high_pass_fc, float low_pass_fc)
  70 {
  71   //  Disaster check: filters in order, low_pass_fc <= Nyquist
  72   assert(high_pass_fc < low_pass_fc);
  73   assert(low_pass_fc < sample_rate / 2.0);
  74
  75   // Calculate omegaT for matched Z transform highpass filters
  76   const float w0t = 2.0 * M_PI * high_pass_fc / sample_rate;
  77
  78   // Calculate k for bilinear transform lowpass filters
  79   const float k = tan(( 2.0 * M_PI * low_pass_fc / sample_rate ) / 2.0);
  80
  81   // precalculate k^2 (makes for a little bit cleaner code)
  82   const float k2 = k * k;
  83
  84   // Calculate biquad coefficients for bandpass filter components
  85   hp1_a1 = -2.0 * exp(-0.3826835 * w0t) * cos(0.9238795 * w0t);
  86   hp1_a2 = exp(2.0 * -0.3826835 * w0t);
  87   hp1_b0 = (1.0 - hp1_a1 + hp1_a2) / 4.0;
  88   hp1_b1 = -2.0 * hp1_b0;
  89   hp1_b2 = hp1_b0;
  90
  91   hp2_a1 = -2.0 * exp(-0.9238795 * w0t) * cos(0.3826835 * w0t);
  92   hp2_a2 = exp(2.0 * -0.9238795 * w0t);
  93   hp2_b0 = (1.0 - hp2_a1 + hp2_a2) / 4.0;
  94   hp2_b1 = -2.0 * hp2_b0;
  95   hp2_b2 = hp2_b0;
  96
  97   lp1_a1 = (2.0 * (k2 - 1.0)) / (k2 + (k / 1.306563) + 1.0);
  98   lp1_a2 = (k2 - (k / 1.306563) + 1.0) / (k2 + (k / 1.306563) + 1.0);
  99   lp1_b0 = k2 / (k2 + (k / 1.306563) + 1.0);
 100   lp1_b1 = 2.0 * lp1_b0;
 101   lp1_b2 = lp1_b0;
 102
 103   lp2_a1 = (2.0 * (k2 - 1.0)) / (k2 + (k / 0.541196) + 1.0);
 104   lp2_a2 = (k2 - (k / 0.541196) + 1.0) / (k2 + (k / 0.541196) + 1.0);
 105   lp2_b0 = k2 / (k2 + (k / 0.541196) + 1.0);
 106   lp2_b1 = 2.0 * lp2_b0;
 107   lp2_b2 = lp2_b0;
 108
 109   // Declare delay line variables for bandpass filter and initialize to zero
 110   hp1w1 = hp1w2 = hp2w1 = hp2w2 = 0.0;
 111   lp1w1 = lp1w2 = lp2w1 = lp2w2 = 0.0;
 112
 113   // Declare delay lines for pink filter network and initialize to zero
 114   lp1 = lp2 = lp3 = lp4 = lp5 = lp6 = 0.0;
 115 }
 116
 117
 118 //
 119 float
 120 ASDCP::PinkFilter::GetNextSample(const float white)
 121 {
 122   // Run pink filter; a parallel network of 1st order LP filters
 123   // Scaled for conventional RNG (need to rescale by sqrt(1/3) for MLS)
 124   lp1 = 0.9994551 * lp1 + 0.00198166688621989 * white;
 125   lp2 = 0.9969859 * lp2 + 0.00263702334184061 * white;
 126   lp3 = 0.9844470 * lp3 + 0.00643213710202331 * white;
 127   lp4 = 0.9161757 * lp4 + 0.01438952538362820 * white;
 128   lp5 = 0.6563399 * lp5 + 0.02698408541064610 * white;
 129   float pink = lp1 + lp2 + lp3 + lp4 + lp5 + lp6 + white * 0.0342675832159306;
 130   lp6 = white * 0.0088766118009356;
 131
 132   // Run bandpass filter; a series network of 4 biquad filters
 133   // Biquad filters implemented in Direct Form II
 134   float w = pink - hp1_a1 * hp1w1 - hp1_a2 * hp1w2;
 135   pink = hp1_b0 * w + hp1_b1 * hp1w1 + hp1_b2 * hp1w2;
 136   hp1w2 = hp1w1;
 137   hp1w1 = w;
 138
 139   w = pink - hp2_a1 * hp2w1 - hp2_a2 * hp2w2;
 140   pink = hp2_b0 * w + hp2_b1 * hp2w1 + hp2_b2 * hp2w2;
 141   hp2w2 = hp2w1;
 142   hp2w1 = w;
 143
 144   w = pink - lp1_a1 * lp1w1 - lp1_a2 * lp1w2;
 145   pink = lp1_b0 * w + lp1_b1 * lp1w1 + lp1_b2 * lp1w2;
 146   lp1w2 = lp1w1;
 147   lp1w1 = w;
 148
 149   w = pink - lp2_a1 * lp2w1 - lp2_a2 * lp2w2;
 150   pink = lp2_b0 * w + lp2_b1 * lp2w1 + lp2_b2 * lp2w2;
 151   lp2w2 = lp2w1;
 152   lp2w1 = w;
 153
 154   // Limit peaks to +/-C_max_amp
 155   if ( pink > C_max_amp )
 156     {
 157       pink = C_max_amp;
 158     }
 159   else if ( pink < -C_max_amp )
 160     {
 161       pink = -C_max_amp;
 162     }
 163
 164   return pink;
 165 }
 166
 167 //
 168 void
 169 ASDCP::ScalePackSample(float sample, byte_t* p, ui32_t word_size)
 170 {
 171   byte_t tmp_buf[4];
 172   Kumu::i2p<i32_t>(KM_i32_LE(sample * C_max_ampl_32), tmp_buf);
 173
 174   switch ( word_size )
 175     {
 176     case 4: *p++ = tmp_buf[0];
 177     case 3: *p++ = tmp_buf[1];
 178     case 2: *p++ = tmp_buf[2];
 179     case 1: *p++ = tmp_buf[3];
 180     }
 181 }
 182 //
 183 // end ST2095_PinkNoise.cpp
 184 //