Topic: libebur128 - (yet another) EBU R 128 implementation (Read 117236 times)
0 Members and 1 Guest are viewing this topic.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #25 – 2011-01-16 17:49:16
I've thought long about this. There are five filter coefficients for a normalized, second-order IIR filter like the ones in BS.1770-1. Such a filter depends on five design parameters: High-pass gain, band-pass gain, low-pass gain, Q factor, and angular frequency (K = tan(pi * f_c / f_s), f_s is the sampling frequency, f_c the frequency "where it happens").

Sounds very similar to what I've found during my Web search:
[blockquote]http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt[/blockquote]
If one knows how the filter was designed, one can solve the resulting system of equations to get the parameters.

Yesterday I had some fun with Maxima and LaTeX, and wrote a little paper that explains what I did.

This paper was very helpful. On page 3 there are formulas for b_0, b_1, b_2, a_1, a_2, obtained with bilinear transform.

Many thanks for providing the links.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #26 – 2011-01-16 19:24:53
On the other hand I'm also interested in Raiden's method for calculating the filter coefficients for various sample rates on the fly. I've searched the whole net several times, ask the question in the forum, unfortunately no answer yet.

It's known as the bilinear transform.  It maps analog filter coefficients to digital ones, warping the frequencies to fit properly within a bandlimited system.  It's a common technique used for filters that require dynamically calculated filter coefficients (such as an EQ).

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #27 – 2011-01-16 19:32:57
On the other hand I'm also interested in Raiden's method for calculating the filter coefficients for various sample rates on the fly. I've searched the whole net several times, ask the question in the forum, unfortunately no answer yet.

It's known as the bilinear transform.  It maps analog filter coefficients to digital ones, warping the frequencies to fit properly within a bandlimited system.  It's a common technique used for filters that require dynamically calculated filter coefficients (such as an EQ).

Thanks a lot. I don't have any experience with DSP so far, it's a long way to go ...

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #28 – 2011-01-16 20:40:14
If you want to have my opinion: I'm pretty sure the true peak values can be estimated quite accurately without resampling, so please don't consider adding resampling in the future.

R128 requires oversampling to compute dBTP. If you don't need dBTP, you don't need to oversample. If you're using this library to do Replay Gain tagging, you want to know sample peak because that is what Replay Gain specifies you should tag.  The value of dBTP metering depends on the application. On one hand, modern DACs don't have a problem decoding signals greater than 0 dBTP. On the other hand, staying below 0 dBTP is good practice; it keeps you from going over full scale through sample rate conversion.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #29 – 2011-01-16 20:46:48
Yesterday I had some fun with Maxima and LaTeX, and wrote a little paper that explains what I did.

Your contributions are appreciated. If you don't have an objection, I may eventually want to use this as a basis generalize the RG filter specification. Right now it only calls out coefficients for 48 and 44.1 kHz implementations.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #30 – 2011-01-16 20:58:26
Quote
now it only calls out coefficients for 48 and 44.1 kHz implementations

wvgain.c from WavPack sources contains coefficients for frequencies from 8 kHz to 192 kHz

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #31 – 2011-01-16 23:18:16
Yesterday I had some fun with Maxima and LaTeX, and wrote a little paper that explains what I did.

Your contributions are appreciated. If you don't have an objection, I may eventually want to use this as a basis generalize the RG filter specification. Right now it only calls out coefficients for 48 and 44.1 kHz implementations.

Please go ahead! If you need the Maxima/tex source code, it's here.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #32 – 2011-01-16 23:25:58
If you want to have my opinion: I'm pretty sure the true peak values can be estimated quite accurately without resampling, so please don't consider adding resampling in the future.

