[SCM] calf/master: + A-Weighting: use correct curve (and calculate it semi-sanely), add comments

[js at users.alioth.debian.org]

The following commit has been merged in the master branch:
commit 80916019790742c6d6a1e02998f3ae18059dea31
Author: Krzysztof Foltman <wdev at foltman.com>
Date:   Mon Nov 3 23:40:42 2008 +0000

    + A-Weighting: use correct curve (and calculate it semi-sanely), add comments

diff --git a/src/benchmark.cpp b/src/benchmark.cpp
index a852836..479c8e1 100644
--- a/src/benchmark.cpp
+++ b/src/benchmark.cpp
@@ -343,10 +343,11 @@ void eq_calc()
 void aweighting_calc()
 {
     aweighter aw;
-    aw.set(48000);
+    float fs = 44100;
+    aw.set(fs);
     for (int i = 10; i < 20000; i += 10)
     {
-        printf("%d %f\n", i, 20*log10(aw.freq_gain(i * 1.0, 48000)));
+        printf("%d %f\n", i, 20*log10(aw.freq_gain(i * 1.0, fs)));
     }
 }
 
diff --git a/src/calf/biquad.h b/src/calf/biquad.h
index d90027a..f98bcc1 100644
--- a/src/calf/biquad.h
+++ b/src/calf/biquad.h
@@ -191,17 +191,21 @@ public:
         if (freq>sr*0.49) freq=(float)(sr*0.49);
         return (float)(tan(M_PI*freq/sr));
     }
-    /// convert analog frequency-corresponding pole value to digital
-    static inline float unwarp(float freq, float sr)
+    /// convert analog angular frequency value to digital
+    static inline float unwarp(float omega, float sr)
     {
         float T = 1.0 / sr;
-        return (2 / T) * atan(freq * T / 2);
+        return (2 / T) * atan(omega * T / 2);
     }
-    /// convert analog frequency-corresponding pole value to digital
-    static inline float unwarpf(float coeff, float sr)
+    /// convert analog filter time constant to digital counterpart
+    static inline float unwarpf(float t, float sr)
     {
-        float T = 1.0 / sr;
-        return 1.0 / ((2 / T) * atan((1.0 / coeff) * T / 2));
+        // this is most likely broken and works by pure accident!
+        float omega = 1.0 / t;
+        omega = unwarp(omega, sr);
+        // I really don't know why does it have to be M_PI and not 2 * M_PI!
+        float f = M_PI / omega;
+        return f / sr;
     }
     /// set digital filter parameters based on given analog filter parameters
     void set_bilinear(float aa0, float aa1, float aa2, float ab0, float ab1, float ab2)
diff --git a/src/calf/loudness.h b/src/calf/loudness.h
index c9dccf3..3dc5a97 100644
--- a/src/calf/loudness.h
+++ b/src/calf/loudness.h
@@ -31,32 +31,37 @@ class aweighter {
 public:
     biquad_d2<float> bq1, bq2, bq3;
     
+    /// Produce one output sample from one input sample
     float process(float sample)
     {
         return bq1.process(bq2.process(bq3.process(sample)));
     }
     
+    /// Set sample rate (updates filter coefficients)
     void set(float sr)
     {
-        // lowpass : H(s) = 1 / (1 + st) = 1 / (1 + s/w) = 1 / (1 + s/(2piF))
-        // wc = 2pi * fc
-
-        // This is not done yet - I need to finish redoing the coeffs properly
-        float f1 = 129.4f / sr; 
-        float f2 = 676.7f / sr; 
-        float f3 = 4636.f / sr; 
-        float f4 = 76655.f / sr;
-        /*
-        float f1 = biquad_coeffs<float>::unwarpf(129.4f, sr);
-        float f2 = biquad_coeffs<float>::unwarpf(676.7f, sr);
-        float f3 = biquad_coeffs<float>::unwarpf(4636.f, sr);
-        float f4 = biquad_coeffs<float>::unwarpf(76655.f, sr);
-        */
+        // analog coeffs taken from: http://www.diracdelta.co.uk/science/source/a/w/aweighting/source.html
+        // first we need to adjust them by doing some obscene sort of reverse pre-warping (a broken one, too!)
+        float f1 = biquad_coeffs<float>::unwarpf(20.6f, sr);
+        float f2 = biquad_coeffs<float>::unwarpf(107.7f, sr);
+        float f3 = biquad_coeffs<float>::unwarpf(738.f, sr);
+        float f4 = biquad_coeffs<float>::unwarpf(12200.f, sr);
+        // then map s domain to z domain using bilinear transform
+        // note: f1 and f4 are double poles
         bq1.set_bilinear(0, 0, 1, f1*f1, 2 * f1, 1);
-        bq2.set_bilinear(4.0, 0, 0, f2*f3, f2 + f3, 1);
+        bq2.set_bilinear(1, 0, 0, f2*f3, f2 + f3, 1);
         bq3.set_bilinear(0, 0, 1, f4*f4, 2 * f4, 1);
+        // the coeffs above give non-normalized value, so it should be normalized to produce 0dB at 1 kHz
+        // find actual gain
+        float gain1kHz = freq_gain(1000.0, sr);
+        // divide one filter's x[n-m] coefficients by that value
+        float gc = 1.0 / gain1kHz;
+        bq1.a0 *= gc;
+        bq1.a1 *= gc;
+        bq1.a2 *= gc;
     }
     
+    /// Reset to zero if at risk of denormals
     void sanitize()
     {
         bq1.sanitize();
@@ -64,6 +69,7 @@ public:
         bq3.sanitize();
     }
     
+    /// Reset state to zero
     void reset()
     {
         bq1.reset();
@@ -71,6 +77,7 @@ public:
         bq3.reset();
     }
     
+    /// Gain and a given frequency
     float freq_gain(float freq, float sr)
     {
         return bq1.freq_gain(freq, sr) * bq2.freq_gain(freq, sr) * bq3.freq_gain(freq, sr);

-- 
calf audio plugins packaging