2 leqm-nrt is a non-real-time implementation
3 of Leq(M) measurement according to ISO 21727:2004(E)
4 "Cinematography -- Method of measurement of perceived
5 loudness of motion-picture audio material"
7 Copyright (C) 2011-2013, 2017 Luca Trisciani
9 This program is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
19 You should have received a copy of the GNU General Public License
20 along with this program. If not, see <http://www.gnu.org/licenses/>.
40 #include <sys/param.h>
41 #include <sys/sysctl.h
45 // Version 0.0.17 (C) Luca Trisciani 2011-2013, 2017
46 // Tool from the DCP-Werkstatt Software Bundle
51 // compile for DEBUG with gcc -g -DEBUG -lsndfile -lfftw3 -lm -lpthread -lrt -o leqm-nrt leqm-nrt.cpp
52 //otherwise gcc -lsndfile -lm -lpthread -lrt -o leqm-nrt leqm-nrt.c
58 double csum; // convolved sum
59 double sum; // flat sum
61 double cmean; //convolved mean
73 struct Sum * ptrtotsum;
76 double * shorttermarray;
80 int equalinterval( double * freqsamples, double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints);
81 int equalinterval2( double freqsamples[], double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile);
82 int convloglin(double * in, double * out, int points);
83 double convlinlog_single(double in);
84 double convloglin_single(double in);
85 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim);
86 double inputcalib (double dbdiffch);
87 int rectify(double * squared, double * inputsamples, int nsamples);
88 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples);
89 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples);
90 int meanoverduration(struct Sum * oldsum);
91 void inversefft1(double * eqfreqresp, double * ir, int npoints);
92 void inversefft2(double * eqfreqresp, double * ir, int npoints);
93 void * worker_function(void * argfunc);
94 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum);
95 double sumandshorttermavrg(double * channelaccumulator, int nsamples);
96 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage);
99 pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
102 int main(int argc, const char ** argv)
104 int npoints = 64; // This value is low for precision. Calibration is done with 32768 point.
105 int origpoints = 21; //number of points in the standard CCIR filter
106 int samplingfreq; // this and the next is defined later taking it from sound file
108 // double normalizer;
110 struct timespec starttime;
111 int fileopenstate = 0;
114 #if defined __unix__ || defined __APPLE__
115 int numCPU = sysconf(_SC_NPROCESSORS_ONLN) - 1;
116 #elif defined _WIN64 || defined _WIN32
118 GetSystemInfo(&sysinfo);
119 int numCPU = sysinfo.dwNumberOfProcessors - 1;
122 double * channelconfcalvector;
123 channelconfcalvector = NULL;
124 printf("leqm-nrt Copyright (C) 2011-2013, 2017 Luca Trisciani\nThis program comes with ABSOLUTELY NO WARRANTY; for details on command line parameters -help\nThis is free software, and you are welcome to redistribute it\nunder the GPL v3 licence.\nProgram will use 1 + %d slave threads.\n", numCPU);
125 //SndfileHandle file;
130 leqm10logfile = NULL;
133 int buffersizems = 850; //ISO 21727:2004 do not contain any indication, TASA seems to indicate 1000, p. 8
134 int buffersizesamples;
135 double tempchcal[128];
137 double * shorttermaveragedarray;
138 shorttermaveragedarray = NULL;
139 int numbershortperiods = 0;
140 int parameterstate = 0;
143 char soundfilename[64];
144 // This is a requirement of sndfile library, do not forget it.
146 memset(&sfinfo, 0, sizeof(sfinfo));
150 { 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-chconfcal <db correction> <db correction> <etc. so many times as channels>\n-logleqm10\n-logleqm\n-buffersize <milliseconds>\n";
152 printf("Please indicate a sound file to be processed.\n");
157 for (int in = 1; in < argc;) {
159 if (!(strncmp(argv[in], "-", 1) == 0)) {
160 if (fileopenstate == 0) {
161 if(! (file = sf_open(argv[in], SFM_READ, &sfinfo))) {
162 printf("Error while opening audio file, could not open %s\n.", argv[in]);
163 puts(sf_strerror(NULL));
167 strcpy(soundfilename, argv[in]);
169 printf("Opened file: %s\n", argv[in]);
170 printf("Sample rate: %d\n", sfinfo.samplerate);
171 printf("Channels: %d\n", sfinfo.channels);
172 printf("Format: %d\n", sfinfo.format);
173 printf("Frames: %d\n", (int) sfinfo.frames);
174 channelconfcalvector = malloc(sizeof(double) * sfinfo.channels);
178 free(channelconfcalvector);
179 channelconfcalvector = NULL;
185 if (strcmp(argv[in], "-chconfcal") == 0) {
186 /* as the order of parameter is free I have to postpone
187 the check for consistency with the number of channels.
188 So first create a temporary array, whose number of element will be checked after
189 the parsing of the command line parameters is finished.
190 The calibration will be expressed in dB on the command line and converted to multiplier
191 here so that it can be stored as a factor in the channelconfcalvector.
197 //if (!(strncmp(argv[in], "-", 1) == 0)) { //changed this to allow negative numbers
198 if (!(strncmp(argv[in], "-", 1) == 0) || isdigit(argv[in][1])) {
199 tempchcal[numcalread++]=atof(argv[in++]);
207 if (strcmp(argv[in], "-convpoints") == 0) {
208 npoints = atoi(argv[in + 1]);
210 printf("Convolution points sets to %d.\n", npoints);
214 if (strcmp(argv[in], "-numcpus") == 0) {
215 numCPU= atoi(argv[in + 1]);
217 printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", numCPU);
221 if (strcmp(argv[in], "-timing") == 0) {
224 printf("Execution time will be measured.\n");
229 if (strcmp(argv[in], "-logleqm10") == 0) {
232 printf("Leq(M)10 data will be logged to the file leqm10.txt\n");
236 if (strcmp(argv[in], "-logleqm") == 0) {
239 printf("Leq(M) data will be logged to the file leqmlog.txt\n");
244 if (strcmp(argv[in], "-leqnw") == 0) {
247 printf("Leq(nW) - unweighted - will be outputted.\n");
252 if (strcmp(argv[in], "-buffersize") == 0) {
253 buffersizems = atoi(argv[in + 1]);
255 printf("Buffersize will be set to %d milliseconds.\n", buffersizems);
260 if (parameterstate==0) {
266 //postprocessing parameters
267 if (numcalread == sfinfo.channels) {
268 for (int cind = 0; cind < sfinfo.channels; cind++) {
269 channelconfcalvector[cind] = convloglin_single(tempchcal[cind]);
273 else if ((numcalread == 0) && (sfinfo.channels == 6)) {
274 double conf51[6] = {0, 0, 0, 0, -3, -3};
275 for (int cind = 0; cind < sfinfo.channels; cind++) {
276 channelconfcalvector[cind] = convloglin_single(conf51[cind]);
281 printf("Either you specified a different number of calibration than number of channels in the file or you do not indicate any calibration and the program cannot infer one from the number of channels. Please specify a channel calibration on the command line.\n");
283 free(channelconfcalvector);
284 channelconfcalvector = NULL;
292 char tempstring[128];
293 strcpy(tempstring, soundfilename);
294 strcat(tempstring, ".leqm10.txt");
295 leqm10logfile = fopen(tempstring, "w");
296 if (leqm10logfile == NULL) {
297 printf("Could not open file to write log leqm10 data!\n");
305 char tempstring[128];
306 strcpy(tempstring, soundfilename);
307 strcat(tempstring, ".leqmlog.txt");
308 leqmlogfile = fopen(tempstring, "w");
309 if (leqmlogfile == NULL) {
310 printf("Could not open file to write log leqm data!\n");
317 clock_gettime(CLOCK_MONOTONIC, &starttime);
324 // reading to a double or float buffer with sndfile take care of normalization
326 static double buffer[BUFFER_LEN]; // it seems this must be static. I don't know why
329 // buffer = new double [BUFFER_LEN];
330 buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
331 if ((sfinfo.samplerate*buffersizems)%1000) {
332 printf("Please fine tune the buffersize according to the sample rate\n");
335 // write a function to do that
339 buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
340 buffer = malloc(sizeof(double)*buffersizesamples);
342 samplingfreq = sfinfo.samplerate;
346 //if duration < 10 mm exit
348 double featdursec = sfinfo.frames / sfinfo.samplerate;
349 if ((featdursec/60.0) < 10.0) {
350 printf("The audio file is too short to measure Leq(m10).\n");
354 //how many short periods in overall duration
355 int remainder = sfinfo.frames % (sfinfo.samplerate*buffersizems/1000);
356 if (remainder == 0) numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000);
357 else numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000) + 1;
360 shorttermaveragedarray = malloc(sizeof(*shorttermaveragedarray)*numbershortperiods);
364 //End opening audio file
368 double freqsamples[] = {31, 63, 100, 200, 400, 800, 1000, 2000, 3150, 4000, 5000, 6300, 7100, 8000, 9000, 10000, 12500, 14000, 16000, 20000, 31500};
369 double freqresp_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};
371 double * eqfreqresp_db;
372 eqfreqresp_db = malloc(sizeof(*eqfreqresp_db)*npoints);
374 double * eqfreqsamples;
375 eqfreqsamples = malloc(sizeof(*eqfreqsamples)*npoints);
377 eqfreqresp = malloc(sizeof(*eqfreqresp)*npoints);
379 ir = malloc(sizeof(*ir)*npoints*2);
382 // And what to do for floating point sample coding?
384 switch(sfinfo.format & SF_FORMAT_SUBMASK) {
385 // all signed bitdepth
399 printf("No known bitdepth! Exiting ...\n");
406 equalinterval2(freqsamples, freqresp_db, eqfreqsamples, eqfreqresp_db, npoints, samplingfreq, origpoints, bitdepth);
407 convloglin(eqfreqresp_db, eqfreqresp, npoints);
410 for (int i=0; i < npoints; i++) {
411 printf("%d\t%.2f\t%.2f\t%.2f\n", i, eqfreqsamples[i], eqfreqresp_db[i], eqfreqresp[i]);
415 inversefft2(eqfreqresp, ir, npoints);
417 // read through the entire file
420 totsum = malloc(sizeof(struct Sum));
423 totsum->nsamples = 0;
425 totsum->mean = 0.0; // Do I write anything here?
428 sf_count_t samples_read = 0;
430 // Main loop through audio file
433 pthread_t tid[numCPU];
434 struct WorkerArgs ** WorkerArgsArray;
435 WorkerArgsArray = malloc(sizeof(struct WorkerArgs *)*numCPU);
436 int staindex = 0; //shorttermarrayindex
439 while((samples_read = sf_read_double(file, buffer, buffersizesamples)) > 0) {
441 WorkerArgsArray[worker_id]=malloc(sizeof(struct WorkerArgs));
442 WorkerArgsArray[worker_id]->nsamples = samples_read;
443 WorkerArgsArray[worker_id]->nch = sfinfo.channels;
444 WorkerArgsArray[worker_id]->npoints=npoints;
445 WorkerArgsArray[worker_id]->ir = ir;
446 WorkerArgsArray[worker_id]->ptrtotsum = totsum;
448 WorkerArgsArray[worker_id]->chconf = channelconfcalvector;
450 WorkerArgsArray[worker_id]->shorttermindex = staindex++;
451 WorkerArgsArray[worker_id]->leqm10flag = 1;
452 WorkerArgsArray[worker_id]->shorttermarray = shorttermaveragedarray;
454 WorkerArgsArray[worker_id]->shorttermindex = 0;
455 WorkerArgsArray[worker_id]->leqm10flag = 0;
458 WorkerArgsArray[worker_id]->argbuffer = malloc(sizeof(double)*buffersizesamples);
459 memcpy(WorkerArgsArray[worker_id]->argbuffer, buffer, samples_read*sizeof(double));
463 pthread_attr_init(&attr);
464 pthread_create(&tid[worker_id], &attr, worker_function, WorkerArgsArray[worker_id]);
467 if (worker_id == numCPU) {
469 //maybe here wait for all cores to output before going on
470 for (int idxcpu = 0; idxcpu < numCPU; idxcpu++) {
471 pthread_join(tid[idxcpu], NULL);
472 free(WorkerArgsArray[idxcpu]);
473 WorkerArgsArray[idxcpu] = NULL;
475 //simply log here your measurement it will be a multiple of your threads and your buffer
477 meanoverduration(totsum); //update leq(m) until now and log it
478 logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
484 //end while worker_id
485 /// End looping cores
486 } // main loop through file
488 //here I should wait for rest worker (< numcpu)
489 //but I need to dispose of thread id.
490 if (worker_id != 0) { // worker_id = 0 means the number of samples was divisible through the number of cpus
491 for (int idxcpu = 0; idxcpu < worker_id; idxcpu++) { //worker_id is at this point one unit more than threads launched
492 pthread_join(tid[idxcpu], NULL);
493 free(WorkerArgsArray[idxcpu]);
494 WorkerArgsArray[idxcpu] = NULL;
496 //also log here for a last value
498 meanoverduration(totsum); //update leq(m) until now and log it
499 logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
502 // mean of scalar sum over duration
504 meanoverduration(totsum);
506 printf("Leq(noW): %.4f\n", totsum->rms); // Leq(no Weighting)
508 printf("Leq(M): %.4f\n", totsum->leqm);
511 struct timespec stoptime;
512 long stoptimenanoseconds;
513 long executionnanoseconds;
514 clock_gettime(CLOCK_MONOTONIC, &stoptime);
516 if (stoptime.tv_nsec < starttime.tv_nsec) {
517 stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
519 stoptimenanoseconds = stoptime.tv_nsec;
521 executionnanoseconds = stoptimenanoseconds - starttime.tv_nsec;
522 printf("Total execution time is %.6f seconds\n", ((double) stoptime.tv_sec) - ((double) starttime.tv_sec) + ((double) executionnanoseconds / 1000000000.00));
528 //Take the array with the short term accumulators
529 double interval = 10.0;
530 //create a rolling average according to rolling interval
531 int rollint; // in short 10*60 = 600 sec 600/0.850
534 //how many element of the array to consider for the rollint?
535 //that is how many buffersizems in the interval - interval could be parameterized(?)
536 double tempint = 60.0 * interval / (((double) buffersizems) /1000.0);
537 rollint = (int) tempint;
538 //dispose of the rest
539 if (tempint - ((double) rollint) > 0) {
545 while(indexlong < (numbershortperiods - rollint)) {
547 double accumulator = 0;
549 double averagedaccumulator = 0;
550 for (int indexshort = 0; indexshort < rollint; indexshort++) {
552 accumulator += shorttermaveragedarray[indexshort+indexlong];
553 } //end internal loop
554 averagedaccumulator = accumulator/((double) rollint);
555 logleqm10(leqm10logfile, ((double) (indexlong+rollint)) * ((double) buffersizems / 1000.0), averagedaccumulator);
557 } //end external loop
559 fclose(leqm10logfile);
560 free(shorttermaveragedarray);
561 shorttermaveragedarray = NULL;
573 eqfreqsamples = NULL;
580 free(channelconfcalvector);
581 channelconfcalvector = NULL;
582 free(WorkerArgsArray);
583 WorkerArgsArray = NULL;
594 void * worker_function(void * argstruct) {
596 struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
598 double * sumandsquarebuffer;
599 double * csumandsquarebuffer;
600 double * chsumaccumulator_norm;
601 double * chsumaccumulator_conv;
604 sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
607 csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
609 chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
611 chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
615 for (int i = 0; i < thisWorkerArgs->nsamples / thisWorkerArgs->nch; i++) {
616 sumandsquarebuffer[i] = 0.0;
617 csumandsquarebuffer[i] = 0.0;
618 chsumaccumulator_norm[i] = 0.0;
619 chsumaccumulator_conv[i] = 0.0;
624 for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
626 double * normalizedbuffer;
627 double * convolvedbuffer;
630 normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
632 convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
635 for (int n=ch, m= 0; n < thisWorkerArgs->nsamples; n += thisWorkerArgs->nch, m++) {
636 // 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
637 //so not normalized but calibrated
638 normalizedbuffer[m] = thisWorkerArgs->argbuffer[n]*thisWorkerArgs->chconf[ch]; //this scale amplitude according to specified calibration
644 convolv_buff(normalizedbuffer, convolvedbuffer, thisWorkerArgs->ir, thisWorkerArgs->nsamples / thisWorkerArgs->nch, thisWorkerArgs->npoints * 2);
645 //rectify, square und sum
646 rectify(csumandsquarebuffer,convolvedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
647 rectify(sumandsquarebuffer,normalizedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
649 accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
650 accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
653 free(normalizedbuffer);
654 normalizedbuffer= NULL;
656 free(convolvedbuffer);
657 convolvedbuffer=NULL;
659 } // loop through channels
661 //Create a function for this also a tag so that the worker know if he has to do this or not
663 if (thisWorkerArgs->leqm10flag) {
664 thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
666 printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
669 pthread_mutex_lock(&mutex);
670 // this should be done under mutex conditions -> shared resources!
671 sumsamples(thisWorkerArgs->ptrtotsum, chsumaccumulator_norm, chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
672 pthread_mutex_unlock(&mutex);
675 free(sumandsquarebuffer);
676 sumandsquarebuffer=NULL;
678 free(csumandsquarebuffer);
679 csumandsquarebuffer=NULL;
681 free(chsumaccumulator_norm);
682 chsumaccumulator_norm=NULL;
684 free(chsumaccumulator_conv);
685 chsumaccumulator_conv=NULL;
687 free(thisWorkerArgs->argbuffer);
688 thisWorkerArgs->argbuffer = NULL;
689 // the memory pointed to by this pointer is freed in main
690 // it is the same memory for all worker
691 // but it is necessary to set pointer to NULL otherwise free will not work later (really?)
692 thisWorkerArgs->chconf = NULL;
698 //to get impulse response frequency response at equally spaced intervals is needed
700 int equalinterval( double * freqsamples, double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
704 double pass = ((double) (samplingfreq >> 1)) / ((double) points);
705 for (int ieq = 0, i = 0; ieq < points; ieq++) {
707 eqfreqsamples[ieq] = freq;
709 if ((freq == 0.0) || (freq < freqsamples[1])) {
710 eqfreqresp[ieq] = freqresp[0];
714 if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
715 eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp[i];
716 } else if (freq >=freqsamples[i+1]) {
717 while(freq >= freqsamples[i+1]) {
719 if ((i + 1) >= origpoints) {
723 if ((i+1) < origpoints) {
724 eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
726 eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
738 //the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
739 //it is also different for the way it interpolates between DC and 31 Hz
740 // Pay attention that also arguments to the functions are changed
741 int equalinterval2( double freqsamples[], double freqresp_db[], double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile) {
745 //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
746 double dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
747 //double dcatt = -90.3;
748 double pass = ((double) (samplingfreq >> 1)) / ((double) points);
749 for (int ieq = 0, i = 0; ieq < points; ieq++) {
751 eqfreqsamples[ieq] = freq;
754 eqfreqresp[ieq] = dcatt;
755 } else if (freq < freqsamples[0]) { // this has a lot of influence on final Leq(M) value
756 eqfreqresp[ieq] = ((freqresp_db[0] - dcatt) / (freqsamples[0] - 0)) * freq + dcatt;
757 //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
761 if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
762 eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp_db[i];
763 } else if (freq >=freqsamples[i+1]) {
764 while(freq >= freqsamples[i+1]) {
766 if ((i + 1) >= origpoints) {
770 if ((i+1) < origpoints) {
771 eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
773 eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
786 int convloglin(double * in, double * out, int points) {
787 for (int i = 0; i < points; i++) {
788 out[i] = powf(10, (in[i]/20.0));
795 double convlinlog_single(double in) {
802 double convloglin_single(double in) {
804 out = powf(10, in/20.0f);
810 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
814 for (int i = 0; i < sigin_dim; i++) {
817 for (int l = impresp_dim - 1; l >=0; l--,m++) {
818 if (m >= sigin_dim) {
821 sum += sigin[m]*impresp[l];
831 void inversefft2(double * eqfreqresp, double * ir, int npoints) {
832 for (int n = 0; n < npoints; n++) {
834 double partial = 0.0;
836 for (int m = 1; m <= npoints -1; m++) {
837 partial = cos(2.0*M_PI*((double) m)*( ( ((double) n) - ( ((double) npoints) * 2.0 -1 ) / 2 ) / ( ((double) npoints) * 2.0) ));
838 parsum = parsum + eqfreqresp[m]*partial;
840 ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
842 printf("%.4f\n", ir[n]);
845 for (int n = 0; n < npoints; n++) {
846 ir[npoints+n] = ir[npoints-(n + 1)];
848 printf("%.4f\n", ir[npoints+n]);
855 // scale input according to required calibration
856 // this could be different for certain digital cinema formats
857 double inputcalib (double dbdiffch) {
859 double coeff = pow(10, dbdiffch/20);
864 //rectify, square and sum
865 int rectify(double * squared, double * inputsamples, int nsamples){
866 for (int i = 0; i < nsamples; i++) {
867 squared[i] = (double) powf(inputsamples[i], 2);
873 int initbuffer(double * buffertoinit, int nsamples) {
874 for (int i = 0; i < nsamples; i++) {
875 buffertoinit[i] = 0.0;
881 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
882 for (int i = 0; i < nsamples; i++) {
883 chaccumulator[i] += inputchannel[i];
888 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples) {
889 ts->nsamples += nsamples;
890 for (int i=0; i < nsamples; i++) {
891 ts->sum += inputsamples[i];
892 ts->csum += cinputsamples[i];
898 int meanoverduration(struct Sum * oldsum) {
899 oldsum->mean = pow(oldsum->sum / ((double) oldsum->nsamples), 0.500);
900 oldsum->cmean = pow(oldsum->csum / ((double) oldsum->nsamples), 0.500);
901 oldsum->rms = 20*log10(oldsum->mean);
902 oldsum->leqm = 20*log10(oldsum->cmean) + 110.600;//
903 // and this must be right because M filter is -5.6 @ 1k Hz that is -25.6 dBFS and to have 85.0 as reference level we must add 25.56 + 85.00 that is 110.6 dB.
904 //this value is obtained calibrating with a -20 dBFS Dolby Tone (RMS) I think this is correct
905 //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
906 //This is only approximate as you should use a separate calibration according to the Dolby Format. Also for SW the tone should be 100Hz (?)
911 double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
913 for (int i=0; i < nsamples; i++) {
914 stsum += channelaccumulator[i];
917 return stsum / (double) nsamples;
920 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
922 fprintf(filehandle, "%.4f", featuretimesec);
923 fprintf(filehandle, "\t");
924 fprintf(filehandle, "%.4f\n", oldsum->leqm);
929 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
930 double leqm10 = 20*log10(pow(longaverage, 0.500)) + 110.600;
931 fprintf(filehandle, "%.4f", featuretimesec);
932 fprintf(filehandle, "\t");
933 fprintf(filehandle, "%.4f\n", leqm10);