I have now implemented sample peak measurement... Maybe I'll add true peak measurement but make it optional. Scanning speed will remain the primary goal.
And if I read your code correctly, you store every block energy in your linked list, even the ones which are quieter than -70 LUFS (which you don't need later). I assume that's a design limitation (buffering, album gain handling, and stuff)?

Thanks for pointing that out! Now those blocks are not saved anymore.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #33 – 2011-01-16 23:35:17
Quote
now it only calls out coefficients for 48 and 44.1 kHz implementations

wvgain.c from WavPack sources contains coefficients for frequencies from 8 kHz to 192 kHz

Yes, for discrete frequencies, not a compact, closed-form representation of the coefficients as a function of sampling frequency that Raiden has presented here.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #34 – 2011-01-16 23:42:36
0.1.4 has been uploaded, with the following new features:
- r128-sndfile now supports MP1/2/3 via libmpg123. Now r128-ffmpeg should only be needed for formats that libsndfile and libmpg123 can't handle, like AC3 etc. If you use the precompiled Win32 binaries, you'll need the libmpg123 dll from here.
- r128-sndfile now supports measurement of peak values for the ReplayGain tags.
- the library now has a function ebur128_loudness_range, that implements the loudness range measurement algorithm described in EBU - TECH 3342. You can test it with r128-sndfile and the command line option "-r". The test cases already work quite well:
Code: [Select]
`seq-3342-1-16bit.wav - LRA: 10.001105seq-3342-2-16bit.wav - LRA: 4.999373seq-3342-3-16bit.wav - LRA: 19.995064seq-3342-4-16bit.wav - LRA: 14.999274seq-3341-7_seq-3342-5-24bit.wav - LRA: 4.974759seq-3341-8_seq-3342-6-24bit.wav - LRA: 14.993651`

Source (tar.gz)
Source (zip)
Win32 build (zip)

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #35 – 2011-01-17 17:22:13
Yesterday I had some fun with Maxima and LaTeX, and wrote a little paper that explains what I did.

This paper was very helpful. On page 3 there are formulas for b_0, b_1, b_2, a_1, a_2, obtained with bilinear transform.

Many thanks again for providing the links. It's just like pointing me to a goldmine.

As far as I understand, you are calculating
[blockquote]K, Q, Vh, Vb, Vl, from K also Fc[/blockquote]from the given a0...b2.

Am I right in assuming that Q, Vh, Vb, Vl and Fc remain constant for all sample rates and only K has to be adapted accordingly?

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #36 – 2011-01-17 18:11:46
Correct.

V are gain values and have no relation to sampling frequency
Q is a "magic number" that effects the shape of the filter.

Fc stays constant - it's the nominal cutoff frequency.  ω is tan(fc/fs *π)  [it has been typo'd as Ω in the paper,  ω/Ω are lowercase/uppercase pairs].  That is, it's the cutoff frequency as a percentage of the sampling rate, and "pre-warped" with tan() to match the frequency warping done by the bilinear transform.  k is often used for ω in source code.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #37 – 2011-01-17 20:36:35
Thanks, works very well for me.

Just thought I'd mention that the python tagging script did not work out of the box on my system (python 2.6.5), I get the following error message:
Code: [Select]
`Traceback (most recent call last):  File "../src/rgtag.py", line 4, in <module>    tg = '{:.2f}'.format(float(sys.argv[2])) + " dB"ValueError: zero length field name in format`

Substituting with an alternative string formatting, it works like a charm:
Code: [Select]
`tg = "%.2f dB" % float(sys.argv[2])tp = "%.8f"    % float(sys.argv[3])ag = "%.2f dB" % float(sys.argv[4])ap = "%.8f"    % float(sys.argv[5])`

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #38 – 2011-01-18 22:53:37
Substituting with an alternative string formatting, it works like a charm:

Thank you for the fix! This is (literally) my first python script.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #39 – 2011-01-18 23:21:18
Version 0.1.5 has arrived!
- I've seperated the scanners. They all support the same options and differ just in the input formats. There are:
-- r128-sndfile (for Vorbis, FLAC, WAV, and many others)
-- r128-mpg123 (for MP1, MP2, MP3)
-- r128-musepack
-- r128-ffmpeg (everything else)
- Added needed DLLs (FFmpeg is compiled as LGPL version)
- Added Musepack support to the tagging script
- Many improvements in the library (thanks to C.R.Helmrich!)
- Fixed loudness range calculation. In 0.1.4, when a new segment was started, calculation stopped. This is now fixed!
- Minor bug fixes.
- A Win64 build! Please let me know if it works. It's 20-30% faster than the Win32 build. This is probably because SSE optimizations are not enabled in the Win32 build. I will try to do this soon.

Source (tar.gz)
Source (zip)
Win32 build
Win64 build

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #40 – 2011-01-19 15:40:00
Correct.

V are gain values and have no relation to sampling frequency
Q is a "magic number" that effects the shape of the filter.

Fc stays constant - it's the nominal cutoff frequency.  ? is tan(fc/fs *?)  [it has been typo'd as ? in the paper,  ?/? are lowercase/uppercase pairs].  That is, it's the cutoff frequency as a percentage of the sampling rate, and "pre-warped" with tan() to match the frequency warping done by the bilinear transform.  k is often used for ? in source code.

