/*
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 .
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#ifdef _WIN32
#include
#elif __APPLE__
#include
#include csum = 0.0;
totsum->sum = 0.0;
totsum->nsamples = 0;
totsum->cmean = 0.0;
totsum->mean = 0.0; // Do I write anything here?
totsum->leqm = 0.0;
totsum->rms = 0.0;
sf_count_t samples_read = 0;
// Main loop through audio file
int worker_id = 0;
pthread_t tid[numCPU];
struct WorkerArgs ** WorkerArgsArray;
WorkerArgsArray = malloc(sizeof(struct WorkerArgs *)*numCPU);
int staindex = 0; //shorttermarrayindex
while((samples_read = sf_read_double(file, buffer, buffersizesamples)) > 0) {
WorkerArgsArray[worker_id]=malloc(sizeof(struct WorkerArgs));
WorkerArgsArray[worker_id]->nsamples = samples_read;
WorkerArgsArray[worker_id]->nch = sfinfo.channels;
WorkerArgsArray[worker_id]->npoints=npoints;
WorkerArgsArray[worker_id]->ir = ir;
WorkerArgsArray[worker_id]->ptrtotsum = totsum;
WorkerArgsArray[worker_id]->chconf = channelconfcalvector;
if (leqm10) {
WorkerArgsArray[worker_id]->shorttermindex = staindex++;
WorkerArgsArray[worker_id]->leqm10flag = 1;
WorkerArgsArray[worker_id]->shorttermarray = shorttermaveragedarray;
} else {
WorkerArgsArray[worker_id]->shorttermindex = 0;
WorkerArgsArray[worker_id]->leqm10flag = 0;
}
WorkerArgsArray[worker_id]->argbuffer = malloc(sizeof(double)*buffersizesamples);
memcpy(WorkerArgsArray[worker_id]->argbuffer, buffer, samples_read*sizeof(double));
pthread_attr_t attr;
pthread_attr_init(&attr);
pthread_create(&tid[worker_id], &attr, worker_function, WorkerArgsArray[worker_id]);
worker_id++;
if (worker_id == numCPU) {
worker_id = 0;
//maybe here wait for all cores to output before going on
for (int idxcpu = 0; idxcpu < numCPU; idxcpu++) {
pthread_join(tid[idxcpu], NULL);
free(WorkerArgsArray[idxcpu]->argbuffer);
WorkerArgsArray[idxcpu]->argbuffer = NULL;
free(WorkerArgsArray[idxcpu]);
WorkerArgsArray[idxcpu] = NULL;
}
//simply log here your measurement it will be a multiple of your threads and your buffer
if (leqmlog) {
meanoverduration(totsum); //update leq(m) until now and log it
logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
} //endlog
}
//end while worker_id
/// End looping cores
} // main loop through file
//here I should wait for rest worker (< numcpu)
//but I need to dispose of thread id.
if (worker_id != 0) { // worker_id = 0 means the number of samples was divisible through the number of cpus
for (int idxcpu = 0; idxcpu < worker_id; idxcpu++) { //worker_id is at this point one unit more than threads launched
pthread_join(tid[idxcpu], NULL);
free(WorkerArgsArray[idxcpu]->argbuffer);
WorkerArgsArray[idxcpu]->argbuffer = NULL;
free(WorkerArgsArray[idxcpu]);
WorkerArgsArray[idxcpu] = NULL;
}
//also log here for a last value
if (leqmlog) {
meanoverduration(totsum); //update leq(m) until now and log it
logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
} //endlog
}
// mean of scalar sum over duration
meanoverduration(totsum);
if (leqnw) {
printf("Leq(nW): %.4f\n", totsum->rms); // Leq(no Weighting)
}
printf("Leq(M): %.4f\n", totsum->leqm);
if(timing) {
struct timespec stoptime;
long stoptimenanoseconds;
long executionnanoseconds;
clock_gettime(CLOCK_MONOTONIC, &stoptime);
if (stoptime.tv_nsec < starttime.tv_nsec) {
stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
} else {
stoptimenanoseconds = stoptime.tv_nsec;
}
executionnanoseconds = stoptimenanoseconds - starttime.tv_nsec;
printf("Total execution time is %.6f seconds\n", ((double) stoptime.tv_sec) - ((double) starttime.tv_sec) + ((double) executionnanoseconds / 1000000000.00));
}
if (leqm10) {
//Take the array with the short term accumulators
double interval = 10.0;
//create a rolling average according to rolling interval
int rollint; // in short 10*60 = 600 sec 600/0.850
//how many element of the array to consider for the rollint?
//that is how many buffersizems in the interval - interval could be parameterized(?)
double tempint = 60.0 * interval / (((double) buffersizems) /1000.0);
rollint = (int) tempint;
//dispose of the rest
if (tempint - ((double) rollint) > 0) {
rollint += 1;
}
//two loops
//external loop
int indexlong = 0;
while(indexlong < (numbershortperiods - rollint)) {
double accumulator = 0;
//internal loop
double averagedaccumulator = 0;
for (int indexshort = 0; indexshort < rollint; indexshort++) {
accumulator += shorttermaveragedarray[indexshort+indexlong];
} //end internal loop
averagedaccumulator = accumulator/((double) rollint);
logleqm10(leqm10logfile, ((double) (indexlong+rollint)) * ((double) buffersizems / 1000.0), averagedaccumulator);
indexlong++;
} //end external loop
fclose(leqm10logfile);
free(shorttermaveragedarray);
shorttermaveragedarray = NULL;
}
if (leqmlog) {
fclose(leqmlogfile);
}
sf_close(file);
free(eqfreqsamples);
eqfreqsamples = NULL;
free(eqfreqresp_db);
eqfreqresp_db=NULL;
free(eqfreqresp);
eqfreqresp = NULL;
free(ir);
ir = NULL;
free(channelconfcalvector);
channelconfcalvector = NULL;
free(WorkerArgsArray);
WorkerArgsArray = NULL;
free(totsum);
totsum = NULL;
free(buffer);
buffer=NULL;
return 0;
}
void * worker_function(void * argstruct) {
struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
double * sumandsquarebuffer;
double * csumandsquarebuffer;
double * chsumaccumulator_norm;
double * chsumaccumulator_conv;
sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
for (int i = 0; i < thisWorkerArgs->nsamples / thisWorkerArgs->nch; i++) {
sumandsquarebuffer[i] = 0.0;
csumandsquarebuffer[i] = 0.0;
chsumaccumulator_norm[i] = 0.0;
chsumaccumulator_conv[i] = 0.0;
}
for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
double * normalizedbuffer;
double * convolvedbuffer;
normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
for (int n=ch, m= 0; n < thisWorkerArgs->nsamples; n += thisWorkerArgs->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
normalizedbuffer[m] = thisWorkerArgs->argbuffer[n]*thisWorkerArgs->chconf[ch]; //this scale amplitude according to specified calibration
}
//convolution
convolv_buff(normalizedbuffer, convolvedbuffer, thisWorkerArgs->ir, thisWorkerArgs->nsamples / thisWorkerArgs->nch, thisWorkerArgs->npoints * 2);
//rectify, square und sum
rectify(csumandsquarebuffer,convolvedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
rectify(sumandsquarebuffer,normalizedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
free(normalizedbuffer);
normalizedbuffer= NULL;
free(convolvedbuffer);
convolvedbuffer=NULL;
} // loop through channels
//Create a function for this also a tag so that the worker know if he has to do this or not
if (thisWorkerArgs->leqm10flag) {
thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
#ifdef DEBUG
printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
#endif
}
pthread_mutex_lock(&mutex);
// this should be done under mutex conditions -> shared resources!
sumsamples(thisWorkerArgs->ptrtotsum, chsumaccumulator_norm, chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
pthread_mutex_unlock(&mutex);
free(sumandsquarebuffer);
sumandsquarebuffer=NULL;
free(csumandsquarebuffer);
csumandsquarebuffer=NULL;
free(chsumaccumulator_norm);
chsumaccumulator_norm=NULL;
free(chsumaccumulator_conv);
chsumaccumulator_conv=NULL;
free(thisWorkerArgs->argbuffer);
thisWorkerArgs->argbuffer = NULL;
// the memory pointed to by this pointer is freed in main
// it is the same memory for all worker
// but it is necessary to set pointer to NULL otherwise free will not work later (really?)
thisWorkerArgs->chconf = NULL;
pthread_exit(0);
}
//to get impulse response frequency response at equally spaced intervals is needed
int equalinterval( double * freqsamples, double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
double freq;
// int findex = 0;
// int rindex = 0;
double pass = ((double) (samplingfreq >> 1)) / ((double) points);
for (int ieq = 0, i = 0; ieq < points; ieq++) {
freq = ieq*pass;
eqfreqsamples[ieq] = freq;
if ((freq == 0.0) || (freq < freqsamples[1])) {
eqfreqresp[ieq] = freqresp[0];
continue;
} else {
if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp[i];
} else if (freq >=freqsamples[i+1]) {
while(freq >= freqsamples[i+1]) {
i++;
if ((i + 1) >= origpoints) {
break;
}
}
if ((i+1) < origpoints) {
eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
} else {
eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
}
}
}
}
return 0;
}
//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
int equalinterval2( double freqsamples[], double freqresp_db[], double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile) {
double freq;
//calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
double dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
//double dcatt = -90.3;
double pass = ((double) (samplingfreq >> 1)) / ((double) points);
for (int ieq = 0, i = 0; ieq < points; ieq++) {
freq = ieq*pass;
eqfreqsamples[ieq] = freq;
if (freq == 0.0) {
eqfreqresp[ieq] = dcatt;
} else if (freq < freqsamples[0]) { // this has a lot of influence on final Leq(M) value
eqfreqresp[ieq] = ((freqresp_db[0] - dcatt) / (freqsamples[0] - 0)) * freq + dcatt;
//eqfreqresp[ieq] = freqresp_db[0]; // Is this meaningful? Shouldn't I interpolate between 0 Hz and 31 Hz? Otherwise for DC I have -35.5 dB
continue;
} else {
if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp_db[i];
} else if (freq >=freqsamples[i+1]) {
while(freq >= freqsamples[i+1]) {
i++;
if ((i + 1) >= origpoints) {
break;
}
}
if ((i+1) < origpoints) {
eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
} else {
eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
}
}
}
}
return 0;
}
int convloglin(double * in, double * out, int points) {
for (int i = 0; i < points; i++) {
out[i] = powf(10, (in[i]/20.0));
}
return 0;
}
double convlinlog_single(double in) {
double out;
out = log(in)*20.0f;
return out;
}
double convloglin_single(double in) {
double out;
out = powf(10, in/20.0f);
return out;
}
// convolution
int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
double sum = 0.0;
for (int i = 0; i < sigin_dim; i++) {
int m = i;
for (int l = impresp_dim - 1; l >=0; l--,m++) {
if (m >= sigin_dim) {
m -= sigin_dim;
}
sum += sigin[m]*impresp[l];
}
sigout[i] = sum;
sum=0.0;
}
return 0;
}
void inversefft2(double * eqfreqresp, double * ir, int npoints) {
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*((double) m)*( ( ((double) n) - ( ((double) npoints) * 2.0 -1 ) / 2 ) / ( ((double) npoints) * 2.0) ));
parsum = parsum + eqfreqresp[m]*partial;
}
ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
#ifdef DEBUG
printf("%.4f\n", ir[n]);
#endif
}
for (int n = 0; n < npoints; n++) {
ir[npoints+n] = ir[npoints-(n + 1)];
#ifdef DEBUG
printf("%.4f\n", ir[npoints+n]);
#endif
}
}
// scale input according to required calibration
// this could be different for certain digital cinema formats
double inputcalib (double dbdiffch) {
double coeff = pow(10, dbdiffch/20);
return coeff;
}
//rectify, square and sum
int rectify(double * squared, double * inputsamples, int nsamples){
for (int i = 0; i < nsamples; i++) {
squared[i] = (double) powf(inputsamples[i], 2);
}
return 0;
}
int initbuffer(double * buffertoinit, int nsamples) {
for (int i = 0; i < nsamples; i++) {
buffertoinit[i] = 0.0;
}
return 0;
}
int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
for (int i = 0; i < nsamples; i++) {
chaccumulator[i] += inputchannel[i];
}
return 0;
}
int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples) {
ts->nsamples += nsamples;
for (int i=0; i < nsamples; i++) {
ts->sum += inputsamples[i];
ts->csum += cinputsamples[i];
}
return 0;
}
int meanoverduration(struct Sum * oldsum) {
oldsum->mean = pow(oldsum->sum / ((double) oldsum->nsamples), 0.500);
oldsum->cmean = pow(oldsum->csum / ((double) oldsum->nsamples), 0.500);
oldsum->rms = 20*log10(oldsum->mean) + 108.010299957;
oldsum->leqm = 20*log10(oldsum->cmean) + 108.010299957;//
/*
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
*/
return 0;
}
double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
double stsum = 0.0;
for (int i=0; i < nsamples; i++) {
stsum += channelaccumulator[i];
}
return stsum / (double) nsamples;
}
void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
fprintf(filehandle, "%.4f", featuretimesec);
fprintf(filehandle, "\t");
fprintf(filehandle, "%.4f\n", oldsum->leqm);
}
void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
double leqm10 = 20*log10(pow(longaverage, 0.500)) + 108.010299957;
fprintf(filehandle, "%.4f", featuretimesec);
fprintf(filehandle, "\t");
fprintf(filehandle, "%.4f\n", leqm10);
}