+++ /dev/null
-prefix=@prefix@
-libdir=@libdir@
-includedir=@includedir@
-
-Name: libleqm-nrt
-Description: LEQ calculation library
-Version: @version@
-Libs: @libs@
-Cflags: -I${includedir}
--- /dev/null
+prefix=@prefix@
+libdir=@libdir@
+includedir=@includedir@
+
+Name: libleqm_nrt
+Description: LEQ calculation library
+Version: @version@
+Libs: @libs@
+Cflags: -I${includedir}
+++ /dev/null
-/*
- leqm-nrt is a non-real-time implementation
- of Leq(M) measurement according to ISO 21727:2004(E)
- "Cinematography -- Method of measurement of perceived
- loudness of motion-picture audio material"
-
- Copyright (C) 2011-2013, 2017-2018 Luca Trisciani
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
- */
-
-#include "leqm-nrt.h"
-
-#include <stdio.h>
-#include <math.h>
-#include <sndfile.h>
-#include <unistd.h>
-#include <pthread.h>
-#include <string.h>
-#include <stdlib.h>
-#include <time.h>
-#include <ctype.h>
-#include <iso646.h>
-
-#include <vector>
-#include <memory>
-#include <iostream>
-#include <thread>
-#include <mutex>
-#include <functional>
-
-// Version 0.0.18 (C) Luca Trisciani 2011-2013, 2017-2018
-// Tool from the DCP-Werkstatt Software Bundle
-
-
-
-// COMPILATION
-// compile for DEBUG with gcc -g -DEBUG -lsndfile -lfftw3 -lm -lpthread -lrt -o leqm-nrt leqm-nrt.cpp
-//otherwise gcc -lsndfile -lm -lpthread -lrt -o leqm-nrt leqm-nrt.c
-
-//#define DEBUG
-
-int main(int argc, const char ** argv)
-{
- int number_of_filter_interpolation_points = 64; // This value is low for precision. Calibration is done with 32768 point.
- int num_cpu = std::thread::hardware_concurrency() - 1;
-
- printf("leqm-nrt Copyright (C) 2011-2013, 2017-2018 Luca Trisciani\nThis program comes with ABSOLUTELY NO WARRANTY.\nThis is free software, and you are welcome to redistribute it\nunder the GPL v3 licence.\nProgram will use 1 + %d slave threads.\n", num_cpu);
- int buffer_size_ms = 850; //ISO 21727:2004 do not contain any indication, TASA seems to indicate 1000, p. 8
- std::vector<double> channel_corrections;
- int parameterstate = 0;
- bool display_leqnw = false;
- std::string sound_filename;
-
- if (argc == 1)
- { const char helptext[] = "Order of parameters is free.\nPossible parameters are:\n-convpoints <integer number> \tNumber of interpolation points for the filter. Default 64\n-numcpus <integer number> \tNumber of slave threads to speed up operation.\n-timing \t\t\tFor benchmarking speed.\n-leqnw\t Print out Leq without M Weighting\n-chconfcal <db correction> <db correction> <etc. so many times as channels>\n-logleqm10\n-logleqm\n-buffersize <milliseconds>\n";
- printf(helptext);
- printf("Please indicate a sound file to be processed.\n");
- return 0;
- }
-
- for (int in = 1; in < argc;) {
-
- if (!(strncmp(argv[in], "-", 1) == 0)) {
- sound_filename = argv[in];
- in++;
- }
-
- if (strcmp(argv[in], "-chconfcal") == 0) {
- /* as the order of parameter is free I have to postpone
- the check for consistency with the number of channels.
- So first create a temporary array, whose number of element will be checked after
- the parsing of the command line parameters is finished.
- The calibration will be expressed in dB on the command line and converted to multiplier
- here so that it can be stored as a factor in the channelconfcalvector.
- */
-
- in++;
- for (;;) {
- if (in < argc) {
- if (!(strncmp(argv[in], "-", 1) == 0) || isdigit(argv[in][1])) {
- channel_corrections.push_back(leqm_nrt::convert_log_to_linear_single(atof(argv[in++])));
- } else break;
- } else break;
-
- } //for
- continue;
- }
-
- if (strcmp(argv[in], "-convpoints") == 0) {
- number_of_filter_interpolation_points = atoi(argv[in + 1]);
- in+=2;
- printf("Convolution points sets to %d.\n", number_of_filter_interpolation_points);
- continue;
-
- }
-
-
- if (strcmp(argv[in], "-version") == 0) {
- in++;
- printf("leqm-nrt version 0.18\n");
- continue;
-
- }
- if (strcmp(argv[in], "-numcpus") == 0) {
- num_cpu= atoi(argv[in + 1]);
- in+=2;
- printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", num_cpu);
- continue;
-
- }
-
- if (strcmp(argv[in], "-leqnw") == 0) {
- display_leqnw = true;
- in++;
- printf("Leq(nW) - unweighted - will be outputted.\n");
- continue;
-
- }
-
- if (strcmp(argv[in], "-buffersize") == 0) {
- buffer_size_ms = atoi(argv[in + 1]);
- in+=2;
- printf("Buffersize will be set to %d milliseconds.\n", buffer_size_ms);
- continue;
-
- }
-
- if (parameterstate==0) {
- break;
- }
- }
-
- auto result = leqm_nrt::calculate_file(sound_filename, channel_corrections, buffer_size_ms, number_of_filter_interpolation_points, num_cpu);
-
- if (display_leqnw) {
- printf("Leq(nW): %.4f\n", result.leq_nw); // Leq(no Weighting)
- }
- printf("Leq(M): %.4f\n", result.leq_m);
-
- return result.status;
-}
-
+++ /dev/null
-/*
- leqm-nrt is a non-real-time implementation
- of Leq(M) measurement according to ISO 21727:2004(E)
- "Cinematography -- Method of measurement of perceived
- loudness of motion-picture audio material"
-
- Copyright (C) 2011-2013, 2017-2018 Luca Trisciani
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
- */
-
-#include "leqm-nrt.h"
-
-#include <stdio.h>
-#include <math.h>
-#include <sndfile.h>
-#include <unistd.h>
-#include <pthread.h>
-#include <string.h>
-#include <stdlib.h>
-#include <time.h>
-#include <ctype.h>
-#include <iso646.h>
-
-#include <vector>
-#include <memory>
-#include <iostream>
-#include <thread>
-#include <mutex>
-#include <functional>
-
-using namespace leqm_nrt;
-
-// Version 0.0.18 (C) Luca Trisciani 2011-2013, 2017-2018
-// Tool from the DCP-Werkstatt Software Bundle
-
-namespace leqm_nrt {
-
-class Worker
-{
-public:
- Worker(std::vector<double> buffer, int nsamples, int nch, int npoints, std::vector<double> const& ir, Sum* sum, std::vector<double> chconf)
- : _buffer(buffer)
- , _nsamples(nsamples)
- , _nch(nch)
- , _npoints(npoints)
- , _ir(ir)
- , _sum(sum)
- , _chconf(chconf)
- {
- _thread = std::thread(&Worker::process, this);
- }
-
- Worker(Worker& other) = delete;
- bool operator=(Worker&) = delete;
- Worker(Worker&& other) = delete;
- bool operator=(Worker&&) = delete;
-
- ~Worker()
- {
- try {
- _thread.join();
- } catch (...)
- {
- }
- }
-
-
-private:
- double sum_and_short_term_avrg(std::vector<double>const & channel_accumulator, int nsamples) const
- {
- double stsum = 0.0;
- for (auto i = 0; i < nsamples; i++) {
- stsum += channel_accumulator[i];
-
- }
- return stsum / nsamples;
- }
-
- void accumulate_ch(std::vector<double>& ch_accumulator, std::vector<double> const& input_channel, int nsamples) const
- {
- for (auto i = 0; i < nsamples; i++) {
- ch_accumulator[i] += input_channel[i];
- }
- }
-
- std::vector<double> rectify(std::vector<double> const& input) const
- {
- std::vector<double> squared;
- for (auto i: input) {
- squared.push_back(pow(i, 2));
- }
- return squared;
- }
-
- std::vector<double> convolve(std::vector<double> const& signal, std::vector<double> const& ir) const
- {
- std::vector<double> result(signal.size());
- double sum = 0.0;
- for (auto i = 0U; i < signal.size(); i++) {
- auto m = i;
- for (int l = ir.size()- 1; l >= 0; l--, m++) {
- if (m >= signal.size()) {
- m -= signal.size();
- }
- sum += signal[m] * ir[l];
- }
- result[i] = sum;
- sum = 0.0;
- }
- return result;
- }
-
- void process()
- {
- int const frames = _nsamples / _nch;
-
- std::vector<double> ch_sum_accumulator_norm(frames);
- std::vector<double> ch_sum_accumulator_conv(frames);
-
- for (int ch = 0; ch < _nch; ch++) {
-
- std::vector<double> normalized_buffer(frames);
-
- for (int n = ch, m = 0; n < _nsamples; n += _nch, m++) {
- // use this for calibration depending on channel config for ex. chconf[6] = {1.0, 1.0, 1.0, 1.0, 0.707945784, 0.707945784} could be the default for 5.1 soundtracks
- //so not normalized but calibrated
- normalized_buffer[m] = _buffer[n] * _chconf[ch]; //this scale amplitude according to specified calibration
- }
-
- auto convolved_buffer = convolve(normalized_buffer, _ir);
- accumulate_ch(ch_sum_accumulator_norm, rectify(normalized_buffer), frames);
- accumulate_ch(ch_sum_accumulator_conv, rectify(convolved_buffer), frames);
- }
-
- _sum->sum_samples(ch_sum_accumulator_norm, ch_sum_accumulator_conv, frames);
- }
-
- std::vector<double> _buffer;
- int _nsamples;
- int _nch;
- int _npoints;
- std::vector<double> const& _ir;
- Sum* _sum;
- std::vector<double> _chconf;
-
- std::thread _thread;
-};
-
-}
-
-
-//the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
-//it is also different for the way it interpolates between DC and 31 Hz
-// Pay attention that also arguments to the functions are changed
-static std::vector<double> equalinterval2(int points, int samplingfreq, int bitdepthsoundfile)
-{
- //ISO 21727:2004(E)
- // M Weighting
- double const frequencies[] = {31, 63, 100, 200, 400, 800, 1000, 2000, 3150, 4000, 5000, 6300, 7100, 8000, 9000, 10000, 12500, 14000, 16000, 20000, 31500};
- double const db[] = {-35.5, -29.5, -25.4, -19.4, -13.4, -7.5, -5.6, 0.0, 3.4, 4.9, 6.1, 6.6, 6.4, 5.8, 4.5, 2.5, -5.6, -10.9, -17.3, -27.8, -48.3};
- int constexpr points_in_standard_ccir_filter = 21;
-
- std::vector<double> freq_response(points);
- //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
- double const dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
- //double dcatt = -90.3;
- double const pass = ((double) (samplingfreq >> 1)) / ((double) points);
- for (int ieq = 0, i = 0; ieq < points; ieq++) {
- double const freq = ieq * pass;
-
- if (freq == 0.0) {
- freq_response[ieq] = dcatt;
- } else if (freq < frequencies[0]) { // this has a lot of influence on final Leq(M) value
- freq_response[ieq] = ((db[0] - dcatt) / (frequencies[0] - 0)) * freq + dcatt;
- continue;
- } else {
-
- if (freq >= frequencies[i] && freq < frequencies[i+1]) {
- freq_response[ieq] = ( (db[i+1] - db[i]) / (frequencies[i+1] - frequencies[i])) * (freq - frequencies[i]) + db[i];
- } else if (freq >= frequencies[i+1]) {
- while (freq >= frequencies[i+1]) {
- i++;
- if ((i + 1) >= points_in_standard_ccir_filter) {
- break;
- }
- }
- if ((i+1) < points_in_standard_ccir_filter) {
- freq_response[ieq] = ( (db[i+1] - db[i]) / (frequencies[i+1] - frequencies[i])) * (freq - frequencies[i]) + db[i];
- } else {
- freq_response[ieq] = ( (1 - db[i]) / ((samplingfreq >> 1) - frequencies[i])) * (freq - frequencies[i]) + db[i];
- }
- }
- }
- }
- return freq_response;
-}
-
-
-static std::vector<double> convert_log_to_linear(std::vector<double> const& in)
-{
- std::vector<double> out(in.size());
- for (auto i = 0U; i < in.size(); i++) {
- out[i] = powf(10, in[i] / 20.0);
- }
- return out;
-}
-
-
-static std::vector<double> inverse_fft(std::vector<double> const& freq_response)
-{
- int const npoints = freq_response.size();
- std::vector<double> ir(npoints * 2);
- for (int n = 0; n < npoints; n++) {
- double parsum = 0.0;
- double partial = 0.0;
-
- for (int m = 1; m <= npoints - 1; m++) {
- partial = cos(2.0 * M_PI * m * ((n - (npoints * 2.0 - 1) / 2) / (npoints * 2.0)));
- parsum = parsum + freq_response[m] * partial;
- }
- ir[n] = (freq_response[0] + 2.0 * parsum) / (npoints * 2.0);
- }
-
- for (int n = 0; n < npoints; n++) {
- ir[npoints + n] = ir[npoints - (n + 1)];
- }
-
- return ir;
-}
-
-
-Result leqm_nrt::calculate_file(
- std::string sound_filename,
- std::vector<double> channel_corrections,
- int buffer_size_ms,
- int number_of_filter_interpolation_points,
- int num_cpu
- )
-{
- SF_INFO sf_info;
- sf_info.format = 0;
- SNDFILE* file = sf_open(sound_filename.c_str(), SFM_READ, &sf_info);
- if (!file) {
- return {-101};
- }
-
- int bitdepth = 0;
- switch(sf_info.format & SF_FORMAT_SUBMASK) {
- // all signed bitdepth
- case 0x0001:
- bitdepth = 8;
- break;
- case 0x0002:
- bitdepth = 16;
- break;
- case 0x0003:
- bitdepth = 24;
- break;
- case 0x0004:
- bitdepth = 32;
- break;
- default:
- printf("No known bitdepth! Exiting ...\n");
- return -1;
- }
-
- Calculator calculator(
- sf_info.channels,
- sf_info.samplerate,
- bitdepth,
- channel_corrections,
- buffer_size_ms,
- number_of_filter_interpolation_points,
- num_cpu
- );
-
- while (true) {
- std::vector<double> buffer(4096);
- auto read = sf_read_double(file, buffer.data(), buffer.size());
- if (read <= 0) {
- break;
- }
-
- buffer.resize(read);
- calculator.add(buffer);
- }
-
- sf_close(file);
-
- return {calculator.leq_m(), calculator.leq_nw()};
-}
-
-
-std::vector<double> calculate_ir(double number_of_filter_interpolation_points, int sample_rate, int bits_per_sample)
-{
- return inverse_fft(convert_log_to_linear(equalinterval2(number_of_filter_interpolation_points, sample_rate, bits_per_sample)));
-}
-
-
-/** @return List of default channel corrections for the given number of channels,
- * or an empty vector if we can't offer a default.
- */
-std::vector<double> default_channel_corrections(int channels)
-{
- std::vector<double> corr;
-
- if (channels == 6) {
- double conf51[] = {0, 0, 0, 0, -3, -3};
- for (auto cind = 0; cind < channels; cind++) {
- corr.push_back(convert_log_to_linear_single(conf51[cind]));
- }
- }
-
- return corr;
-}
-
-
-double leqm_nrt::convert_log_to_linear_single(double in)
-{
- return powf(10, in / 20.0f);
-}
-
-
-Calculator::Calculator(
- int channels,
- int sample_rate,
- int bits_per_sample,
- std::vector<double> channel_corrections,
- int buffer_size_ms,
- int number_of_filter_interpolation_points,
- int num_cpu
- )
- : _channels(channels)
- , _channel_corrections(channel_corrections)
- , _number_of_filter_interpolation_points(number_of_filter_interpolation_points)
- , _num_cpu(num_cpu)
-{
- if ((sample_rate * buffer_size_ms) % 1000) {
- throw BadBufferSizeError();
- }
-
- if (static_cast<int>(channel_corrections.size()) != channels) {
- channel_corrections = default_channel_corrections(channels);
- }
- if (static_cast<int>(channel_corrections.size()) != channels) {
- throw BadChannelCorrectionsError();
- }
-
- _ir = calculate_ir(number_of_filter_interpolation_points, sample_rate, bits_per_sample);
- _buffer.resize((sample_rate * channels * buffer_size_ms) / 1000);
-}
-
-
-void Calculator::process_buffer()
-{
- if (_buffer_free_offset == 0) {
- return;
- }
-
- _workers.push_back(
- std::make_shared<Worker>(
- _buffer,
- _buffer_free_offset,
- _channels,
- _number_of_filter_interpolation_points,
- _ir,
- &_sum,
- _channel_corrections
- )
- );
-
- if (static_cast<int>(_workers.size()) == _num_cpu) {
- _workers.clear();
- }
-
- _buffer_free_offset = 0;
-}
-
-
-void Calculator::add(std::vector<double> samples)
-{
- size_t samples_offset = 0;
-
- while (samples_offset < samples.size()) {
- /* Copy some data into our buffer */
- auto to_copy = std::min(samples.size() - samples_offset, _buffer.size() - _buffer_free_offset);
- memcpy(_buffer.data() + _buffer_free_offset, samples.data() + samples_offset, to_copy * sizeof(double));
- samples_offset += to_copy;
- _buffer_free_offset += to_copy;
-
- /* Process the buffer if it's full */
- if (_buffer_free_offset == _buffer.size()) {
- process_buffer();
- }
- }
-}
-
-
-double Calculator::leq_m()
-{
- process_buffer();
- _workers.clear();
-
- return _sum.leqm();
-}
-
-
-double Calculator::leq_nw()
-{
- process_buffer();
- _workers.clear();
-
- return _sum.rms();
-}
+++ /dev/null
-#include <string>
-#include <vector>
-#include <functional>
-#include <stdexcept>
-#include <mutex>
-#include <cmath>
-#include <memory>
-
-namespace leqm_nrt {
-
-class Sum
-{
-public:
- void sum_samples(std::vector<double> const& input_samples, std::vector<double> const& c_input_samples, int nsamples)
- {
- _mutex.lock();
- _nsamples += nsamples;
- for (auto i = 0; i < nsamples; i++) {
- _sum += input_samples[i];
- _csum += c_input_samples[i];
- }
- _mutex.unlock();
- }
-
- int nsamples() const
- {
- return _nsamples;
- }
-
- /*
- How the final offset is calculated without reference to a test tone:
- P0 is the SPL reference 20 uPa
- Reference SPL is RMS ! So 85 SPL over 20 uPa is 10^4.25 x 0.000020 = 0.355655882 Pa (RMS),
- but Peak value is 0.355655882 x sqr(2) = 0.502973372 that is 20 x log ( 0.502973372 / 0.000020) = 88.010299957
- To that one has to add the 20 dB offset of the reference -20dBFS: 88.010299957 + 20.00 = 108.010299957
-
- But ISO 21727:2004(E) ask for a reference level "measured using an average responding meter". So reference level is not 0.707, but 0.637 = 2/pi
- */
-
- double rms() const
- {
- return 20 * log10(mean()) + 108.010299957;
- }
-
- double leqm() const
- {
- return 20 * log10(cmean()) + 108.010299957;
- }
-
-private:
- double mean() const
- {
- return pow(_sum / _nsamples, 0.500);
- }
-
- double cmean() const
- {
- return pow(_csum / _nsamples, 0.500);
- }
-
- double _csum = 0.0; // convolved sum
- double _sum = 0.0; // flat sum
- int _nsamples = 0;
- std::mutex _mutex;
-};
-
-
-struct Result
-{
- Result(int status_)
- : status(status_)
- {}
-
- Result(double leq_m_, double leq_nw_)
- : status(0)
- , leq_m(leq_m_)
- , leq_nw(leq_nw_)
- {}
-
- /** 0 on success, or
- *
- * -100: Either channel_corrections contained a different number of
- * calibrations than; number of channels in the file or it was empty and the
- * program cannot infer one from the number of channels. Please specify a
- * values in channel_corrections.
- *
- * -101: Failed to open the sound file.
- *
- * -102: buffer_size_ms is not an integer number of samples at the sound file's rate.
- */
- int status;
-
- double leq_m;
- double leq_nw;
-};
-
-
-Result calculate_file(
- std::string sound_filename,
- std::vector<double> channel_corrections,
- int buffer_size_ms,
- int number_of_filter_interpolation_points,
- int num_cpu
- );
-
-
-double convert_log_to_linear_single(double in);
-
-
-class BadBufferSizeError : public std::runtime_error
-{
-public:
- BadBufferSizeError()
- : std::runtime_error("Buffer size does not correspond to an integer number of samples")
- {}
-};
-
-
-class BadChannelCorrectionsError : public std::runtime_error
-{
-public:
- BadChannelCorrectionsError()
- : std::runtime_error("Incorrect number of channel corrections given, and no defaults are available")
- {}
-};
-
-
-class Worker;
-
-
-class Calculator
-{
-public:
- Calculator(
- int channels,
- int sample_rate,
- int bits_per_sample,
- std::vector<double> channel_corrections,
- int buffer_size_ms,
- int number_of_filter_interpolation_points,
- int num_cpu
- );
-
- Calculator(Calculator&) = delete;
- Calculator(Calculator&&) = delete;
- bool operator=(Calculator&) = delete;
- bool operator=(Calculator&&) = delete;
-
- void add(std::vector<double> samples);
-
- double leq_m();
- double leq_nw();
-
-private:
- void process_buffer();
-
- int _channels;
- std::vector<double> _channel_corrections;
- int _number_of_filter_interpolation_points;
- int _num_cpu;
- std::vector<std::shared_ptr<Worker>> _workers;
- Sum _sum;
- std::vector<double> _ir;
- std::vector<double> _buffer;
- size_t _buffer_free_offset = 0;
-};
-
-}
-
--- /dev/null
+/*
+ leqm-nrt is a non-real-time implementation
+ of Leq(M) measurement according to ISO 21727:2004(E)
+ "Cinematography -- Method of measurement of perceived
+ loudness of motion-picture audio material"
+
+ Copyright (C) 2011-2013, 2017-2018 Luca Trisciani
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ */
+
+#include "leqm-nrt.h"
+
+#include <stdio.h>
+#include <math.h>
+#include <sndfile.h>
+#include <unistd.h>
+#include <pthread.h>
+#include <string.h>
+#include <stdlib.h>
+#include <time.h>
+#include <ctype.h>
+#include <iso646.h>
+
+#include <vector>
+#include <memory>
+#include <iostream>
+#include <thread>
+#include <mutex>
+#include <functional>
+
+using namespace leqm_nrt;
+
+// Version 0.0.18 (C) Luca Trisciani 2011-2013, 2017-2018
+// Tool from the DCP-Werkstatt Software Bundle
+
+namespace leqm_nrt {
+
+class Worker
+{
+public:
+ Worker(std::vector<double> buffer, int nsamples, int nch, int npoints, std::vector<double> const& ir, Sum* sum, std::vector<double> chconf)
+ : _buffer(buffer)
+ , _nsamples(nsamples)
+ , _nch(nch)
+ , _npoints(npoints)
+ , _ir(ir)
+ , _sum(sum)
+ , _chconf(chconf)
+ {
+ _thread = std::thread(&Worker::process, this);
+ }
+
+ Worker(Worker& other) = delete;
+ bool operator=(Worker&) = delete;
+ Worker(Worker&& other) = delete;
+ bool operator=(Worker&&) = delete;
+
+ ~Worker()
+ {
+ try {
+ _thread.join();
+ } catch (...)
+ {
+ }
+ }
+
+
+private:
+ double sum_and_short_term_avrg(std::vector<double>const & channel_accumulator, int nsamples) const
+ {
+ double stsum = 0.0;
+ for (auto i = 0; i < nsamples; i++) {
+ stsum += channel_accumulator[i];
+
+ }
+ return stsum / nsamples;
+ }
+
+ void accumulate_ch(std::vector<double>& ch_accumulator, std::vector<double> const& input_channel, int nsamples) const
+ {
+ for (auto i = 0; i < nsamples; i++) {
+ ch_accumulator[i] += input_channel[i];
+ }
+ }
+
+ std::vector<double> rectify(std::vector<double> const& input) const
+ {
+ std::vector<double> squared;
+ for (auto i: input) {
+ squared.push_back(pow(i, 2));
+ }
+ return squared;
+ }
+
+ std::vector<double> convolve(std::vector<double> const& signal, std::vector<double> const& ir) const
+ {
+ std::vector<double> result(signal.size());
+ double sum = 0.0;
+ for (auto i = 0U; i < signal.size(); i++) {
+ auto m = i;
+ for (int l = ir.size()- 1; l >= 0; l--, m++) {
+ if (m >= signal.size()) {
+ m -= signal.size();
+ }
+ sum += signal[m] * ir[l];
+ }
+ result[i] = sum;
+ sum = 0.0;
+ }
+ return result;
+ }
+
+ void process()
+ {
+ int const frames = _nsamples / _nch;
+
+ std::vector<double> ch_sum_accumulator_norm(frames);
+ std::vector<double> ch_sum_accumulator_conv(frames);
+
+ for (int ch = 0; ch < _nch; ch++) {
+
+ std::vector<double> normalized_buffer(frames);
+
+ for (int n = ch, m = 0; n < _nsamples; n += _nch, m++) {
+ // use this for calibration depending on channel config for ex. chconf[6] = {1.0, 1.0, 1.0, 1.0, 0.707945784, 0.707945784} could be the default for 5.1 soundtracks
+ //so not normalized but calibrated
+ normalized_buffer[m] = _buffer[n] * _chconf[ch]; //this scale amplitude according to specified calibration
+ }
+
+ auto convolved_buffer = convolve(normalized_buffer, _ir);
+ accumulate_ch(ch_sum_accumulator_norm, rectify(normalized_buffer), frames);
+ accumulate_ch(ch_sum_accumulator_conv, rectify(convolved_buffer), frames);
+ }
+
+ _sum->sum_samples(ch_sum_accumulator_norm, ch_sum_accumulator_conv, frames);
+ }
+
+ std::vector<double> _buffer;
+ int _nsamples;
+ int _nch;
+ int _npoints;
+ std::vector<double> const& _ir;
+ Sum* _sum;
+ std::vector<double> _chconf;
+
+ std::thread _thread;
+};
+
+}
+
+
+//the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
+//it is also different for the way it interpolates between DC and 31 Hz
+// Pay attention that also arguments to the functions are changed
+static std::vector<double> equalinterval2(int points, int samplingfreq, int bitdepthsoundfile)
+{
+ //ISO 21727:2004(E)
+ // M Weighting
+ double const frequencies[] = {31, 63, 100, 200, 400, 800, 1000, 2000, 3150, 4000, 5000, 6300, 7100, 8000, 9000, 10000, 12500, 14000, 16000, 20000, 31500};
+ double const db[] = {-35.5, -29.5, -25.4, -19.4, -13.4, -7.5, -5.6, 0.0, 3.4, 4.9, 6.1, 6.6, 6.4, 5.8, 4.5, 2.5, -5.6, -10.9, -17.3, -27.8, -48.3};
+ int constexpr points_in_standard_ccir_filter = 21;
+
+ std::vector<double> freq_response(points);
+ //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
+ double const dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
+ //double dcatt = -90.3;
+ double const pass = ((double) (samplingfreq >> 1)) / ((double) points);
+ for (int ieq = 0, i = 0; ieq < points; ieq++) {
+ double const freq = ieq * pass;
+
+ if (freq == 0.0) {
+ freq_response[ieq] = dcatt;
+ } else if (freq < frequencies[0]) { // this has a lot of influence on final Leq(M) value
+ freq_response[ieq] = ((db[0] - dcatt) / (frequencies[0] - 0)) * freq + dcatt;
+ continue;
+ } else {
+
+ if (freq >= frequencies[i] && freq < frequencies[i+1]) {
+ freq_response[ieq] = ( (db[i+1] - db[i]) / (frequencies[i+1] - frequencies[i])) * (freq - frequencies[i]) + db[i];
+ } else if (freq >= frequencies[i+1]) {
+ while (freq >= frequencies[i+1]) {
+ i++;
+ if ((i + 1) >= points_in_standard_ccir_filter) {
+ break;
+ }
+ }
+ if ((i+1) < points_in_standard_ccir_filter) {
+ freq_response[ieq] = ( (db[i+1] - db[i]) / (frequencies[i+1] - frequencies[i])) * (freq - frequencies[i]) + db[i];
+ } else {
+ freq_response[ieq] = ( (1 - db[i]) / ((samplingfreq >> 1) - frequencies[i])) * (freq - frequencies[i]) + db[i];
+ }
+ }
+ }
+ }
+ return freq_response;
+}
+
+
+static std::vector<double> convert_log_to_linear(std::vector<double> const& in)
+{
+ std::vector<double> out(in.size());
+ for (auto i = 0U; i < in.size(); i++) {
+ out[i] = powf(10, in[i] / 20.0);
+ }
+ return out;
+}
+
+
+static std::vector<double> inverse_fft(std::vector<double> const& freq_response)
+{
+ int const npoints = freq_response.size();
+ std::vector<double> ir(npoints * 2);
+ for (int n = 0; n < npoints; n++) {
+ double parsum = 0.0;
+ double partial = 0.0;
+
+ for (int m = 1; m <= npoints - 1; m++) {
+ partial = cos(2.0 * M_PI * m * ((n - (npoints * 2.0 - 1) / 2) / (npoints * 2.0)));
+ parsum = parsum + freq_response[m] * partial;
+ }
+ ir[n] = (freq_response[0] + 2.0 * parsum) / (npoints * 2.0);
+ }
+
+ for (int n = 0; n < npoints; n++) {
+ ir[npoints + n] = ir[npoints - (n + 1)];
+ }
+
+ return ir;
+}
+
+
+Result leqm_nrt::calculate_file(
+ std::string sound_filename,
+ std::vector<double> channel_corrections,
+ int buffer_size_ms,
+ int number_of_filter_interpolation_points,
+ int num_cpu
+ )
+{
+ SF_INFO sf_info;
+ sf_info.format = 0;
+ SNDFILE* file = sf_open(sound_filename.c_str(), SFM_READ, &sf_info);
+ if (!file) {
+ return {-101};
+ }
+
+ int bitdepth = 0;
+ switch(sf_info.format & SF_FORMAT_SUBMASK) {
+ // all signed bitdepth
+ case 0x0001:
+ bitdepth = 8;
+ break;
+ case 0x0002:
+ bitdepth = 16;
+ break;
+ case 0x0003:
+ bitdepth = 24;
+ break;
+ case 0x0004:
+ bitdepth = 32;
+ break;
+ default:
+ printf("No known bitdepth! Exiting ...\n");
+ return -1;
+ }
+
+ Calculator calculator(
+ sf_info.channels,
+ sf_info.samplerate,
+ bitdepth,
+ channel_corrections,
+ buffer_size_ms,
+ number_of_filter_interpolation_points,
+ num_cpu
+ );
+
+ while (true) {
+ std::vector<double> buffer(4096);
+ auto read = sf_read_double(file, buffer.data(), buffer.size());
+ if (read <= 0) {
+ break;
+ }
+
+ buffer.resize(read);
+ calculator.add(buffer);
+ }
+
+ sf_close(file);
+
+ return {calculator.leq_m(), calculator.leq_nw()};
+}
+
+
+std::vector<double> calculate_ir(double number_of_filter_interpolation_points, int sample_rate, int bits_per_sample)
+{
+ return inverse_fft(convert_log_to_linear(equalinterval2(number_of_filter_interpolation_points, sample_rate, bits_per_sample)));
+}
+
+
+/** @return List of default channel corrections for the given number of channels,
+ * or an empty vector if we can't offer a default.
+ */
+std::vector<double> default_channel_corrections(int channels)
+{
+ std::vector<double> corr;
+
+ if (channels == 6) {
+ double conf51[] = {0, 0, 0, 0, -3, -3};
+ for (auto cind = 0; cind < channels; cind++) {
+ corr.push_back(convert_log_to_linear_single(conf51[cind]));
+ }
+ }
+
+ return corr;
+}
+
+
+double leqm_nrt::convert_log_to_linear_single(double in)
+{
+ return powf(10, in / 20.0f);
+}
+
+
+Calculator::Calculator(
+ int channels,
+ int sample_rate,
+ int bits_per_sample,
+ std::vector<double> channel_corrections,
+ int buffer_size_ms,
+ int number_of_filter_interpolation_points,
+ int num_cpu
+ )
+ : _channels(channels)
+ , _channel_corrections(channel_corrections)
+ , _number_of_filter_interpolation_points(number_of_filter_interpolation_points)
+ , _num_cpu(num_cpu)
+{
+ if ((sample_rate * buffer_size_ms) % 1000) {
+ throw BadBufferSizeError();
+ }
+
+ if (static_cast<int>(channel_corrections.size()) != channels) {
+ channel_corrections = default_channel_corrections(channels);
+ }
+ if (static_cast<int>(channel_corrections.size()) != channels) {
+ throw BadChannelCorrectionsError();
+ }
+
+ _ir = calculate_ir(number_of_filter_interpolation_points, sample_rate, bits_per_sample);
+ _buffer.resize((sample_rate * channels * buffer_size_ms) / 1000);
+}
+
+
+void Calculator::process_buffer()
+{
+ if (_buffer_free_offset == 0) {
+ return;
+ }
+
+ _workers.push_back(
+ std::make_shared<Worker>(
+ _buffer,
+ _buffer_free_offset,
+ _channels,
+ _number_of_filter_interpolation_points,
+ _ir,
+ &_sum,
+ _channel_corrections
+ )
+ );
+
+ if (static_cast<int>(_workers.size()) == _num_cpu) {
+ _workers.clear();
+ }
+
+ _buffer_free_offset = 0;
+}
+
+
+void Calculator::add(std::vector<double> samples)
+{
+ size_t samples_offset = 0;
+
+ while (samples_offset < samples.size()) {
+ /* Copy some data into our buffer */
+ auto to_copy = std::min(samples.size() - samples_offset, _buffer.size() - _buffer_free_offset);
+ memcpy(_buffer.data() + _buffer_free_offset, samples.data() + samples_offset, to_copy * sizeof(double));
+ samples_offset += to_copy;
+ _buffer_free_offset += to_copy;
+
+ /* Process the buffer if it's full */
+ if (_buffer_free_offset == _buffer.size()) {
+ process_buffer();
+ }
+ }
+}
+
+
+double Calculator::leq_m()
+{
+ process_buffer();
+ _workers.clear();
+
+ return _sum.leqm();
+}
+
+
+double Calculator::leq_nw()
+{
+ process_buffer();
+ _workers.clear();
+
+ return _sum.rms();
+}
--- /dev/null
+#include <string>
+#include <vector>
+#include <functional>
+#include <stdexcept>
+#include <mutex>
+#include <cmath>
+#include <memory>
+
+namespace leqm_nrt {
+
+class Sum
+{
+public:
+ void sum_samples(std::vector<double> const& input_samples, std::vector<double> const& c_input_samples, int nsamples)
+ {
+ _mutex.lock();
+ _nsamples += nsamples;
+ for (auto i = 0; i < nsamples; i++) {
+ _sum += input_samples[i];
+ _csum += c_input_samples[i];
+ }
+ _mutex.unlock();
+ }
+
+ int nsamples() const
+ {
+ return _nsamples;
+ }
+
+ /*
+ How the final offset is calculated without reference to a test tone:
+ P0 is the SPL reference 20 uPa
+ Reference SPL is RMS ! So 85 SPL over 20 uPa is 10^4.25 x 0.000020 = 0.355655882 Pa (RMS),
+ but Peak value is 0.355655882 x sqr(2) = 0.502973372 that is 20 x log ( 0.502973372 / 0.000020) = 88.010299957
+ To that one has to add the 20 dB offset of the reference -20dBFS: 88.010299957 + 20.00 = 108.010299957
+
+ But ISO 21727:2004(E) ask for a reference level "measured using an average responding meter". So reference level is not 0.707, but 0.637 = 2/pi
+ */
+
+ double rms() const
+ {
+ return 20 * log10(mean()) + 108.010299957;
+ }
+
+ double leqm() const
+ {
+ return 20 * log10(cmean()) + 108.010299957;
+ }
+
+private:
+ double mean() const
+ {
+ return pow(_sum / _nsamples, 0.500);
+ }
+
+ double cmean() const
+ {
+ return pow(_csum / _nsamples, 0.500);
+ }
+
+ double _csum = 0.0; // convolved sum
+ double _sum = 0.0; // flat sum
+ int _nsamples = 0;
+ std::mutex _mutex;
+};
+
+
+struct Result
+{
+ Result(int status_)
+ : status(status_)
+ {}
+
+ Result(double leq_m_, double leq_nw_)
+ : status(0)
+ , leq_m(leq_m_)
+ , leq_nw(leq_nw_)
+ {}
+
+ /** 0 on success, or
+ *
+ * -100: Either channel_corrections contained a different number of
+ * calibrations than; number of channels in the file or it was empty and the
+ * program cannot infer one from the number of channels. Please specify a
+ * values in channel_corrections.
+ *
+ * -101: Failed to open the sound file.
+ *
+ * -102: buffer_size_ms is not an integer number of samples at the sound file's rate.
+ */
+ int status;
+
+ double leq_m;
+ double leq_nw;
+};
+
+
+Result calculate_file(
+ std::string sound_filename,
+ std::vector<double> channel_corrections,
+ int buffer_size_ms,
+ int number_of_filter_interpolation_points,
+ int num_cpu
+ );
+
+
+double convert_log_to_linear_single(double in);
+
+
+class BadBufferSizeError : public std::runtime_error
+{
+public:
+ BadBufferSizeError()
+ : std::runtime_error("Buffer size does not correspond to an integer number of samples")
+ {}
+};
+
+
+class BadChannelCorrectionsError : public std::runtime_error
+{
+public:
+ BadChannelCorrectionsError()
+ : std::runtime_error("Incorrect number of channel corrections given, and no defaults are available")
+ {}
+};
+
+
+class Worker;
+
+
+class Calculator
+{
+public:
+ Calculator(
+ int channels,
+ int sample_rate,
+ int bits_per_sample,
+ std::vector<double> channel_corrections,
+ int buffer_size_ms,
+ int number_of_filter_interpolation_points,
+ int num_cpu
+ );
+
+ Calculator(Calculator&) = delete;
+ Calculator(Calculator&&) = delete;
+ bool operator=(Calculator&) = delete;
+ bool operator=(Calculator&&) = delete;
+
+ void add(std::vector<double> samples);
+
+ double leq_m();
+ double leq_nw();
+
+private:
+ void process_buffer();
+
+ int _channels;
+ std::vector<double> _channel_corrections;
+ int _number_of_filter_interpolation_points;
+ int _num_cpu;
+ std::vector<std::shared_ptr<Worker>> _workers;
+ Sum _sum;
+ std::vector<double> _ir;
+ std::vector<double> _buffer;
+ size_t _buffer_free_offset = 0;
+};
+
+}
+
--- /dev/null
+/*
+ leqm-nrt is a non-real-time implementation
+ of Leq(M) measurement according to ISO 21727:2004(E)
+ "Cinematography -- Method of measurement of perceived
+ loudness of motion-picture audio material"
+
+ Copyright (C) 2011-2013, 2017-2018 Luca Trisciani
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ */
+
+#include "leqm-nrt.h"
+
+#include <stdio.h>
+#include <math.h>
+#include <sndfile.h>
+#include <unistd.h>
+#include <pthread.h>
+#include <string.h>
+#include <stdlib.h>
+#include <time.h>
+#include <ctype.h>
+#include <iso646.h>
+
+#include <vector>
+#include <memory>
+#include <iostream>
+#include <thread>
+#include <mutex>
+#include <functional>
+
+// Version 0.0.18 (C) Luca Trisciani 2011-2013, 2017-2018
+// Tool from the DCP-Werkstatt Software Bundle
+
+
+
+// COMPILATION
+// compile for DEBUG with gcc -g -DEBUG -lsndfile -lfftw3 -lm -lpthread -lrt -o leqm-nrt leqm-nrt.cpp
+//otherwise gcc -lsndfile -lm -lpthread -lrt -o leqm-nrt leqm-nrt.c
+
+//#define DEBUG
+
+int main(int argc, const char ** argv)
+{
+ int number_of_filter_interpolation_points = 64; // This value is low for precision. Calibration is done with 32768 point.
+ int num_cpu = std::thread::hardware_concurrency() - 1;
+
+ printf("leqm-nrt Copyright (C) 2011-2013, 2017-2018 Luca Trisciani\nThis program comes with ABSOLUTELY NO WARRANTY.\nThis is free software, and you are welcome to redistribute it\nunder the GPL v3 licence.\nProgram will use 1 + %d slave threads.\n", num_cpu);
+ int buffer_size_ms = 850; //ISO 21727:2004 do not contain any indication, TASA seems to indicate 1000, p. 8
+ std::vector<double> channel_corrections;
+ int parameterstate = 0;
+ bool display_leqnw = false;
+ std::string sound_filename;
+
+ if (argc == 1)
+ { const char helptext[] = "Order of parameters is free.\nPossible parameters are:\n-convpoints <integer number> \tNumber of interpolation points for the filter. Default 64\n-numcpus <integer number> \tNumber of slave threads to speed up operation.\n-timing \t\t\tFor benchmarking speed.\n-leqnw\t Print out Leq without M Weighting\n-chconfcal <db correction> <db correction> <etc. so many times as channels>\n-logleqm10\n-logleqm\n-buffersize <milliseconds>\n";
+ printf(helptext);
+ printf("Please indicate a sound file to be processed.\n");
+ return 0;
+ }
+
+ for (int in = 1; in < argc;) {
+
+ if (!(strncmp(argv[in], "-", 1) == 0)) {
+ sound_filename = argv[in];
+ in++;
+ }
+
+ if (strcmp(argv[in], "-chconfcal") == 0) {
+ /* as the order of parameter is free I have to postpone
+ the check for consistency with the number of channels.
+ So first create a temporary array, whose number of element will be checked after
+ the parsing of the command line parameters is finished.
+ The calibration will be expressed in dB on the command line and converted to multiplier
+ here so that it can be stored as a factor in the channelconfcalvector.
+ */
+
+ in++;
+ for (;;) {
+ if (in < argc) {
+ if (!(strncmp(argv[in], "-", 1) == 0) || isdigit(argv[in][1])) {
+ channel_corrections.push_back(leqm_nrt::convert_log_to_linear_single(atof(argv[in++])));
+ } else break;
+ } else break;
+
+ } //for
+ continue;
+ }
+
+ if (strcmp(argv[in], "-convpoints") == 0) {
+ number_of_filter_interpolation_points = atoi(argv[in + 1]);
+ in+=2;
+ printf("Convolution points sets to %d.\n", number_of_filter_interpolation_points);
+ continue;
+
+ }
+
+
+ if (strcmp(argv[in], "-version") == 0) {
+ in++;
+ printf("leqm-nrt version 0.18\n");
+ continue;
+
+ }
+ if (strcmp(argv[in], "-numcpus") == 0) {
+ num_cpu= atoi(argv[in + 1]);
+ in+=2;
+ printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", num_cpu);
+ continue;
+
+ }
+
+ if (strcmp(argv[in], "-leqnw") == 0) {
+ display_leqnw = true;
+ in++;
+ printf("Leq(nW) - unweighted - will be outputted.\n");
+ continue;
+
+ }
+
+ if (strcmp(argv[in], "-buffersize") == 0) {
+ buffer_size_ms = atoi(argv[in + 1]);
+ in+=2;
+ printf("Buffersize will be set to %d milliseconds.\n", buffer_size_ms);
+ continue;
+
+ }
+
+ if (parameterstate==0) {
+ break;
+ }
+ }
+
+ auto result = leqm_nrt::calculate_file(sound_filename, channel_corrections, buffer_size_ms, number_of_filter_interpolation_points, num_cpu);
+
+ if (display_leqnw) {
+ printf("Leq(nW): %.4f\n", result.leq_nw); // Leq(no Weighting)
+ }
+ printf("Leq(M): %.4f\n", result.leq_m);
+
+ return result.status;
+}
+
def build(bld):
obj = bld(features='cxx cxxshlib')
- obj.name = 'libleqm-nrt'
- obj.target = 'leqm-nrt'
+ obj.name = 'libleqm_nrt'
+ obj.target = 'leqm_nrt'
obj.export_includes = ['.']
- obj.source = 'leqm-nrt.cc'
+ obj.source = 'leqm_nrt.cc'
obj = bld(features='cxx cxxprogram')
- obj.name = 'leqm-nrt-cli'
- obj.target = 'leqm-nrt'
- obj.source = 'leqm-nrt-cli.cc'
- obj.use = 'libleqm-nrt'
+ obj.name = 'leqm_nrt_cli'
+ obj.target = 'leqm_nrt'
+ obj.source = 'leqm_nrt_cli.cc'
+ obj.use = 'libleqm_nrt'
obj.uselib = 'SNDFILE'
- bld.install_files('${PREFIX}/include', 'leqm-nrt.h')
+ bld.install_files('${PREFIX}/include', 'leqm_nrt.h')
conf.check_cfg(package='sndfile', args='--cflags --libs', uselib_store='SNDFILE')
def build(bld):
- bld(source='libleqm-nrt.pc.in',
+ bld(source='libleqm_nrt.pc.in',
version=VERSION,
includedir='%s/include/libdcp%s' % (bld.env.PREFIX, bld.env.API_VERSION),
- libs="-L${libdir} -lleqm-nrt",
+ libs="-L${libdir} -lleqm_nrt",
install_path='${LIBDIR}/pkgconfig')
bld.recurse('src')
projects / leqm-nrt.git / commitdiff