I've finally managed (thanks to the pointers provided by Raiden) to find a closed solution to the re-quantization problem for digital biquad filters as it appears in the context of (but not restricted to) BS.1770.

Assume youv'e given the coefficients a1, a2, b0, b1, b2, b3 of a digital biquad filter for a particular sample frequency Fs (cf. e.g. http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt).

The re-quantization problem consists in calculating the coefficents a1', a2', b0', b1', b2', and b3' of a digital biquad filter with the same characteristics as the given one, but for another sample frequency Fs'.

The key to the solution are the pointers provided by Raiden:
According to [1] the coeffients a1, a2, b0, b1, b2 of the given digital biquad filter can be mapped to the parameters Fc, Q, Vl, Vb, and Vh of the corresponding analog filter by the following equations (using the notation from [2]):

Code: [Select]
`(1)    (1 + K / Q + K^2) * a1  =  2 * (K^2 - 1)(2)    (1 + K / Q + K^2) * a2  =  1 - K / Q + K^2(3)    (1 + K / Q + K^2) * b0  =  Vh + Vb * K / Q + Vl * K^2(4)    (1 + K / Q + K^2) * b1  =  2 * (Vl * K^2 - Vh)(5)    (1 + K / Q + K^2) * b2  =  Vh - Vb * K / Q + Vl * K^2`

with

Code: [Select]
`(6) K = tan(pi * Fc / Fs).`

In order to solve the above stated re-quantization problem we do the following:
• Solve the eqs (1) - (5) for Fc, Q, Vl, Vb, and Vh.
• Use the equations

Code: [Select]
`(1')    (1 + K' / Q + K'^2) * a1'  =  2 * (K'^2 - 1)(2')    (1 + K' / Q + K'^2) * a2'  =  1 - K' / Q + K'^2(3')    (1 + K' / Q + K'^2) * b0'  =  Vh' + Vb' * K' / Q + Vl * K'^2(4')    (1 + K' / Q + K'^2) * b1'  =  2 * (Vl * K'^2 - Vh)(5')    (1 + K' / Q + K'^2) * b2'  =  Vh - Vb * K' / Q + Vl * K^2`

with

Code: [Select]
`(6') K' = tan(pi * Fc / Fs')`

in order to calculate the coefficents a1', a2', b0', b1', and b2' for a digital biquad filter with the same characteristcs Fc, Q, Vl, Vb, and Vh as the given one but for a different sample frequency Fs'.
We start by solving eqs. (1) and (2) for K^2 and K/Q, respectively. By substituting

Code: [Select]
`x11  =  a1 - 2x12  =  a1x1  =  -a1 - 2x21  =  a2 - 1x22  =  a2 + 1x2  =  -a2 + 1DX = x22 * x11 - x12 * x21`

and using well known methods we arrive at

Code: [Select]
`(6)  K^2  =  (x22 * x1 - x12 * x2) / DX(7)  K/Q  =  (x11 * x2 - x21 * x1) / DX`

Next we solve eqs. (3), (4), and (5) for Vh, Vb, and Vl. Introducing

Code: [Select]
`(8) a0  =  1 + K / Q + K^2`

and reordering eqs. (3), (4), and (5) they read

Code: [Select]
`(3a)       Vh + K/Q * Vb +       K^2 * Vl  =  b0 * a0(4a)  -2 * Vh            + (2 * K^2) * Vl  =  b1 * a0(5a)       Vh - K/Q * Vb +       K^2 * Vl  =  b2 * a0`

Now it's not hard any longer to find the solutions:

Code: [Select]
`            a0 * (b0 - b2)(9)  Vb  =  --------------               2 * K/Q            a0 * (b0 + b1 + b2)(10) Vl  =  -------------------                  4 * K^2             a0 * (b0 - b1 + b2)(11) Vh  =  -------------------                    4`

