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., 59 Temple Place, Suite 330,
  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<float> 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<float>::unwarpf(20.6f, sr);
  46         float f2 = biquad_coeffs<float>::unwarpf(107.7f, sr);
  47         float f3 = biquad_coeffs<float>::unwarpf(738.f, sr);
  48         float f4 = biquad_coeffs<float>::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 };
  89
  90 #endif