src/calf/loudness.h

   1 /* Calf DSP Library
   2  * A-weighting filter for
   3  * Copyright (C) 2001-2007 Krzysztof Foltman
   4  *
   5  * Most of code in this file is based on freely
   6  * available other work of other people (filter equations).
   7  *
   8  * This program is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Lesser General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2 of the License, or (at your option) any later version.
  12  *
  13  * This program is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Lesser General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU Lesser General
  19  * Public License along with this program; if not, write to the
  20  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  21  * Boston, MA 02111-1307, USA.
  22  */
  23 #ifndef __CALF_LOUDNESS_H
  24 #define __CALF_LOUDNESS_H
  25
  26 #include "biquad.h"
  27
  28 namespace dsp {
  29
  30 class aweighter {
  31 public:
  32     biquad_d2 bq1, bq2, bq3;
  33
  34     /// Produce one output sample from one input sample
  35     float process(float sample)
  36     {
  37         return bq1.process(bq2.process(bq3.process(sample)));
  38     }
  39
  40     /// Set sample rate (updates filter coefficients)
  41     void set(float sr)
  42     {
  43         // analog coeffs taken from: http://www.diracdelta.co.uk/science/source/a/w/aweighting/source.html
  44         // first we need to adjust them by doing some obscene sort of reverse pre-warping (a broken one, too!)
  45         float f1 = biquad_coeffs::unwarpf(20.6f, sr);
  46         float f2 = biquad_coeffs::unwarpf(107.7f, sr);
  47         float f3 = biquad_coeffs::unwarpf(738.f, sr);
  48         float f4 = biquad_coeffs::unwarpf(12200.f, sr);
  49         // then map s domain to z domain using bilinear transform
  50         // note: f1 and f4 are double poles
  51         bq1.set_bilinear(0, 0, 1, f1*f1, 2 * f1, 1);
  52         bq2.set_bilinear(1, 0, 0, f2*f3, f2 + f3, 1);
  53         bq3.set_bilinear(0, 0, 1, f4*f4, 2 * f4, 1);
  54         // the coeffs above give non-normalized value, so it should be normalized to produce 0dB at 1 kHz
  55         // find actual gain
  56         float gain1kHz = freq_gain(1000.0, sr);
  57         // divide one filter's x[n-m] coefficients by that value
  58         float gc = 1.0 / gain1kHz;
  59         bq1.a0 *= gc;
  60         bq1.a1 *= gc;
  61         bq1.a2 *= gc;
  62     }
  63
  64     /// Reset to zero if at risk of denormals
  65     void sanitize()
  66     {
  67         bq1.sanitize();
  68         bq2.sanitize();
  69         bq3.sanitize();
  70     }
  71
  72     /// Reset state to zero
  73     void reset()
  74     {
  75         bq1.reset();
  76         bq2.reset();
  77         bq3.reset();
  78     }
  79
  80     /// Gain and a given frequency
  81     float freq_gain(float freq, float sr)
  82     {
  83         return bq1.freq_gain(freq, sr) * bq2.freq_gain(freq, sr) * bq3.freq_gain(freq, sr);
  84     }
  85
  86 };
  87
  88 class riaacurve {
  89 public:
  90     biquad_d2 r1;
  91     biquad_d2 brickw;
  92
  93     /// Produce one output sample from one input sample
  94     float process(float sample)
  95     {
  96         return brickw.process(r1.process(sample));
  97     }
  98
  99     /// Set sample rate (updates filter coefficients)
 100     void set(float sr, int mode, int type)
 101     {
 102         float i,j,k,g,a0,a1,a2,b1,b2,tau1,tau2,tau3;
 103         switch(type) {
 104             case 0: //"Columbia"
 105                 i = 100.f;
 106                 j = 500.f;
 107                 k = 1590.f;
 108                 break;
 109             case 1: //"EMI"
 110                 i = 70.f;
 111                 j = 500.f;
 112                 k = 2500.f;
 113                 break;
 114             case 2: //"BSI(78rpm)"
 115                 i = 50.f;
 116                 j = 353.f;
 117                 k = 3180.f;
 118                 break;
 119             case 4: //"CD Mastering"
 120                 tau1 = 0.000050f;
 121                 tau2 = 0.000015f;
 122                 tau3 = 0.0000001f;// 1.6MHz out of audible range for null impact
 123                 i = 1.f / (2.f * M_PI * tau1);
 124                 j = 1.f / (2.f * M_PI * tau2);
 125                 k = 1.f / (2.f * M_PI * tau3);
 126             break;
 127             case 3: //"RIAA"
 128             default:
 129                 tau1 = 0.003180f;
 130                 tau2 = 0.000318f;
 131                 tau3 = 0.000075f;
 132                 i = 1.f / (2.f * M_PI * tau1);
 133                 j = 1.f / (2.f * M_PI * tau2);
 134                 k = 1.f / (2.f * M_PI * tau3);
 135             break;
 136         }
 137
 138         i *= 2.f * M_PI;
 139         j *= 2.f * M_PI;
 140         k *= 2.f * M_PI;
 141
 142         float t = 1.f / sr;
 143
 144         if (mode == 0) { //Reproduction
 145             g = 1.f / (4.f+2.f*i*t+2.f*k*t+i*k*t*t);
 146             a0 = (2.f*t+j*t*t)*g;
 147             a1 = (2.f*j*t*t)*g;
 148             a2 = (-2.f*t+j*t*t)*g;
 149             b1 = (-8.f+2.f*i*k*t*t)*g;
 150             b2 = (4.f-2.f*i*t-2.f*k*t+i*k*t*t)*g;
 151         } else {  //Production
 152             g = 1.f / (2.f*t+j*t*t);
 153             a0 = (4.f+2.f*i*t+2.f*k*t+i*k*t*t)*g;
 154             a1 = (-8.f+2.f*i*k*t*t)*g;
 155             a2 = (4.f-2.f*i*t-2.f*k*t+i*k*t*t)*g;
 156             b1 = (2.f*j*t*t)*g;
 157             b2 = (-2.f*t+j*t*t)*g;
 158         }
 159
 160         r1.sanitize();
 161
 162         //swap a1 b1, a2 b2
 163         r1.set_bilinear_direct(a0, a1, a2, b1, b2);
 164
 165         // the coeffs above give non-normalized value, so it should be normalized to produce 0dB at 1 kHz
 166         // find actual gain
 167         float gain1kHz = freq_gain(1000.0, sr);
 168         // divide one filter's x[n-m] coefficients by that value
 169         float gc = 1.0 / gain1kHz;
 170         r1.a0 *= gc;
 171         r1.a1 *= gc;
 172         r1.a2 *= gc;
 173         r1.sanitize();
 174
 175         float cutfreq = std::min(0.45f * sr, 21000.f);
 176         brickw.set_lp_rbj(cutfreq, 1.f, sr, 1.f);
 177         brickw.sanitize();
 178     }
 179
 180     /// Reset to zero if at risk of denormals
 181     void sanitize()
 182     {
 183         r1.sanitize();
 184     }
 185
 186     /// Reset state to zero
 187     void reset()
 188     {
 189         r1.reset();
 190     }
 191
 192     /// Gain and a given frequency
 193     float freq_gain(float freq, float sr) const
 194     {
 195         return r1.freq_gain(freq, sr);
 196     }
 197
 198 };
 199
 200 };
 201
 202 #endif