Finally we observe from eqs. (6) and (6') the following:

Code: [Select]
`(6)     K   =  tan(pi * Fc / Fs)(6)'    K'  =  tan(pi * Fc / Fs')(12)   K'  =  tan(atan(K) * Fs / Fs')`

Eqs. (6) and (7) along with (9), (10), (11), and (12) provide everything we need for re-quantizing any given digital biquad filter. We demonstrate this by the following C code (please note that this code re-quantizes digital biqad filters on the fly, no pre-processing by an external algebra system ist needed).

Code: [Select]
`typedef struct biquad {  double fs;  double a1, a2;  double b0, b1, b2;} biquad_t;typedef struct biquad_ps {  double k;  double q;  double vb;  double vl;  double vh;} biquad_ps_t;void biquad_get_ps(biquad_t *biquad, biquad_ps_t *ps){  double x11 = biquad->a1 - 2;  double x12 = biquad->a1;  double x1 = -biquad->a1 - 2;  double x21 = biquad->a2 - 1;  double x22 = biquad->a2 + 1;  double x2 = -biquad->a2 + 1;  double dx = x22*x11 - x12*x21;  double k_sq = (x22*x1 - x12*x2)/dx;  double k_by_q = (x11*x2 - x21*x1)/dx;  double a0 = 1.0 + k_by_q + k_sq;  ps->k = sqrt(k_sq);  ps->q = ps->k/k_by_q;  ps->vb = 0.5*a0*(biquad->b0 - biquad->b2)/k_by_q;  ps->vl = 0.25*a0*(biquad->b0 + biquad->b1 + biquad->b2)/k_sq;  ps->vh = 0.25*a0*(biquad->b0 - biquad->b1 + biquad->b2);}biquad_t *biquad_requantize(biquad_t *in, biquad_t *out){  if (in->fs==out->fs)    return in;  else {    biquad_ps_t ps;    double k, k_sq, k_by_q, a0;    biquad_get_ps(in, &ps);    k=tan((in->fs/out->fs)*atan(ps.k));    k_sq = k*k;    k_by_q = k/ps.q;    a0 = 1.0 + k_by_q + k_sq;    out->a1 = (2.0*(k_sq - 1.0))/a0;    out->a2 = (1.0 - k_by_q + k_sq)/a0;    out->b0 = (ps.vh + ps.vb*k_by_q + ps.vl*k_sq)/a0;    out->b1 = (2.0 * (ps.vl*k_sq - ps.vh))/a0;    out->b2 = (ps.vh - ps.vb*k_by_q + ps.vl*k_sq)/a0;    return out;  }}`

The following code demonstrates how to re-quantize the 48 kHz BS.1770 pre-filter to it's 44.1 kHz equivalent using the above functions.

Code: [Select]
`int main(int argc, char **argv){  int i;  biquad_t pre48000={    .fs=48000,    .a1=-1.69065929318241,    .a2=0.73248077421585,    .b0=1.53512485958697,    .b1=-2.69169618940638,    .b2=1.19839281085285  };  biquad_t pre44100={ .fs=44100 };  biquad_requantize(&pre48000, &pre44100);  printf("a1: %f, %f\n", pre48000.a1, pre44100.a1);  printf("a2: %f, %f\n", pre48000.a2, pre44100.a2);  printf("b1: %f, %f\n", pre48000.b0, pre44100.b0);  printf("b2: %f, %f\n", pre48000.b1, pre44100.b1);  printf("b3: %f, %f\n", pre48000.b2, pre44100.b2);}`

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #41 – 2011-01-23 15:49:29
The following code demonstrates how to re-quantize the 48 kHz BS.1770 pre-filter to it's 44.1 kHz equivalent using the above functions.

I get exactly the same coefficients. Nice work!

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #42 – 2011-01-23 16:00:12
- The main new feature is a scanning/tagging script written in Python. It takes one or more directory(ies) as input and scans them for music files. If there are many directories, it will scan them in parallel to use all available cores:
Code: [Select]
`usage: r128-tag [-r] <directory(ies)> ...       r128-tag <file(s)> ...    -h: show this help    -r: scan directory(ies) recursively  r128-tag scans music files and tags them with ReplayGain compatible tags.  If more than one directory is given, all available cores will be used.`

- optimized the LRA calculations
- some minor library improvements
- Linux builds!

Source (tar.gz)
Source (zip)
Win32 build (zip)
Win64 build (zip)
Linux32 build (tar.gz)
Linux64 build (tar.gz)

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #43 – 2011-01-23 17:17:21
Any idea why r128-mpg123 does not work for me on Windows 7 x64?
Code: [Select]
`C:\Users\frost\Downloads\libebur128-0.1.6-win64>python ./r128-tag.py -r "C:\Users\frost\Downloads\Motörhead\2008 - Motörizer"Note: Illegal Audio-MPEG-Header 0x41504554 at offset 5074432.Note: Trying to resync...Note: Skipped 1024 bytes in input.Internal MPG123 error!Note: Illegal Audio-MPEG-Header 0x41504554 at offset 5106947.Note: Trying to resync...Note: Skipped 1024 bytes in input.Internal MPG123 error!Note: Illegal Audio-MPEG-Header 0x41504554 at offset 5802950.Note: Trying to resync...Note: Skipped 1024 bytes in input.Internal MPG123 error!`

These are normal MP3 files encoded with Lame 3.93. If i strip all tags from the file it doesn't complain anymore. Files are tagged with ID3 and APEv2.
Blubb

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #44 – 2011-01-23 19:00:36
Thanks for the new version. When I check the same file twice by calling r128-sndfile {-t} file.wav file.wav, I sometimes get slightly different results (by 0.01 LU or so). Could this indicate an initialization bug? The same happens on version 0.1.5.

Btw, did you manage to get the Win32 build SSE optimized?

Chris

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #45 – 2011-01-23 19:28:27
Thanks for the new version. When I check the same file twice by calling r128-sndfile {-t} file.wav file.wav, I sometimes get slightly different results (by 0.01 LU or so). Could this indicate an initialization bug? The same happens on version 0.1.5.

This is intentional (for now). When starting a new segment, the last unfinished block of the last segment does not get dropped, but becomes the first block of the new segment. Only the last unfinished block of the whole programme gets dropped. I wonder how other implementations and ReplayGain scanners handle this edge-case.

Btw, did you manage to get the Win32 build SSE optimized?

Not yet... But it should just be "-msse2 -mfpmath=sse" so I will do it in the next version.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #46 – 2011-01-23 19:42:05
Any idea why r128-mpg123 does not work for me on Windows 7 x64?
These are normal MP3 files encoded with Lame 3.93. If i strip all tags from the file it doesn't complain anymore. Files are tagged with ID3 and APEv2.

Thanks for the report! I could reproduce with APEv2 tagged MP3's. It looks like MPG123 sees APE tags as errors in the stream (which they are in some way, I guess ). I will try to handle this more gracefully.

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #48 – 2011-01-24 10:55:20
Thanks. Could you fix the build script to run on FreeBSD? It has troubles finding libraries. I got r128-mpg123 to compile by installing mpg123 and making a symlink:
Code: [Select]
`ln -s /usr/local/lib/libmpg123.so /lib/libmpg123.so`

Same for r128-sndfile:
Code: [Select]
`ln -s /usr/local/lib/libsndfile.so /lib/libsndfile.so`

And for r128-ffmpeg you have to create 2 symlinks:
Code: [Select]
`ln -s /usr/local/lib/libavcodec.so /lib/libavcodec.soln -s /usr/local/lib/libavformat.so /lib/libavformat.so`

After this the software is running fine on FreeBSD, unfortunately, I couldn't use the python script:
Code: [Select]
`Traceback (most recent call last):  File "./r128-tag.py", line 112, in <module>    pool = multiprocessing.Pool(processes=number_threads)  File "/usr/local/lib/python2.6/multiprocessing/__init__.py", line 227, in Pool    return Pool(processes, initializer, initargs)  File "/usr/local/lib/python2.6/multiprocessing/pool.py", line 84, in __init__    self._setup_queues()  File "/usr/local/lib/python2.6/multiprocessing/pool.py", line 130, in _setup_queues    from .queues import SimpleQueue  File "/usr/local/lib/python2.6/multiprocessing/queues.py", line 22, in <module>    from multiprocessing.synchronize import Lock, BoundedSemaphore, Semaphore, Condition  File "/usr/local/lib/python2.6/multiprocessing/synchronize.py", line 33, in <module>    " function, see issue 3770.")ImportError: This platform lacks a functioning sem_open implementation, therefore, the required synchronization primitives needed will not function, see issue 3770.`

Solution is to compile lang/python27 or whatever version you prefer and enable SEM and PTH options. Then forget to install mutagen and uncomment try-expect block to learn about it, install mutagen and voilá: Runs on FreeBSD .
Blubb

## libebur128 - (yet another) EBU R 128 implementation

##### Reply #49 – 2011-01-24 11:45:15
Thanks. Could you fix the build script to run on FreeBSD? It has troubles finding libraries. I got r128-mpg123 to compile by installing mpg123 and making a symlink:

Yay! Another BSD user  Development so far is done on Linux, but I will try to build on FreeBSD later today.
After this the software is running fine on FreeBSD, unfortunately, I couldn't use the python script:
Solution is to compile lang/python27 or whatever version you prefer and enable SEM and PTH options. Then forget to install mutagen and uncomment try-expect block to learn about it, install mutagen and voilá: Runs on FreeBSD .

Awesome. I've already made the multiprocessing module optional and improved the tagging error messages.