src/leqm-nrt.c

   1 /*
   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"
   6
   7     Copyright (C) 2011-2013, 2017-2018 Luca Trisciani
   8
   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.
  13
  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.
  18
  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/>.
  21
  22  */
  23
  24
  25
  26 #include <stdio.h>
  27 #include <math.h>
  28 #include <sndfile.h>
  29 #include <unistd.h>
  30 #include <pthread.h>
  31 #include <string.h>
  32 #include <stdlib.h>
  33 #include <time.h>
  34 #include <ctype.h>
  35 #include <iso646.h>
  36
  37 #ifdef _WIN32
  38 #include <windows.h>
  39 #elif __APPLE__
  40 #include <sys/param.h>
  41 #include <sys/sysctl.h
  42 #endif
  43
  44
  45 // Version 0.0.18 (C) Luca Trisciani 2011-2013, 2017-2018
  46 // Tool from the DCP-Werkstatt Software Bundle
  47
  48
  49
  50 // COMPILATION
  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
  53
  54 //#define DEBUG
  55
  56
  57 struct Sum {
  58   double csum; // convolved sum
  59   double sum; // flat sum
  60     int nsamples;
  61   double cmean; //convolved mean
  62     double mean;
  63     double leqm;
  64   double rms;
  65 };
  66
  67 struct WorkerArgs {
  68   double * argbuffer;
  69   int nsamples;
  70   int nch;
  71   int npoints;
  72   double * ir;
  73   struct Sum * ptrtotsum;
  74   double * chconf;
  75   int shorttermindex;
  76   double * shorttermarray;
  77   int leqm10flag;
  78 };
  79
  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);
  97
  98
  99 pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
 100
 101
 102 int main(int argc, const char ** argv)
 103 {
 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
 107   int bitdepth;
 108   // double normalizer;
 109   int timing = 0;
 110   struct timespec starttime;
 111   int fileopenstate = 0;
 112   int leqm10 = 0;
 113   int leqmlog = 0;
 114   #if defined __unix__ || defined  __APPLE__
 115   int numCPU = sysconf(_SC_NPROCESSORS_ONLN) - 1;
 116   #elif defined _WIN64 || defined _WIN32
 117   SYSTEM_INFO sysinfo;
 118   GetSystemInfo(&sysinfo);
 119   int numCPU = sysinfo.dwNumberOfProcessors - 1;
 120   #endif
 121
 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;
 126   SNDFILE *file;
 127   file=NULL;
 128   SF_INFO sfinfo;
 129   FILE *leqm10logfile;
 130   leqm10logfile = NULL;
 131   FILE *leqmlogfile;
 132   leqmlogfile = 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];
 136         int numcalread = 0;
 137         double * shorttermaveragedarray;
 138         shorttermaveragedarray = NULL;
 139         int numbershortperiods = 0;
 140         int parameterstate = 0;
 141         int leqnw = 0;
 142
 143         char soundfilename[512];
 144         // This is a requirement of sndfile library, do not forget it.
 145
 146         memset(&sfinfo, 0, sizeof(sfinfo));
 147
 148
 149   if (argc == 1)
 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";
 151       printf(helptext);
 152       printf("Please indicate a sound file to be processed.\n");
 153       return 0;
 154   }
 155
 156
 157     for (int in = 1; in < argc;) {
 158
 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));
 164             return 1;
 165           }
 166
 167           strcpy(soundfilename, argv[in]);
 168              fileopenstate = 1;
 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);
 175              in++;
 176              continue;
 177         } else {
 178           free(channelconfcalvector);
 179           channelconfcalvector = NULL;
 180           return 0;
 181         }
 182       }
 183
 184
 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.
 192         */
 193
 194         in++;
 195         for (;;)  {
 196         if (in < argc) {
 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++]);
 200           } else break;
 201         } else break;
 202
 203         } //for
 204         continue;
 205       }
 206
 207          if (strcmp(argv[in], "-convpoints") == 0) {
 208              npoints = atoi(argv[in + 1]);
 209              in+=2;
 210              printf("Convolution points sets to %d.\n", npoints);
 211              continue;
 212
 213       }
 214
 215
 216                 if (strcmp(argv[in], "-version") == 0) {
 217              in++;
 218              printf("leqm-nrt version 0.18\n");
 219              continue;
 220
 221       }
 222               if (strcmp(argv[in], "-numcpus") == 0) {
 223                 numCPU= atoi(argv[in + 1]);
 224              in+=2;
 225              printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", numCPU);
 226              continue;
 227
 228       }
 229               if (strcmp(argv[in], "-timing") == 0) {
 230                 timing = 1;
 231              in++;
 232              printf("Execution time will be measured.\n");
 233              continue;
 234
 235       }
 236
 237                       if (strcmp(argv[in], "-logleqm10") == 0) {
 238                 leqm10 = 1;
 239              in++;
 240              printf("Leq(M)10 data will be logged to the file leqm10.txt\n");
 241              continue;
 242
 243       }
 244                                       if (strcmp(argv[in], "-logleqm") == 0) {
 245                 leqmlog = 1;
 246              in++;
 247              printf("Leq(M) data will be logged to the file leqmlog.txt\n");
 248              continue;
 249
 250       }
 251
 252                                       if (strcmp(argv[in], "-leqnw") == 0) {
 253                 leqnw = 1;
 254              in++;
 255              printf("Leq(nW) - unweighted -  will be outputted.\n");
 256              continue;
 257
 258       }
 259
 260                                         if (strcmp(argv[in], "-buffersize") == 0) {
 261                 buffersizems = atoi(argv[in + 1]);
 262              in+=2;
 263              printf("Buffersize will be set to %d milliseconds.\n", buffersizems);
 264              continue;
 265
 266       }
 267
 268                                         if (parameterstate==0) {
 269                                           break;
 270                                         }
 271     }
 272 // Open audio file
 273
 274 //postprocessing parameters
 275     if (numcalread == sfinfo.channels) {
 276       for (int cind = 0; cind < sfinfo.channels; cind++) {
 277         channelconfcalvector[cind] = convloglin_single(tempchcal[cind]);
 278
 279       }
 280     }
 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]);
 285         }
 286
 287     } else {
 288
 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");
 290
 291       free(channelconfcalvector);
 292       channelconfcalvector = NULL;
 293       return 0;
 294     }
 295
 296
 297
 298
 299     if (leqm10) {
 300       char tempstring[128];
 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");
 306       }
 307     }
 308
 309
 310
 311
 312     if (leqmlog) {
 313       char tempstring[128];
 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");
 319       }
 320     }
 321
 322
 323   if (timing) {
 324
 325   clock_gettime(CLOCK_MONOTONIC, &starttime);
 326   }
 327
 328
 329
 330
 331
 332   // reading to a double or float buffer with sndfile take care of normalization
 333  /*
 334  static double  buffer[BUFFER_LEN]; // it seems this must be static. I don't know why
 335  */
 336  double * buffer;
 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");
 341    //close file
 342    // free memory
 343    // write a function to do that
 344    return 1;
 345  }
 346
 347   buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
 348   buffer = malloc(sizeof(double)*buffersizesamples);
 349
 350  samplingfreq = sfinfo.samplerate;
 351
 352  if(leqm10) {
 353
 354    //if duration < 10 mm exit
 355
 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");
 359      return 0;
 360    }
 361
 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;
 366
 367    //allocate array
 368    shorttermaveragedarray = malloc(sizeof(*shorttermaveragedarray)*numbershortperiods);
 369  }
 370
 371
 372  //End opening audio file
 373
 374   //ISO 21727:2004(E)
 375   // M Weighting
 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};
 378
 379   double * eqfreqresp_db;
 380   eqfreqresp_db = malloc(sizeof(*eqfreqresp_db)*npoints);
 381
 382   double * eqfreqsamples;
 383   eqfreqsamples = malloc(sizeof(*eqfreqsamples)*npoints);
 384   double * eqfreqresp;
 385   eqfreqresp = malloc(sizeof(*eqfreqresp)*npoints);
 386   double * ir;
 387   ir = malloc(sizeof(*ir)*npoints*2);
 388
 389
 390 // And what to do for floating point sample coding?
 391
 392    switch(sfinfo.format & SF_FORMAT_SUBMASK) {
 393    // all signed bitdepth
 394  case 0x0001:
 395    bitdepth = 8;
 396    break;
 397  case 0x0002:
 398    bitdepth = 16;
 399    break;
 400  case 0x0003:
 401    bitdepth = 24;
 402    break;
 403  case 0x0004:
 404    bitdepth = 32;
 405    break;
 406  default:
 407    printf("No known bitdepth! Exiting ...\n");
 408    return -1;
 409    }
 410
 411
 412
 413
 414    equalinterval2(freqsamples, freqresp_db, eqfreqsamples, eqfreqresp_db, npoints, samplingfreq, origpoints, bitdepth);
 415   convloglin(eqfreqresp_db, eqfreqresp, npoints);
 416
 417     #ifdef DEBUG
 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]);
 420     }
 421     #endif
 422
 423     inversefft2(eqfreqresp, ir, npoints);
 424
 425 // read through the entire file
 426
 427    struct Sum * totsum;
 428     totsum = malloc(sizeof(struct Sum));
 429     totsum->csum = 0.0;
 430     totsum->sum = 0.0;
 431     totsum->nsamples = 0;
 432     totsum->cmean = 0.0;
 433     totsum->mean = 0.0; // Do I write anything here?
 434     totsum->leqm = 0.0;
 435     totsum->rms = 0.0;
 436     sf_count_t samples_read = 0;
 437
 438  // Main loop through audio file
 439
 440  int worker_id = 0;
 441  pthread_t tid[numCPU];
 442  struct WorkerArgs ** WorkerArgsArray;
 443  WorkerArgsArray = malloc(sizeof(struct WorkerArgs *)*numCPU);
 444  int staindex = 0; //shorttermarrayindex
 445
 446
 447  while((samples_read = sf_read_double(file, buffer, buffersizesamples)) > 0) {
 448
 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;
 455
 456    WorkerArgsArray[worker_id]->chconf = channelconfcalvector;
 457    if (leqm10) {
 458    WorkerArgsArray[worker_id]->shorttermindex = staindex++;
 459    WorkerArgsArray[worker_id]->leqm10flag = 1;
 460    WorkerArgsArray[worker_id]->shorttermarray = shorttermaveragedarray;
 461    } else {
 462      WorkerArgsArray[worker_id]->shorttermindex = 0;
 463      WorkerArgsArray[worker_id]->leqm10flag = 0;
 464    }
 465
 466    WorkerArgsArray[worker_id]->argbuffer = malloc(sizeof(double)*buffersizesamples);
 467    memcpy(WorkerArgsArray[worker_id]->argbuffer, buffer, samples_read*sizeof(double));
 468
 469
 470    pthread_attr_t attr;
 471    pthread_attr_init(&attr);
 472    pthread_create(&tid[worker_id], &attr, worker_function, WorkerArgsArray[worker_id]);
 473
 474    worker_id++;
 475    if (worker_id == numCPU) {
 476        worker_id = 0;
 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;
 484        }
 485               //simply log here your measurement it will be a multiple of your threads and your buffer
 486        if (leqmlog) {
 487          meanoverduration(totsum); //update leq(m) until now and log it
 488        logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
 489                } //endlog
 490    }
 491
 492
 493
 494     //end while worker_id
 495  /// End looping cores
 496   } // main loop through file
 497
 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;
 507    }
 508         //also log here for a last value
 509        if (leqmlog) {
 510          meanoverduration(totsum); //update leq(m) until now and log it
 511        logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
 512                } //endlog
 513         }
 514  // mean of scalar sum over duration
 515
 516  meanoverduration(totsum);
 517  if (leqnw) {
 518  printf("Leq(nW): %.4f\n", totsum->rms); // Leq(no Weighting)
 519  }
 520  printf("Leq(M): %.4f\n", totsum->leqm);
 521
 522   if(timing) {
 523    struct timespec stoptime;
 524    long stoptimenanoseconds;
 525    long executionnanoseconds;
 526    clock_gettime(CLOCK_MONOTONIC, &stoptime);
 527
 528    if (stoptime.tv_nsec < starttime.tv_nsec) {
 529      stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
 530    } else {
 531      stoptimenanoseconds = stoptime.tv_nsec;
 532    }
 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));
 535  }
 536
 537
 538  if (leqm10) {
 539
 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
 544
 545
 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) {
 552      rollint += 1;
 553    }
 554    //two loops
 555    //external loop
 556    int indexlong = 0;
 557    while(indexlong < (numbershortperiods - rollint)) {
 558
 559      double accumulator = 0;
 560      //internal loop
 561      double averagedaccumulator = 0;
 562      for (int indexshort = 0; indexshort < rollint; indexshort++) {
 563
 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);
 568      indexlong++;
 569    } //end external loop
 570
 571      fclose(leqm10logfile);
 572    free(shorttermaveragedarray);
 573    shorttermaveragedarray = NULL;
 574  }
 575
 576
 577  if (leqmlog) {
 578
 579    fclose(leqmlogfile);
 580  }
 581
 582  sf_close(file);
 583
 584  free(eqfreqsamples);
 585  eqfreqsamples = NULL;
 586   free(eqfreqresp_db);
 587  eqfreqresp_db=NULL;
 588  free(eqfreqresp);
 589  eqfreqresp = NULL;
 590  free(ir);
 591  ir = NULL;
 592  free(channelconfcalvector);
 593  channelconfcalvector = NULL;
 594  free(WorkerArgsArray);
 595  WorkerArgsArray = NULL;
 596
 597  free(totsum);
 598  totsum = NULL;
 599  free(buffer);
 600  buffer=NULL;
 601    return 0;
 602 }
 603
 604
 605
 606 void * worker_function(void * argstruct) {
 607
 608   struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
 609
 610    double * sumandsquarebuffer;
 611    double * csumandsquarebuffer;
 612   double * chsumaccumulator_norm;
 613   double * chsumaccumulator_conv;
 614
 615
 616   sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 617
 618
 619   csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 620
 621   chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 622
 623   chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 624
 625
 626
 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;
 632   }
 633
 634
 635
 636   for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
 637
 638     double * normalizedbuffer;
 639     double * convolvedbuffer;
 640
 641
 642     normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 643
 644     convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 645
 646
 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
 651
 652
 653  }
 654
 655  //convolution
 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);
 660
 661  accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 662  accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 663
 664
 665  free(normalizedbuffer);
 666  normalizedbuffer= NULL;
 667
 668  free(convolvedbuffer);
 669  convolvedbuffer=NULL;
 670
 671  } // loop through channels
 672
 673     //Create a function for this also a tag so that the worker know if he has to do this or not
 674
 675   if (thisWorkerArgs->leqm10flag) {
 676     thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 677     #ifdef DEBUG
 678     printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
 679     #endif
 680   }
 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);
 685
 686
 687   free(sumandsquarebuffer);
 688   sumandsquarebuffer=NULL;
 689
 690   free(csumandsquarebuffer);
 691   csumandsquarebuffer=NULL;
 692
 693   free(chsumaccumulator_norm);
 694   chsumaccumulator_norm=NULL;
 695
 696   free(chsumaccumulator_conv);
 697   chsumaccumulator_conv=NULL;
 698
 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;
 705  pthread_exit(0);
 706
 707 }
 708
 709
 710 //to get impulse response frequency response at equally spaced intervals is needed
 711
 712 int equalinterval( double * freqsamples, double  * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
 713     double freq;
 714     // int findex = 0;
 715     // int rindex = 0;
 716     double pass = ((double) (samplingfreq >> 1)) / ((double) points);
 717     for (int ieq = 0, i = 0; ieq < points; ieq++) {
 718         freq = ieq*pass;
 719         eqfreqsamples[ieq] = freq;
 720
 721         if ((freq == 0.0) || (freq < freqsamples[1])) {
 722           eqfreqresp[ieq] = freqresp[0];
 723             continue;
 724     } else {
 725
 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]) {
 730                 i++;
 731                 if ((i + 1) >= origpoints) {
 732                   break;
 733                 }
 734             }
 735             if ((i+1) < origpoints) {
 736             eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
 737             } else {
 738               eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
 739             }
 740         }
 741         }
 742     }
 743     return 0;
 744 }
 745
 746
 747
 748
 749
 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) {
 754     double freq;
 755
 756
 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++) {
 762         freq = ieq*pass;
 763         eqfreqsamples[ieq] = freq;
 764
 765         if (freq == 0.0) {
 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
 770             continue;
 771     } else {
 772
 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]) {
 777                 i++;
 778                 if ((i + 1) >= origpoints) {
 779                   break;
 780                 }
 781             }
 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];
 784             } else {
 785               eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
 786             }
 787         }
 788         }
 789     }
 790     return 0;
 791 }
 792
 793
 794
 795
 796
 797
 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));
 801   }
 802
 803   return 0;
 804 }
 805
 806
 807 double convlinlog_single(double in) {
 808   double out;
 809     out = log(in)*20.0f;
 810   return out;
 811 }
 812
 813
 814 double convloglin_single(double in) {
 815   double out;
 816   out = powf(10, in/20.0f);
 817   return out;
 818 }
 819
 820 // convolution
 821
 822 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
 823
 824
 825   double  sum = 0.0;
 826   for (int i = 0; i < sigin_dim; i++) {
 827
 828     int m = i;
 829     for (int l = impresp_dim - 1; l >=0; l--,m++) {
 830       if (m >= sigin_dim) {
 831         m -= sigin_dim;
 832       }
 833       sum += sigin[m]*impresp[l];
 834     }
 835     sigout[i] = sum;
 836     sum=0.0;
 837     }
 838   return 0;
 839
 840 }
 841
 842
 843 void  inversefft2(double * eqfreqresp, double * ir, int npoints) {
 844   for (int n = 0; n < npoints; n++) {
 845     double parsum = 0.0;
 846     double partial = 0.0;
 847
 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;
 851     }
 852     ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
 853     #ifdef DEBUG
 854     printf("%.4f\n", ir[n]);
 855     #endif
 856   }
 857   for (int n = 0; n < npoints; n++) {
 858     ir[npoints+n] = ir[npoints-(n + 1)];
 859     #ifdef DEBUG
 860     printf("%.4f\n", ir[npoints+n]);
 861     #endif
 862   }
 863
 864
 865 }
 866
 867 // scale input according to required calibration
 868 // this could be different for certain digital cinema formats
 869 double inputcalib (double dbdiffch) {
 870
 871     double coeff = pow(10, dbdiffch/20);
 872     return coeff;
 873
 874 }
 875
 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);
 880     }
 881     return 0;
 882
 883 }
 884
 885 int initbuffer(double * buffertoinit, int nsamples) {
 886   for (int i = 0; i < nsamples; i++) {
 887     buffertoinit[i] = 0.0;
 888
 889   }
 890   return 0;
 891 }
 892
 893 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
 894   for (int i = 0; i < nsamples; i++) {
 895     chaccumulator[i] += inputchannel[i];
 896   }
 897   return 0;
 898 }
 899
 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];
 905   }
 906   return 0;
 907
 908 }
 909
 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.0851;
 914    oldsum->leqm = 20*log10(oldsum->cmean) +  108.0851;//
 915      //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.
 916    //this value is obtained calibrating with a -20 dBFS.
 917
 918 return 0;
 919 }
 920
 921 double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
 922   double stsum = 0.0;
 923   for (int i=0; i < nsamples; i++) {
 924     stsum += channelaccumulator[i];
 925
 926   }
 927   return stsum / (double) nsamples;
 928 }
 929
 930 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
 931
 932   fprintf(filehandle, "%.4f", featuretimesec);
 933   fprintf(filehandle, "\t");
 934   fprintf(filehandle, "%.4f\n", oldsum->leqm);
 935
 936
 937 }
 938
 939 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
 940   double leqm10 = 20*log10(pow(longaverage, 0.500)) +  108.0851;
 941   fprintf(filehandle, "%.4f", featuretimesec);
 942   fprintf(filehandle, "\t");
 943   fprintf(filehandle, "%.4f\n", leqm10);
 944
 945 }