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-2018 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.18 (C) Luca Trisciani 2011-2013, 2017-2018
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-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", 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[1024];
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-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";
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);
216 if (strcmp(argv[in], "-version") == 0) {
218 printf("leqm-nrt version 0.18\n");
222 if (strcmp(argv[in], "-numcpus") == 0) {
223 numCPU= atoi(argv[in + 1]);
225 printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", numCPU);
229 if (strcmp(argv[in], "-timing") == 0) {
232 printf("Execution time will be measured.\n");
237 if (strcmp(argv[in], "-logleqm10") == 0) {
240 printf("Leq(M)10 data will be logged to the file leqm10.txt\n");
244 if (strcmp(argv[in], "-logleqm") == 0) {
247 printf("Leq(M) data will be logged to the file leqmlog.txt\n");
252 if (strcmp(argv[in], "-leqnw") == 0) {
255 printf("Leq(nW) - unweighted - will be outputted.\n");
260 if (strcmp(argv[in], "-buffersize") == 0) {
261 buffersizems = atoi(argv[in + 1]);
263 printf("Buffersize will be set to %d milliseconds.\n", buffersizems);
268 if (parameterstate==0) {
274 //postprocessing parameters
275 if (numcalread == sfinfo.channels) {
276 for (int cind = 0; cind < sfinfo.channels; cind++) {
277 channelconfcalvector[cind] = convloglin_single(tempchcal[cind]);
281 else if ((numcalread == 0) && (sfinfo.channels == 6)) {
282 double conf51[6] = {0, 0, 0, 0, -3, -3};
283 for (int cind = 0; cind < sfinfo.channels; cind++) {
284 channelconfcalvector[cind] = convloglin_single(conf51[cind]);
289 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");
291 free(channelconfcalvector);
292 channelconfcalvector = NULL;
300 char tempstring[1536];
301 strcpy(tempstring, soundfilename);
302 strcat(tempstring, ".leqm10.txt");
303 leqm10logfile = fopen(tempstring, "w");
304 if (leqm10logfile == NULL) {
305 printf("Could not open file to write log leqm10 data!\n");
313 char tempstring[1536];
314 strcpy(tempstring, soundfilename);
315 strcat(tempstring, ".leqmlog.txt");
316 leqmlogfile = fopen(tempstring, "w");
317 if (leqmlogfile == NULL) {
318 printf("Could not open file to write log leqm data!\n");
325 clock_gettime(CLOCK_MONOTONIC, &starttime);
332 // reading to a double or float buffer with sndfile take care of normalization
334 static double buffer[BUFFER_LEN]; // it seems this must be static. I don't know why
337 // buffer = new double [BUFFER_LEN];
338 //buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
339 if ((sfinfo.samplerate*buffersizems)%1000) {
340 printf("Please fine tune the buffersize according to the sample rate\n");
343 // write a function to do that
347 buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
348 buffer = malloc(sizeof(double)*buffersizesamples);
350 samplingfreq = sfinfo.samplerate;
354 //if duration < 10 mm exit
356 double featdursec = sfinfo.frames / sfinfo.samplerate;
357 if ((featdursec/60.0) < 10.0) {
358 printf("The audio file is too short to measure Leq(m10).\n");
362 //how many short periods in overall duration
363 int remainder = sfinfo.frames % (sfinfo.samplerate*buffersizems/1000);
364 if (remainder == 0) numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000);
365 else numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000) + 1;
368 shorttermaveragedarray = malloc(sizeof(*shorttermaveragedarray)*numbershortperiods);
372 //End opening audio file
376 double freqsamples[] = {31, 63, 100, 200, 400, 800, 1000, 2000, 3150, 4000, 5000, 6300, 7100, 8000, 9000, 10000, 12500, 14000, 16000, 20000, 31500};
377 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};
379 double * eqfreqresp_db;
380 eqfreqresp_db = malloc(sizeof(*eqfreqresp_db)*npoints);
382 double * eqfreqsamples;
383 eqfreqsamples = malloc(sizeof(*eqfreqsamples)*npoints);
385 eqfreqresp = malloc(sizeof(*eqfreqresp)*npoints);
387 ir = malloc(sizeof(*ir)*npoints*2);
390 // And what to do for floating point sample coding?
392 switch(sfinfo.format & SF_FORMAT_SUBMASK) {
393 // all signed bitdepth
407 printf("No known bitdepth! Exiting ...\n");
414 equalinterval2(freqsamples, freqresp_db, eqfreqsamples, eqfreqresp_db, npoints, samplingfreq, origpoints, bitdepth);
415 convloglin(eqfreqresp_db, eqfreqresp, npoints);
418 for (int i=0; i < npoints; i++) {
419 printf("%d\t%.2f\t%.2f\t%.6f\n", i, eqfreqsamples[i], eqfreqresp_db[i], eqfreqresp[i]);
423 inversefft2(eqfreqresp, ir, npoints);
425 // read through the entire file
428 totsum = malloc(sizeof(struct Sum));
431 totsum->nsamples = 0;
433 totsum->mean = 0.0; // Do I write anything here?
436 sf_count_t samples_read = 0;
438 // Main loop through audio file
441 pthread_t tid[numCPU];
442 struct WorkerArgs ** WorkerArgsArray;
443 WorkerArgsArray = malloc(sizeof(struct WorkerArgs *)*numCPU);
444 int staindex = 0; //shorttermarrayindex
447 while((samples_read = sf_read_double(file, buffer, buffersizesamples)) > 0) {
449 WorkerArgsArray[worker_id]=malloc(sizeof(struct WorkerArgs));
450 WorkerArgsArray[worker_id]->nsamples = samples_read;
451 WorkerArgsArray[worker_id]->nch = sfinfo.channels;
452 WorkerArgsArray[worker_id]->npoints=npoints;
453 WorkerArgsArray[worker_id]->ir = ir;
454 WorkerArgsArray[worker_id]->ptrtotsum = totsum;
456 WorkerArgsArray[worker_id]->chconf = channelconfcalvector;
458 WorkerArgsArray[worker_id]->shorttermindex = staindex++;
459 WorkerArgsArray[worker_id]->leqm10flag = 1;
460 WorkerArgsArray[worker_id]->shorttermarray = shorttermaveragedarray;
462 WorkerArgsArray[worker_id]->shorttermindex = 0;
463 WorkerArgsArray[worker_id]->leqm10flag = 0;
466 WorkerArgsArray[worker_id]->argbuffer = malloc(sizeof(double)*buffersizesamples);
467 memcpy(WorkerArgsArray[worker_id]->argbuffer, buffer, samples_read*sizeof(double));
471 pthread_attr_init(&attr);
472 pthread_create(&tid[worker_id], &attr, worker_function, WorkerArgsArray[worker_id]);
475 if (worker_id == numCPU) {
477 //maybe here wait for all cores to output before going on
478 for (int idxcpu = 0; idxcpu < numCPU; idxcpu++) {
479 pthread_join(tid[idxcpu], NULL);
480 free(WorkerArgsArray[idxcpu]->argbuffer);
481 WorkerArgsArray[idxcpu]->argbuffer = NULL;
482 free(WorkerArgsArray[idxcpu]);
483 WorkerArgsArray[idxcpu] = NULL;
485 //simply log here your measurement it will be a multiple of your threads and your buffer
487 meanoverduration(totsum); //update leq(m) until now and log it
488 logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
494 //end while worker_id
495 /// End looping cores
496 } // main loop through file
498 //here I should wait for rest worker (< numcpu)
499 //but I need to dispose of thread id.
500 if (worker_id != 0) { // worker_id = 0 means the number of samples was divisible through the number of cpus
501 for (int idxcpu = 0; idxcpu < worker_id; idxcpu++) { //worker_id is at this point one unit more than threads launched
502 pthread_join(tid[idxcpu], NULL);
503 free(WorkerArgsArray[idxcpu]->argbuffer);
504 WorkerArgsArray[idxcpu]->argbuffer = NULL;
505 free(WorkerArgsArray[idxcpu]);
506 WorkerArgsArray[idxcpu] = NULL;
508 //also log here for a last value
510 meanoverduration(totsum); //update leq(m) until now and log it
511 logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
514 // mean of scalar sum over duration
516 meanoverduration(totsum);
518 printf("Leq(nW): %.4f\n", totsum->rms); // Leq(no Weighting)
520 printf("Leq(M): %.4f\n", totsum->leqm);
523 struct timespec stoptime;
524 long stoptimenanoseconds;
525 long executionnanoseconds;
526 clock_gettime(CLOCK_MONOTONIC, &stoptime);
528 if (stoptime.tv_nsec < starttime.tv_nsec) {
529 stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
531 stoptimenanoseconds = stoptime.tv_nsec;
533 executionnanoseconds = stoptimenanoseconds - starttime.tv_nsec;
534 printf("Total execution time is %.6f seconds\n", ((double) stoptime.tv_sec) - ((double) starttime.tv_sec) + ((double) executionnanoseconds / 1000000000.00));
540 //Take the array with the short term accumulators
541 double interval = 10.0;
542 //create a rolling average according to rolling interval
543 int rollint; // in short 10*60 = 600 sec 600/0.850
546 //how many element of the array to consider for the rollint?
547 //that is how many buffersizems in the interval - interval could be parameterized(?)
548 double tempint = 60.0 * interval / (((double) buffersizems) /1000.0);
549 rollint = (int) tempint;
550 //dispose of the rest
551 if (tempint - ((double) rollint) > 0) {
557 while(indexlong < (numbershortperiods - rollint)) {
559 double accumulator = 0;
561 double averagedaccumulator = 0;
562 for (int indexshort = 0; indexshort < rollint; indexshort++) {
564 accumulator += shorttermaveragedarray[indexshort+indexlong];
565 } //end internal loop
566 averagedaccumulator = accumulator/((double) rollint);
567 logleqm10(leqm10logfile, ((double) (indexlong+rollint)) * ((double) buffersizems / 1000.0), averagedaccumulator);
569 } //end external loop
571 fclose(leqm10logfile);
572 free(shorttermaveragedarray);
573 shorttermaveragedarray = NULL;
585 eqfreqsamples = NULL;
592 free(channelconfcalvector);
593 channelconfcalvector = NULL;
594 free(WorkerArgsArray);
595 WorkerArgsArray = NULL;
606 void * worker_function(void * argstruct) {
608 struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
610 double * sumandsquarebuffer;
611 double * csumandsquarebuffer;
612 double * chsumaccumulator_norm;
613 double * chsumaccumulator_conv;
616 sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
619 csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
621 chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
623 chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
627 for (int i = 0; i < thisWorkerArgs->nsamples / thisWorkerArgs->nch; i++) {
628 sumandsquarebuffer[i] = 0.0;
629 csumandsquarebuffer[i] = 0.0;
630 chsumaccumulator_norm[i] = 0.0;
631 chsumaccumulator_conv[i] = 0.0;
636 for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
638 double * normalizedbuffer;
639 double * convolvedbuffer;
642 normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
644 convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
647 for (int n=ch, m= 0; n < thisWorkerArgs->nsamples; n += thisWorkerArgs->nch, m++) {
648 // 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
649 //so not normalized but calibrated
650 normalizedbuffer[m] = thisWorkerArgs->argbuffer[n]*thisWorkerArgs->chconf[ch]; //this scale amplitude according to specified calibration
656 convolv_buff(normalizedbuffer, convolvedbuffer, thisWorkerArgs->ir, thisWorkerArgs->nsamples / thisWorkerArgs->nch, thisWorkerArgs->npoints * 2);
657 //rectify, square und sum
658 rectify(csumandsquarebuffer,convolvedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
659 rectify(sumandsquarebuffer,normalizedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
661 accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
662 accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
665 free(normalizedbuffer);
666 normalizedbuffer= NULL;
668 free(convolvedbuffer);
669 convolvedbuffer=NULL;
671 } // loop through channels
673 //Create a function for this also a tag so that the worker know if he has to do this or not
675 if (thisWorkerArgs->leqm10flag) {
676 thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
678 printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
681 pthread_mutex_lock(&mutex);
682 // this should be done under mutex conditions -> shared resources!
683 sumsamples(thisWorkerArgs->ptrtotsum, chsumaccumulator_norm, chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
684 pthread_mutex_unlock(&mutex);
687 free(sumandsquarebuffer);
688 sumandsquarebuffer=NULL;
690 free(csumandsquarebuffer);
691 csumandsquarebuffer=NULL;
693 free(chsumaccumulator_norm);
694 chsumaccumulator_norm=NULL;
696 free(chsumaccumulator_conv);
697 chsumaccumulator_conv=NULL;
699 free(thisWorkerArgs->argbuffer);
700 thisWorkerArgs->argbuffer = NULL;
701 // the memory pointed to by this pointer is freed in main
702 // it is the same memory for all worker
703 // but it is necessary to set pointer to NULL otherwise free will not work later (really?)
704 thisWorkerArgs->chconf = NULL;
710 //to get impulse response frequency response at equally spaced intervals is needed
712 int equalinterval( double * freqsamples, double * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
716 double pass = ((double) (samplingfreq >> 1)) / ((double) points);
717 for (int ieq = 0, i = 0; ieq < points; ieq++) {
719 eqfreqsamples[ieq] = freq;
721 if ((freq == 0.0) || (freq < freqsamples[1])) {
722 eqfreqresp[ieq] = freqresp[0];
726 if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
727 eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp[i];
728 } else if (freq >=freqsamples[i+1]) {
729 while(freq >= freqsamples[i+1]) {
731 if ((i + 1) >= origpoints) {
735 if ((i+1) < origpoints) {
736 eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
738 eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
750 //the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
751 //it is also different for the way it interpolates between DC and 31 Hz
752 // Pay attention that also arguments to the functions are changed
753 int equalinterval2( double freqsamples[], double freqresp_db[], double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile) {
757 //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
758 double dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
759 //double dcatt = -90.3;
760 double pass = ((double) (samplingfreq >> 1)) / ((double) points);
761 for (int ieq = 0, i = 0; ieq < points; ieq++) {
763 eqfreqsamples[ieq] = freq;
766 eqfreqresp[ieq] = dcatt;
767 } else if (freq < freqsamples[0]) { // this has a lot of influence on final Leq(M) value
768 eqfreqresp[ieq] = ((freqresp_db[0] - dcatt) / (freqsamples[0] - 0)) * freq + dcatt;
769 //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
773 if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
774 eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp_db[i];
775 } else if (freq >=freqsamples[i+1]) {
776 while(freq >= freqsamples[i+1]) {
778 if ((i + 1) >= origpoints) {
782 if ((i+1) < origpoints) {
783 eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
785 eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
798 int convloglin(double * in, double * out, int points) {
799 for (int i = 0; i < points; i++) {
800 out[i] = powf(10, (in[i]/20.0));
807 double convlinlog_single(double in) {
814 double convloglin_single(double in) {
816 out = powf(10, in/20.0f);
822 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
826 for (int i = 0; i < sigin_dim; i++) {
829 for (int l = impresp_dim - 1; l >=0; l--,m++) {
830 if (m >= sigin_dim) {
833 sum += sigin[m]*impresp[l];
843 void inversefft2(double * eqfreqresp, double * ir, int npoints) {
844 for (int n = 0; n < npoints; n++) {
846 double partial = 0.0;
848 for (int m = 1; m <= npoints -1; m++) {
849 partial = cos(2.0*M_PI*((double) m)*( ( ((double) n) - ( ((double) npoints) * 2.0 -1 ) / 2 ) / ( ((double) npoints) * 2.0) ));
850 parsum = parsum + eqfreqresp[m]*partial;
852 ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
854 printf("%.4f\n", ir[n]);
857 for (int n = 0; n < npoints; n++) {
858 ir[npoints+n] = ir[npoints-(n + 1)];
860 printf("%.4f\n", ir[npoints+n]);
867 // scale input according to required calibration
868 // this could be different for certain digital cinema formats
869 double inputcalib (double dbdiffch) {
871 double coeff = pow(10, dbdiffch/20);
876 //rectify, square and sum
877 int rectify(double * squared, double * inputsamples, int nsamples){
878 for (int i = 0; i < nsamples; i++) {
879 squared[i] = (double) powf(inputsamples[i], 2);
885 int initbuffer(double * buffertoinit, int nsamples) {
886 for (int i = 0; i < nsamples; i++) {
887 buffertoinit[i] = 0.0;
893 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
894 for (int i = 0; i < nsamples; i++) {
895 chaccumulator[i] += inputchannel[i];
900 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples) {
901 ts->nsamples += nsamples;
902 for (int i=0; i < nsamples; i++) {
903 ts->sum += inputsamples[i];
904 ts->csum += cinputsamples[i];
910 int meanoverduration(struct Sum * oldsum) {
911 oldsum->mean = pow(oldsum->sum / ((double) oldsum->nsamples), 0.500);
912 oldsum->cmean = pow(oldsum->csum / ((double) oldsum->nsamples), 0.500);
913 oldsum->rms = 20*log10(oldsum->mean) + 108.010299957;
914 oldsum->leqm = 20*log10(oldsum->cmean) + 108.010299957;//
917 How the final offset is calculated without reference to a test tone:
918 P0 is the SPL reference 20 uPa
919 Reference SPL is RMS ! So 85 SPL over 20 uPa is 10^4.25 x 0.000020 = 0.355655882 Pa (RMS),
920 but Peak value is 0.355655882 x sqr(2) = 0.502973372 that is 20 x log ( 0.502973372 / 0.000020) = 88.010299957
921 To that one has to add the 20 dB offset of the reference -20dBFS: 88.010299957 + 20.00 = 108.010299957
923 /*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
928 double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
930 for (int i=0; i < nsamples; i++) {
931 stsum += channelaccumulator[i];
934 return stsum / (double) nsamples;
937 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
939 fprintf(filehandle, "%.4f", featuretimesec);
940 fprintf(filehandle, "\t");
941 fprintf(filehandle, "%.4f\n", oldsum->leqm);
946 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
947 double leqm10 = 20*log10(pow(longaverage, 0.500)) + 108.010299957;
948 fprintf(filehandle, "%.4f", featuretimesec);
949 fprintf(filehandle, "\t");
950 fprintf(filehandle, "%.4f\n", leqm10);