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/>.
37 // Version 0.0.17 (C) Luca Trisciani 2011-2013, 2017
38 // Tool from the DCP-Werkstatt Software Bundle
43 // compile for DEBUG with gcc -g -DEBUG -lsndfile -lfftw3 -lm -lpthread -lrt -o leqm-nrt leqm-nrt.cpp
44 //otherwise gcc -lsndfile -lm -lpthread -lrt -o leqm-nrt leqm-nrt.c
50 double csum; // convolved sum
51 double sum; // flat sum
53 double cmean; //convolved mean
65 struct Sum * ptrtotsum;
68 double * shorttermarray;
72 int equalinterval( double * freqsamples, double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints);
73 int equalinterval2( double freqsamples[], double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile);
74 int convloglin(double * in, double * out, int points);
75 double convlinlog_single(double in);
76 double convloglin_single(double in);
77 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim);
78 double inputcalib (double dbdiffch);
79 int rectify(double * squared, double * inputsamples, int nsamples);
80 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples);
81 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples);
82 int meanoverduration(struct Sum * oldsum);
83 void inversefft1(double * eqfreqresp, double * ir, int npoints);
84 void inversefft2(double * eqfreqresp, double * ir, int npoints);
85 void * worker_function(void * argfunc);
86 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum);
87 double sumandshorttermavrg(double * channelaccumulator, int nsamples);
88 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage);
91 pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
94 int main(int argc, const char ** argv)
96 int npoints = 64; // This value is low for precision. Calibration is done with 32768 point.
97 int origpoints = 21; //number of points in the standard CCIR filter
98 int samplingfreq; // this and the next is defined later taking it from sound file
100 // double normalizer;
102 struct timespec starttime;
103 int fileopenstate = 0;
106 #if defined __unix__ || defined __APPLE__
107 int numCPU = sysconf(_SC_NPROCESSORS_ONLN) - 1;
108 #elif defined _WIN64 || defined _WIN32
110 GetSystemInfo(&sysinfo);
111 int numCPU = sysinfo.dwNumberOfProcessors - 1;
114 double * channelconfcalvector;
115 channelconfcalvector = NULL;
116 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);
117 //SndfileHandle file;
122 leqm10logfile = NULL;
125 int buffersizems = 850; //ISO 21727:2004 do not contain any indication, TASA seems to indicate 1000, p. 8
126 int buffersizesamples;
127 double tempchcal[128];
129 double * shorttermaveragedarray;
130 shorttermaveragedarray = NULL;
131 int numbershortperiods = 0;
132 int parameterstate = 0;
135 char soundfilename[64];
136 // This is a requirement of sndfile library, do not forget it.
138 memset(&sfinfo, 0, sizeof(sfinfo));
142 { 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";
144 printf("Please indicate a sound file to be processed.\n");
149 for (int in = 1; in < argc;) {
151 if (!(strncmp(argv[in], "-", 1) == 0)) {
152 if (fileopenstate == 0) {
153 if(! (file = sf_open(argv[in], SFM_READ, &sfinfo))) {
154 printf("Error while opening audio file, could not open %s\n.", argv[in]);
155 puts(sf_strerror(NULL));
159 strcpy(soundfilename, argv[in]);
161 printf("Opened file: %s\n", argv[in]);
162 printf("Sample rate: %d\n", sfinfo.samplerate);
163 printf("Channels: %d\n", sfinfo.channels);
164 printf("Format: %d\n", sfinfo.format);
165 printf("Frames: %d\n", (int) sfinfo.frames);
166 channelconfcalvector = malloc(sizeof(double) * sfinfo.channels);
170 free(channelconfcalvector);
171 channelconfcalvector = NULL;
177 if (strcmp(argv[in], "-chconfcal") == 0) {
178 /* as the order of parameter is free I have to postpone
179 the check for consistency with the number of channels.
180 So first create a temporary array, whose number of element will be checked after
181 the parsing of the command line parameters is finished.
182 The calibration will be expressed in dB on the command line and converted to multiplier
183 here so that it can be stored as a factor in the channelconfcalvector.
189 //if (!(strncmp(argv[in], "-", 1) == 0)) { //changed this to allow negative numbers
190 if (!(strncmp(argv[in], "-", 1) == 0) || isdigit(argv[in][1])) {
191 tempchcal[numcalread++]=atof(argv[in++]);
199 if (strcmp(argv[in], "-convpoints") == 0) {
200 npoints = atoi(argv[in + 1]);
202 printf("Convolution points sets to %d.\n", npoints);
206 if (strcmp(argv[in], "-numcpus") == 0) {
207 numCPU= atoi(argv[in + 1]);
209 printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", numCPU);
213 if (strcmp(argv[in], "-timing") == 0) {
216 printf("Execution time will be measured.\n");
221 if (strcmp(argv[in], "-logleqm10") == 0) {
224 printf("Leq(M)10 data will be logged to the file leqm10.txt\n");
228 if (strcmp(argv[in], "-logleqm") == 0) {
231 printf("Leq(M) data will be logged to the file leqmlog.txt\n");
236 if (strcmp(argv[in], "-leqnw") == 0) {
239 printf("Leq(nW) - unweighted - will be outputted.\n");
244 if (strcmp(argv[in], "-buffersize") == 0) {
245 buffersizems = atoi(argv[in + 1]);
247 printf("Buffersize will be set to %d milliseconds.\n", buffersizems);
252 if (parameterstate==0) {
258 //postprocessing parameters
259 if (numcalread == sfinfo.channels) {
260 for (int cind = 0; cind < sfinfo.channels; cind++) {
261 channelconfcalvector[cind] = convloglin_single(tempchcal[cind]);
265 else if ((numcalread == 0) && (sfinfo.channels == 6)) {
266 double conf51[6] = {0, 0, 0, 0, -3, -3};
267 for (int cind = 0; cind < sfinfo.channels; cind++) {
268 channelconfcalvector[cind] = convloglin_single(conf51[cind]);
273 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");
275 free(channelconfcalvector);
276 channelconfcalvector = NULL;
284 char tempstring[128];
285 strcpy(tempstring, soundfilename);
286 strcat(tempstring, ".leqm10.txt");
287 leqm10logfile = fopen(tempstring, "w");
288 if (leqm10logfile == NULL) {
289 printf("Could not open file to write log leqm10 data!\n");
297 char tempstring[128];
298 strcpy(tempstring, soundfilename);
299 strcat(tempstring, ".leqmlog.txt");
300 leqmlogfile = fopen(tempstring, "w");
301 if (leqmlogfile == NULL) {
302 printf("Could not open file to write log leqm data!\n");
309 clock_gettime(CLOCK_MONOTONIC, &starttime);
316 // reading to a double or float buffer with sndfile take care of normalization
318 static double buffer[BUFFER_LEN]; // it seems this must be static. I don't know why
321 // buffer = new double [BUFFER_LEN];
322 buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
323 if ((sfinfo.samplerate*buffersizems)%1000) {
324 printf("Please fine tune the buffersize according to the sample rate\n");
327 // write a function to do that
331 buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
332 buffer = malloc(sizeof(double)*buffersizesamples);
334 samplingfreq = sfinfo.samplerate;
338 //if duration < 10 mm exit
340 double featdursec = sfinfo.frames / sfinfo.samplerate;
341 if ((featdursec/60.0) < 10.0) {
342 printf("The audio file is too short to measure Leq(m10).\n");
346 //how many short periods in overall duration
347 int remainder = sfinfo.frames % (sfinfo.samplerate*buffersizems/1000);
348 if (remainder == 0) numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000);
349 else numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000) + 1;
352 shorttermaveragedarray = malloc(sizeof(*shorttermaveragedarray)*numbershortperiods);
356 //End opening audio file
360 double freqsamples[] = {31, 63, 100, 200, 400, 800, 1000, 2000, 3150, 4000, 5000, 6300, 7100, 8000, 9000, 10000, 12500, 14000, 16000, 20000, 31500};
361 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};
363 double * eqfreqresp_db;
364 eqfreqresp_db = malloc(sizeof(*eqfreqresp_db)*npoints);
366 double * eqfreqsamples;
367 eqfreqsamples = malloc(sizeof(*eqfreqsamples)*npoints);
369 eqfreqresp = malloc(sizeof(*eqfreqresp)*npoints);
371 ir = malloc(sizeof(*ir)*npoints*2);
374 // And what to do for floating point sample coding?
376 switch(sfinfo.format & SF_FORMAT_SUBMASK) {
377 // all signed bitdepth
391 printf("No known bitdepth! Exiting ...\n");
398 equalinterval2(freqsamples, freqresp_db, eqfreqsamples, eqfreqresp_db, npoints, samplingfreq, origpoints, bitdepth);
399 convloglin(eqfreqresp_db, eqfreqresp, npoints);
402 for (int i=0; i < npoints; i++) {
403 printf("%d\t%.2f\t%.2f\t%.2f\n", i, eqfreqsamples[i], eqfreqresp_db[i], eqfreqresp[i]);
407 inversefft2(eqfreqresp, ir, npoints);
409 // read through the entire file
412 totsum = malloc(sizeof(struct Sum));
415 totsum->nsamples = 0;
417 totsum->mean = 0.0; // Do I write anything here?
420 sf_count_t samples_read = 0;
422 // Main loop through audio file
425 pthread_t tid[numCPU];
426 struct WorkerArgs ** WorkerArgsArray;
427 WorkerArgsArray = malloc(sizeof(struct WorkerArgs *)*numCPU);
428 int staindex = 0; //shorttermarrayindex
431 while((samples_read = sf_read_double(file, buffer, buffersizesamples)) > 0) {
433 WorkerArgsArray[worker_id]=malloc(sizeof(struct WorkerArgs));
434 WorkerArgsArray[worker_id]->nsamples = samples_read;
435 WorkerArgsArray[worker_id]->nch = sfinfo.channels;
436 WorkerArgsArray[worker_id]->npoints=npoints;
437 WorkerArgsArray[worker_id]->ir = ir;
438 WorkerArgsArray[worker_id]->ptrtotsum = totsum;
440 WorkerArgsArray[worker_id]->chconf = channelconfcalvector;
442 WorkerArgsArray[worker_id]->shorttermindex = staindex++;
443 WorkerArgsArray[worker_id]->leqm10flag = 1;
444 WorkerArgsArray[worker_id]->shorttermarray = shorttermaveragedarray;
446 WorkerArgsArray[worker_id]->shorttermindex = 0;
447 WorkerArgsArray[worker_id]->leqm10flag = 0;
450 WorkerArgsArray[worker_id]->argbuffer = malloc(sizeof(double)*buffersizesamples);
451 memcpy(WorkerArgsArray[worker_id]->argbuffer, buffer, samples_read*sizeof(double));
455 pthread_attr_init(&attr);
456 pthread_create(&tid[worker_id], &attr, worker_function, WorkerArgsArray[worker_id]);
459 if (worker_id == numCPU) {
461 //maybe here wait for all cores to output before going on
462 for (int idxcpu = 0; idxcpu < numCPU; idxcpu++) {
463 pthread_join(tid[idxcpu], NULL);
464 free(WorkerArgsArray[idxcpu]);
465 WorkerArgsArray[idxcpu] = NULL;
467 //simply log here your measurement it will be a multiple of your threads and your buffer
469 meanoverduration(totsum); //update leq(m) until now and log it
470 logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
476 //end while worker_id
477 /// End looping cores
478 } // main loop through file
480 //here I should wait for rest worker (< numcpu)
481 //but I need to dispose of thread id.
482 if (worker_id != 0) { // worker_id = 0 means the number of samples was divisible through the number of cpus
483 for (int idxcpu = 0; idxcpu < worker_id; idxcpu++) { //worker_id is at this point one unit more than threads launched
484 pthread_join(tid[idxcpu], NULL);
485 free(WorkerArgsArray[idxcpu]);
486 WorkerArgsArray[idxcpu] = NULL;
488 //also log here for a last value
490 meanoverduration(totsum); //update leq(m) until now and log it
491 logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
494 // mean of scalar sum over duration
496 meanoverduration(totsum);
498 printf("Leq(noW): %.4f\n", totsum->rms); // Leq(no Weighting)
500 printf("Leq(M): %.4f\n", totsum->leqm);
503 struct timespec stoptime;
504 long stoptimenanoseconds;
505 long executionnanoseconds;
506 clock_gettime(CLOCK_MONOTONIC, &stoptime);
508 if (stoptime.tv_nsec < starttime.tv_nsec) {
509 stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
511 stoptimenanoseconds = stoptime.tv_nsec;
513 executionnanoseconds = stoptimenanoseconds - starttime.tv_nsec;
514 printf("Total execution time is %.6f seconds\n", ((double) stoptime.tv_sec) - ((double) starttime.tv_sec) + ((double) executionnanoseconds / 1000000000.00));
520 //Take the array with the short term accumulators
521 double interval = 10.0;
522 //create a rolling average according to rolling interval
523 int rollint; // in short 10*60 = 600 sec 600/0.850
526 //how many element of the array to consider for the rollint?
527 //that is how many buffersizems in the interval - interval could be parameterized(?)
528 double tempint = 60.0 * interval / (((double) buffersizems) /1000.0);
529 rollint = (int) tempint;
530 //dispose of the rest
531 if (tempint - ((double) rollint) > 0) {
537 while(indexlong < (numbershortperiods - rollint)) {
539 double accumulator = 0;
541 double averagedaccumulator = 0;
542 for (int indexshort = 0; indexshort < rollint; indexshort++) {
544 accumulator += shorttermaveragedarray[indexshort+indexlong];
545 } //end internal loop
546 averagedaccumulator = accumulator/((double) rollint);
547 logleqm10(leqm10logfile, ((double) (indexlong+rollint)) * ((double) buffersizems / 1000.0), averagedaccumulator);
549 } //end external loop
551 fclose(leqm10logfile);
552 free(shorttermaveragedarray);
553 shorttermaveragedarray = NULL;
565 eqfreqsamples = NULL;
572 free(channelconfcalvector);
573 channelconfcalvector = NULL;
574 free(WorkerArgsArray);
575 WorkerArgsArray = NULL;
586 void * worker_function(void * argstruct) {
588 struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
590 double * sumandsquarebuffer;
591 double * csumandsquarebuffer;
592 double * chsumaccumulator_norm;
593 double * chsumaccumulator_conv;
596 sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
599 csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
601 chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
603 chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
607 for (int i = 0; i < thisWorkerArgs->nsamples / thisWorkerArgs->nch; i++) {
608 sumandsquarebuffer[i] = 0.0;
609 csumandsquarebuffer[i] = 0.0;
610 chsumaccumulator_norm[i] = 0.0;
611 chsumaccumulator_conv[i] = 0.0;
616 for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
618 double * normalizedbuffer;
619 double * convolvedbuffer;
622 normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
624 convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
627 for (int n=ch, m= 0; n < thisWorkerArgs->nsamples; n += thisWorkerArgs->nch, m++) {
628 // 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
629 //so not normalized but calibrated
630 normalizedbuffer[m] = thisWorkerArgs->argbuffer[n]*thisWorkerArgs->chconf[ch]; //this scale amplitude according to specified calibration
636 convolv_buff(normalizedbuffer, convolvedbuffer, thisWorkerArgs->ir, thisWorkerArgs->nsamples / thisWorkerArgs->nch, thisWorkerArgs->npoints * 2);
637 //rectify, square und sum
638 rectify(csumandsquarebuffer,convolvedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
639 rectify(sumandsquarebuffer,normalizedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
641 accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
642 accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
645 free(normalizedbuffer);
646 normalizedbuffer= NULL;
648 free(convolvedbuffer);
649 convolvedbuffer=NULL;
651 } // loop through channels
653 //Create a function for this also a tag so that the worker know if he has to do this or not
655 if (thisWorkerArgs->leqm10flag) {
656 thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
658 printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
661 pthread_mutex_lock(&mutex);
662 // this should be done under mutex conditions -> shared resources!
663 sumsamples(thisWorkerArgs->ptrtotsum, chsumaccumulator_norm, chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
664 pthread_mutex_unlock(&mutex);
667 free(sumandsquarebuffer);
668 sumandsquarebuffer=NULL;
670 free(csumandsquarebuffer);
671 csumandsquarebuffer=NULL;
673 free(chsumaccumulator_norm);
674 chsumaccumulator_norm=NULL;
676 free(chsumaccumulator_conv);
677 chsumaccumulator_conv=NULL;
679 free(thisWorkerArgs->argbuffer);
680 thisWorkerArgs->argbuffer = NULL;
681 // the memory pointed to by this pointer is freed in main
682 // it is the same memory for all worker
683 // but it is necessary to set pointer to NULL otherwise free will not work later (really?)
684 thisWorkerArgs->chconf = NULL;
690 //to get impulse response frequency response at equally spaced intervals is needed
692 int equalinterval( double * freqsamples, double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
696 double pass = ((double) (samplingfreq >> 1)) / ((double) points);
697 for (int ieq = 0, i = 0; ieq < points; ieq++) {
699 eqfreqsamples[ieq] = freq;
701 if ((freq == 0.0) || (freq < freqsamples[1])) {
702 eqfreqresp[ieq] = freqresp[0];
706 if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
707 eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp[i];
708 } else if (freq >=freqsamples[i+1]) {
709 while(freq >= freqsamples[i+1]) {
711 if ((i + 1) >= origpoints) {
715 if ((i+1) < origpoints) {
716 eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
718 eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
730 //the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
731 //it is also different for the way it interpolates between DC and 31 Hz
732 // Pay attention that also arguments to the functions are changed
733 int equalinterval2( double freqsamples[], double freqresp_db[], double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile) {
737 //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
738 double dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
739 //double dcatt = -90.3;
740 double pass = ((double) (samplingfreq >> 1)) / ((double) points);
741 for (int ieq = 0, i = 0; ieq < points; ieq++) {
743 eqfreqsamples[ieq] = freq;
746 eqfreqresp[ieq] = dcatt;
747 } else if (freq < freqsamples[0]) { // this has a lot of influence on final Leq(M) value
748 eqfreqresp[ieq] = ((freqresp_db[0] - dcatt) / (freqsamples[0] - 0)) * freq + dcatt;
749 //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
753 if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
754 eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp_db[i];
755 } else if (freq >=freqsamples[i+1]) {
756 while(freq >= freqsamples[i+1]) {
758 if ((i + 1) >= origpoints) {
762 if ((i+1) < origpoints) {
763 eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
765 eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
778 int convloglin(double * in, double * out, int points) {
779 for (int i = 0; i < points; i++) {
780 out[i] = powf(10, (in[i]/20.0));
787 double convlinlog_single(double in) {
794 double convloglin_single(double in) {
796 out = powf(10, in/20.0f);
802 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
806 for (int i = 0; i < sigin_dim; i++) {
809 for (int l = impresp_dim - 1; l >=0; l--,m++) {
810 if (m >= sigin_dim) {
813 sum += sigin[m]*impresp[l];
823 void inversefft2(double * eqfreqresp, double * ir, int npoints) {
824 for (int n = 0; n < npoints; n++) {
826 double partial = 0.0;
828 for (int m = 1; m <= npoints -1; m++) {
829 partial = cos(2.0*M_PI*((double) m)*( ( ((double) n) - ( ((double) npoints) * 2.0 -1 ) / 2 ) / ( ((double) npoints) * 2.0) ));
830 parsum = parsum + eqfreqresp[m]*partial;
832 ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
834 printf("%.4f\n", ir[n]);
837 for (int n = 0; n < npoints; n++) {
838 ir[npoints+n] = ir[npoints-(n + 1)];
840 printf("%.4f\n", ir[npoints+n]);
847 // scale input according to required calibration
848 // this could be different for certain digital cinema formats
849 double inputcalib (double dbdiffch) {
851 double coeff = pow(10, dbdiffch/20);
856 //rectify, square and sum
857 int rectify(double * squared, double * inputsamples, int nsamples){
858 for (int i = 0; i < nsamples; i++) {
859 squared[i] = (double) powf(inputsamples[i], 2);
865 int initbuffer(double * buffertoinit, int nsamples) {
866 for (int i = 0; i < nsamples; i++) {
867 buffertoinit[i] = 0.0;
873 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
874 for (int i = 0; i < nsamples; i++) {
875 chaccumulator[i] += inputchannel[i];
880 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples) {
881 ts->nsamples += nsamples;
882 for (int i=0; i < nsamples; i++) {
883 ts->sum += inputsamples[i];
884 ts->csum += cinputsamples[i];
890 int meanoverduration(struct Sum * oldsum) {
891 oldsum->mean = pow(oldsum->sum / ((double) oldsum->nsamples), 0.500);
892 oldsum->cmean = pow(oldsum->csum / ((double) oldsum->nsamples), 0.500);
893 oldsum->rms = 20*log10(oldsum->mean);
894 oldsum->leqm = 20*log10(oldsum->cmean) + 110.600;//
895 // 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.
896 //this value is obtained calibrating with a -20 dBFS Dolby Tone (RMS) I think this is correct
897 //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
898 //This is only approximate as you should use a separate calibration according to the Dolby Format. Also for SW the tone should be 100Hz (?)
903 double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
905 for (int i=0; i < nsamples; i++) {
906 stsum += channelaccumulator[i];
909 return stsum / (double) nsamples;
912 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
914 fprintf(filehandle, "%.4f", featuretimesec);
915 fprintf(filehandle, "\t");
916 fprintf(filehandle, "%.4f\n", oldsum->leqm);
921 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
922 double leqm10 = 20*log10(pow(longaverage, 0.500)) + 110.600;
923 fprintf(filehandle, "%.4f", featuretimesec);
924 fprintf(filehandle, "\t");
925 fprintf(filehandle, "%.4f\n", leqm10);