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 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
  36
  37 // Version 0.0.17 (C) Luca Trisciani 2011-2013, 2017
  38 // Tool from the DCP-Werkstatt Software Bundle
  39
  40
  41
  42 // COMPILATION
  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
  45
  46 //#define DEBUG
  47
  48
  49 struct Sum {
  50   double csum; // convolved sum
  51   double sum; // flat sum
  52     int nsamples;
  53   double cmean; //convolved mean
  54     double mean;
  55     double leqm;
  56   double rms;
  57 };
  58
  59 struct WorkerArgs {
  60   double * argbuffer;
  61   int nsamples;
  62   int nch;
  63   int npoints;
  64   double * ir;
  65   struct Sum * ptrtotsum;
  66   double * chconf;
  67   int shorttermindex;
  68   double * shorttermarray;
  69   int leqm10flag;
  70 };
  71
  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);
  89
  90
  91 pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
  92
  93
  94 int main(int argc, const char ** argv)
  95 {
  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
  99   int bitdepth;
 100   // double normalizer;
 101   int timing = 0;
 102   struct timespec starttime;
 103   int fileopenstate = 0;
 104   int leqm10 = 0;
 105   int leqmlog = 0;
 106   int numCPU = sysconf(_SC_NPROCESSORS_ONLN) - 1;
 107   double * channelconfcalvector;
 108   channelconfcalvector = NULL;
 109   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);
 110   //SndfileHandle file;
 111   SNDFILE *file;
 112   file=NULL;
 113   SF_INFO sfinfo;
 114   FILE *leqm10logfile;
 115   leqm10logfile = NULL;
 116   FILE *leqmlogfile;
 117   leqmlogfile = NULL;
 118   int buffersizems = 850; //ISO 21727:2004 do not contain any indication, TASA seems to indicate 1000, p. 8
 119   int buffersizesamples;
 120         double tempchcal[128];
 121         int numcalread = 0;
 122         double * shorttermaveragedarray;
 123         shorttermaveragedarray = NULL;
 124         int numbershortperiods = 0;
 125         int parameterstate = 0;
 126         int leqnw = 0;
 127
 128         char soundfilename[64];
 129         // This is a requirement of sndfile library, do not forget it.
 130
 131         memset(&sfinfo, 0, sizeof(sfinfo));
 132
 133
 134   if (argc == 1)
 135     { 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";
 136       printf(helptext);
 137       printf("Please indicate a sound file to be processed.\n");
 138       return 0;
 139   }
 140
 141
 142     for (int in = 1; in < argc;) {
 143
 144       if (!(strncmp(argv[in], "-", 1) == 0)) {
 145         if (fileopenstate == 0) {
 146           if(! (file = sf_open(argv[in], SFM_READ, &sfinfo))) {
 147             printf("Error while opening audio file, could not open  %s\n.", argv[in]);
 148             puts(sf_strerror(NULL));
 149             return 1;
 150           }
 151
 152           strcpy(soundfilename, argv[in]);
 153              fileopenstate = 1;
 154              printf("Opened file: %s\n", argv[in]);
 155              printf("Sample rate: %d\n", sfinfo.samplerate);
 156              printf("Channels: %d\n", sfinfo.channels);
 157              printf("Format: %d\n", sfinfo.format);
 158              printf("Frames: %d\n", (int) sfinfo.frames);
 159              channelconfcalvector = malloc(sizeof(double) * sfinfo.channels);
 160              in++;
 161              continue;
 162         } else {
 163           free(channelconfcalvector);
 164           channelconfcalvector = NULL;
 165           return 0;
 166         }
 167       }
 168
 169
 170       if (strcmp(argv[in], "-chconfcal") == 0) {
 171         /* as the order of parameter is free I have to postpone
 172            the check for consistency with the number of channels.
 173            So first create a temporary array, whose number of element will be checked after
 174            the parsing of the command line parameters is finished.
 175            The calibration will be expressed in dB on the command line and converted to multiplier
 176            here so that it can be stored as a factor in the channelconfcalvector.
 177         */
 178
 179         in++;
 180         for (;;)  {
 181         if (in < argc) {
 182           //if (!(strncmp(argv[in], "-", 1) == 0)) { //changed this to allow negative numbers
 183             if (!(strncmp(argv[in], "-", 1) == 0) || isdigit(argv[in][1])) {
 184           tempchcal[numcalread++]=atof(argv[in++]);
 185           } else break;
 186         } else break;
 187
 188         } //for
 189         continue;
 190       }
 191
 192          if (strcmp(argv[in], "-convpoints") == 0) {
 193              npoints = atoi(argv[in + 1]);
 194              in+=2;
 195              printf("Convolution points sets to %d.\n", npoints);
 196              continue;
 197
 198       }
 199               if (strcmp(argv[in], "-numcpus") == 0) {
 200                 numCPU= atoi(argv[in + 1]);
 201              in+=2;
 202              printf("Number of threads manually set to %d. Default is number of cores in the system minus one.\n", numCPU);
 203              continue;
 204
 205       }
 206               if (strcmp(argv[in], "-timing") == 0) {
 207                 timing = 1;
 208              in++;
 209              printf("Execution time will be measured.\n");
 210              continue;
 211
 212       }
 213
 214                       if (strcmp(argv[in], "-logleqm10") == 0) {
 215                 leqm10 = 1;
 216              in++;
 217              printf("Leq(M)10 data will be logged to the file leqm10.txt\n");
 218              continue;
 219
 220       }
 221                                       if (strcmp(argv[in], "-logleqm") == 0) {
 222                 leqmlog = 1;
 223              in++;
 224              printf("Leq(M) data will be logged to the file leqmlog.txt\n");
 225              continue;
 226
 227       }
 228
 229                                       if (strcmp(argv[in], "-leqnw") == 0) {
 230                 leqnw = 1;
 231              in++;
 232              printf("Leq(nW) - unweighted -  will be outputted.\n");
 233              continue;
 234
 235       }
 236
 237                                         if (strcmp(argv[in], "-buffersize") == 0) {
 238                 buffersizems = atoi(argv[in + 1]);
 239              in+=2;
 240              printf("Buffersize will be set to %d milliseconds.\n", buffersizems);
 241              continue;
 242
 243       }
 244
 245                                         if (parameterstate==0) {
 246                                           break;
 247                                         }
 248     }
 249 // Open audio file
 250
 251 //postprocessing parameters
 252     if (numcalread == sfinfo.channels) {
 253       for (int cind = 0; cind < sfinfo.channels; cind++) {
 254         channelconfcalvector[cind] = convloglin_single(tempchcal[cind]);
 255
 256       }
 257     }
 258     else if ((numcalread == 0) && (sfinfo.channels == 6)) {
 259         double conf51[6] = {0, 0, 0, 0, -3, -3};
 260         for (int cind = 0; cind < sfinfo.channels; cind++) {
 261           channelconfcalvector[cind] = convloglin_single(conf51[cind]);
 262         }
 263
 264     } else {
 265
 266       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");
 267
 268       free(channelconfcalvector);
 269       channelconfcalvector = NULL;
 270       return 0;
 271     }
 272
 273
 274
 275
 276     if (leqm10) {
 277       char tempstring[128];
 278       strcpy(tempstring, soundfilename);
 279       strcat(tempstring, ".leqm10.txt");
 280       leqm10logfile = fopen(tempstring, "w");
 281       if (leqm10logfile == NULL) {
 282         printf("Could not open file to write log leqm10 data!\n");
 283       }
 284     }
 285
 286
 287
 288
 289     if (leqmlog) {
 290       char tempstring[128];
 291       strcpy(tempstring, soundfilename);
 292       strcat(tempstring, ".leqmlog.txt");
 293       leqmlogfile = fopen(tempstring, "w");
 294       if (leqmlogfile == NULL) {
 295         printf("Could not open file to write log leqm data!\n");
 296       }
 297     }
 298
 299
 300   if (timing) {
 301
 302   clock_gettime(CLOCK_MONOTONIC, &starttime);
 303   }
 304
 305
 306
 307
 308
 309   // reading to a double or float buffer with sndfile take care of normalization
 310  /*
 311  static double  buffer[BUFFER_LEN]; // it seems this must be static. I don't know why
 312  */
 313  double * buffer;
 314  // buffer = new double [BUFFER_LEN];
 315  buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
 316  if ((sfinfo.samplerate*buffersizems)%1000) {
 317    printf("Please fine tune the buffersize according to the sample rate\n");
 318    //close file
 319    // free memory
 320    // write a function to do that
 321    return 1;
 322  }
 323
 324   buffersizesamples = (sfinfo.samplerate*sfinfo.channels*buffersizems)/1000;
 325   buffer = malloc(sizeof(double)*buffersizesamples);
 326
 327  samplingfreq = sfinfo.samplerate;
 328
 329  if(leqm10) {
 330
 331    //if duration < 10 mm exit
 332
 333    double featdursec = sfinfo.frames / sfinfo.samplerate;
 334    if ((featdursec/60.0) < 10.0) {
 335      printf("The audio file is too short to measure Leq(m10).\n");
 336      return 0;
 337    }
 338
 339    //how many short periods in overall duration
 340    int remainder = sfinfo.frames % (sfinfo.samplerate*buffersizems/1000);
 341    if (remainder == 0)  numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000);
 342    else  numbershortperiods = sfinfo.frames/(sfinfo.samplerate*buffersizems/1000) + 1;
 343
 344    //allocate array
 345    shorttermaveragedarray = malloc(sizeof(*shorttermaveragedarray)*numbershortperiods);
 346  }
 347
 348
 349  //End opening audio file
 350
 351   //ISO 21727:2004(E)
 352   // M Weighting
 353   double freqsamples[] = {31, 63, 100, 200, 400, 800, 1000, 2000, 3150, 4000, 5000, 6300, 7100, 8000, 9000, 10000, 12500, 14000, 16000, 20000, 31500};
 354   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};
 355
 356   double * eqfreqresp_db;
 357   eqfreqresp_db = malloc(sizeof(*eqfreqresp_db)*npoints);
 358
 359   double * eqfreqsamples;
 360   eqfreqsamples = malloc(sizeof(*eqfreqsamples)*npoints);
 361   double * eqfreqresp;
 362   eqfreqresp = malloc(sizeof(*eqfreqresp)*npoints);
 363   double * ir;
 364   ir = malloc(sizeof(*ir)*npoints*2);
 365
 366
 367 // And what to do for floating point sample coding?
 368
 369    switch(sfinfo.format & SF_FORMAT_SUBMASK) {
 370    // all signed bitdepth
 371  case 0x0001:
 372    bitdepth = 8;
 373    break;
 374  case 0x0002:
 375    bitdepth = 16;
 376    break;
 377  case 0x0003:
 378    bitdepth = 24;
 379    break;
 380  case 0x0004:
 381    bitdepth = 32;
 382    break;
 383  default:
 384    printf("No known bitdepth! Exiting ...\n");
 385    return -1;
 386    }
 387
 388
 389
 390
 391    equalinterval2(freqsamples, freqresp_db, eqfreqsamples, eqfreqresp_db, npoints, samplingfreq, origpoints, bitdepth);
 392   convloglin(eqfreqresp_db, eqfreqresp, npoints);
 393
 394     #ifdef DEBUG
 395     for (int i=0; i < npoints; i++) {
 396       printf("%d\t%.2f\t%.2f\t%.2f\n", i, eqfreqsamples[i], eqfreqresp_db[i], eqfreqresp[i]);
 397     }
 398     #endif
 399
 400     inversefft2(eqfreqresp, ir, npoints);
 401
 402 // read through the entire file
 403
 404    struct Sum * totsum;
 405     totsum = malloc(sizeof(struct Sum));
 406     totsum->csum = 0.0;
 407     totsum->sum = 0.0;
 408     totsum->nsamples = 0;
 409     totsum->cmean = 0.0;
 410     totsum->mean = 0.0; // Do I write anything here?
 411     totsum->leqm = 0.0;
 412     totsum->rms = 0.0;
 413     sf_count_t samples_read = 0;
 414
 415  // Main loop through audio file
 416
 417  int worker_id = 0;
 418  pthread_t tid[numCPU];
 419  struct WorkerArgs ** WorkerArgsArray;
 420  WorkerArgsArray = malloc(sizeof(struct WorkerArgs *)*numCPU);
 421  int staindex = 0; //shorttermarrayindex
 422
 423
 424  while((samples_read = sf_read_double(file, buffer, buffersizesamples)) > 0) {
 425
 426    WorkerArgsArray[worker_id]=malloc(sizeof(struct WorkerArgs));
 427    WorkerArgsArray[worker_id]->nsamples = samples_read;
 428    WorkerArgsArray[worker_id]->nch = sfinfo.channels;
 429    WorkerArgsArray[worker_id]->npoints=npoints;
 430    WorkerArgsArray[worker_id]->ir = ir;
 431    WorkerArgsArray[worker_id]->ptrtotsum = totsum;
 432
 433    WorkerArgsArray[worker_id]->chconf = channelconfcalvector;
 434    if (leqm10) {
 435    WorkerArgsArray[worker_id]->shorttermindex = staindex++;
 436    WorkerArgsArray[worker_id]->leqm10flag = 1;
 437    WorkerArgsArray[worker_id]->shorttermarray = shorttermaveragedarray;
 438    } else {
 439      WorkerArgsArray[worker_id]->shorttermindex = 0;
 440      WorkerArgsArray[worker_id]->leqm10flag = 0;
 441    }
 442
 443    WorkerArgsArray[worker_id]->argbuffer = malloc(sizeof(double)*buffersizesamples);
 444    memcpy(WorkerArgsArray[worker_id]->argbuffer, buffer, samples_read*sizeof(double));
 445
 446
 447    pthread_attr_t attr;
 448    pthread_attr_init(&attr);
 449    pthread_create(&tid[worker_id], &attr, worker_function, WorkerArgsArray[worker_id]);
 450
 451    worker_id++;
 452    if (worker_id == numCPU) {
 453        worker_id = 0;
 454        //maybe here wait for all cores to output before going on
 455        for (int idxcpu = 0; idxcpu < numCPU; idxcpu++) {
 456        pthread_join(tid[idxcpu], NULL);
 457        free(WorkerArgsArray[idxcpu]);
 458        WorkerArgsArray[idxcpu] = NULL;
 459        }
 460               //simply log here your measurement it will be a multiple of your threads and your buffer
 461        if (leqmlog) {
 462          meanoverduration(totsum); //update leq(m) until now and log it
 463        logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
 464                } //endlog
 465    }
 466
 467
 468
 469     //end while worker_id
 470  /// End looping cores
 471   } // main loop through file
 472
 473  //here I should wait for rest worker (< numcpu)
 474  //but I need to dispose of thread id.
 475  if (worker_id != 0) { // worker_id = 0 means the number of samples was divisible through the number of cpus
 476    for (int idxcpu = 0; idxcpu < worker_id; idxcpu++) { //worker_id is at this point one unit more than threads launched
 477      pthread_join(tid[idxcpu], NULL);
 478      free(WorkerArgsArray[idxcpu]);
 479      WorkerArgsArray[idxcpu] = NULL;
 480    }
 481         //also log here for a last value
 482        if (leqmlog) {
 483          meanoverduration(totsum); //update leq(m) until now and log it
 484        logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
 485                } //endlog
 486         }
 487  // mean of scalar sum over duration
 488
 489  meanoverduration(totsum);
 490  if (leqnw) {
 491  printf("Leq(noW): %.4f\n", totsum->rms); // Leq(no Weighting)
 492  }
 493  printf("Leq(M): %.4f\n", totsum->leqm);
 494
 495   if(timing) {
 496    struct timespec stoptime;
 497    long stoptimenanoseconds;
 498    long executionnanoseconds;
 499    clock_gettime(CLOCK_MONOTONIC, &stoptime);
 500
 501    if (stoptime.tv_nsec < starttime.tv_nsec) {
 502      stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
 503    } else {
 504      stoptimenanoseconds = stoptime.tv_nsec;
 505    }
 506    executionnanoseconds = stoptimenanoseconds - starttime.tv_nsec;
 507    printf("Total execution time is %.6f seconds\n", ((double) stoptime.tv_sec) - ((double) starttime.tv_sec) + ((double) executionnanoseconds / 1000000000.00));
 508  }
 509
 510
 511  if (leqm10) {
 512
 513    //Take the array with the short term accumulators
 514    double interval = 10.0;
 515    //create a rolling average according to rolling interval
 516    int rollint; // in short 10*60 = 600 sec 600/0.850
 517
 518
 519    //how many element of the array to consider for the rollint?
 520    //that is how many buffersizems in the interval - interval could be parameterized(?)
 521    double tempint = 60.0 * interval / (((double) buffersizems) /1000.0);
 522    rollint = (int) tempint;
 523    //dispose of the rest
 524    if (tempint - ((double) rollint) > 0) {
 525      rollint += 1;
 526    }
 527    //two loops
 528    //external loop
 529    int indexlong = 0;
 530    while(indexlong < (numbershortperiods - rollint)) {
 531
 532      double accumulator = 0;
 533      //internal loop
 534      double averagedaccumulator = 0;
 535      for (int indexshort = 0; indexshort < rollint; indexshort++) {
 536
 537        accumulator += shorttermaveragedarray[indexshort+indexlong];
 538      } //end internal loop
 539      averagedaccumulator = accumulator/((double) rollint);
 540      logleqm10(leqm10logfile, ((double) (indexlong+rollint)) * ((double) buffersizems / 1000.0), averagedaccumulator);
 541      indexlong++;
 542    } //end external loop
 543
 544      fclose(leqm10logfile);
 545    free(shorttermaveragedarray);
 546    shorttermaveragedarray = NULL;
 547  }
 548
 549
 550  if (leqmlog) {
 551
 552    fclose(leqmlogfile);
 553  }
 554
 555  sf_close(file);
 556
 557  free(eqfreqsamples);
 558  eqfreqsamples = NULL;
 559   free(eqfreqresp_db);
 560  eqfreqresp_db=NULL;
 561  free(eqfreqresp);
 562  eqfreqresp = NULL;
 563  free(ir);
 564  ir = NULL;
 565  free(channelconfcalvector);
 566  channelconfcalvector = NULL;
 567  free(WorkerArgsArray);
 568  WorkerArgsArray = NULL;
 569
 570  free(totsum);
 571  totsum = NULL;
 572  free(buffer);
 573  buffer=NULL;
 574    return 0;
 575 }
 576
 577
 578
 579 void * worker_function(void * argstruct) {
 580
 581   struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
 582
 583    double * sumandsquarebuffer;
 584    double * csumandsquarebuffer;
 585   double * chsumaccumulator_norm;
 586   double * chsumaccumulator_conv;
 587
 588
 589   sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 590
 591
 592   csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 593
 594   chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 595
 596   chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 597
 598
 599
 600   for (int i = 0; i < thisWorkerArgs->nsamples / thisWorkerArgs->nch; i++) {
 601     sumandsquarebuffer[i] = 0.0;
 602     csumandsquarebuffer[i] = 0.0;
 603   chsumaccumulator_norm[i] = 0.0;
 604   chsumaccumulator_conv[i] = 0.0;
 605   }
 606
 607
 608
 609   for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
 610
 611     double * normalizedbuffer;
 612     double * convolvedbuffer;
 613
 614
 615     normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 616
 617     convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 618
 619
 620     for (int n=ch, m= 0; n < thisWorkerArgs->nsamples; n += thisWorkerArgs->nch, m++) {
 621      // 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
 622       //so not normalized but calibrated
 623    normalizedbuffer[m] = thisWorkerArgs->argbuffer[n]*thisWorkerArgs->chconf[ch]; //this scale amplitude according to specified calibration
 624
 625
 626  }
 627
 628  //convolution
 629  convolv_buff(normalizedbuffer, convolvedbuffer, thisWorkerArgs->ir, thisWorkerArgs->nsamples / thisWorkerArgs->nch, thisWorkerArgs->npoints * 2);
 630  //rectify, square und sum
 631  rectify(csumandsquarebuffer,convolvedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 632  rectify(sumandsquarebuffer,normalizedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 633
 634  accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 635  accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 636
 637
 638  free(normalizedbuffer);
 639  normalizedbuffer= NULL;
 640
 641  free(convolvedbuffer);
 642  convolvedbuffer=NULL;
 643
 644  } // loop through channels
 645
 646     //Create a function for this also a tag so that the worker know if he has to do this or not
 647
 648   if (thisWorkerArgs->leqm10flag) {
 649     thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 650     #ifdef DEBUG
 651     printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
 652     #endif
 653   }
 654   pthread_mutex_lock(&mutex);
 655   // this should be done under mutex conditions -> shared resources!
 656   sumsamples(thisWorkerArgs->ptrtotsum, chsumaccumulator_norm, chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 657   pthread_mutex_unlock(&mutex);
 658
 659
 660   free(sumandsquarebuffer);
 661   sumandsquarebuffer=NULL;
 662
 663   free(csumandsquarebuffer);
 664   csumandsquarebuffer=NULL;
 665
 666   free(chsumaccumulator_norm);
 667   chsumaccumulator_norm=NULL;
 668
 669   free(chsumaccumulator_conv);
 670   chsumaccumulator_conv=NULL;
 671
 672   free(thisWorkerArgs->argbuffer);
 673   thisWorkerArgs->argbuffer = NULL;
 674   // the memory pointed to by this pointer is freed in main
 675   // it is the same memory for all worker
 676   // but it is necessary to set pointer to NULL otherwise free will not work later (really?)
 677   thisWorkerArgs->chconf = NULL;
 678  pthread_exit(0);
 679
 680 }
 681
 682
 683 //to get impulse response frequency response at equally spaced intervals is needed
 684
 685 int equalinterval( double * freqsamples, double  * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
 686     double freq;
 687     // int findex = 0;
 688     // int rindex = 0;
 689     double pass = ((double) (samplingfreq >> 1)) / ((double) points);
 690     for (int ieq = 0, i = 0; ieq < points; ieq++) {
 691         freq = ieq*pass;
 692         eqfreqsamples[ieq] = freq;
 693
 694         if ((freq == 0.0) || (freq < freqsamples[1])) {
 695           eqfreqresp[ieq] = freqresp[0];
 696             continue;
 697     } else {
 698
 699         if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
 700           eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp[i];
 701         } else if (freq >=freqsamples[i+1]) {
 702             while(freq >= freqsamples[i+1]) {
 703                 i++;
 704                 if ((i + 1) >= origpoints) {
 705                   break;
 706                 }
 707             }
 708             if ((i+1) < origpoints) {
 709             eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
 710             } else {
 711               eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
 712             }
 713         }
 714         }
 715     }
 716     return 0;
 717 }
 718
 719
 720
 721
 722
 723 //the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
 724 //it is also different for the way it interpolates between DC and 31 Hz
 725 // Pay attention that also arguments to the functions are changed
 726 int equalinterval2( double freqsamples[], double  freqresp_db[], double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile) {
 727     double freq;
 728
 729
 730     //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
 731     double dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
 732     //double dcatt = -90.3;
 733     double pass = ((double) (samplingfreq >> 1)) / ((double) points);
 734     for (int ieq = 0, i = 0; ieq < points; ieq++) {
 735         freq = ieq*pass;
 736         eqfreqsamples[ieq] = freq;
 737
 738         if (freq == 0.0) {
 739           eqfreqresp[ieq] = dcatt;
 740         } else if (freq < freqsamples[0]) { // this has a lot of influence on final Leq(M) value
 741           eqfreqresp[ieq] = ((freqresp_db[0] - dcatt) / (freqsamples[0] - 0)) * freq + dcatt;
 742           //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
 743             continue;
 744     } else {
 745
 746         if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
 747           eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp_db[i];
 748         } else if (freq >=freqsamples[i+1]) {
 749             while(freq >= freqsamples[i+1]) {
 750                 i++;
 751                 if ((i + 1) >= origpoints) {
 752                   break;
 753                 }
 754             }
 755             if ((i+1) < origpoints) {
 756             eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
 757             } else {
 758               eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
 759             }
 760         }
 761         }
 762     }
 763     return 0;
 764 }
 765
 766
 767
 768
 769
 770
 771 int convloglin(double * in, double * out, int points) {
 772   for (int i = 0; i < points; i++) {
 773     out[i] = powf(10, (in[i]/20.0));
 774   }
 775
 776   return 0;
 777 }
 778
 779
 780 double convlinlog_single(double in) {
 781   double out;
 782     out = log(in)*20.0f;
 783   return out;
 784 }
 785
 786
 787 double convloglin_single(double in) {
 788   double out;
 789   out = powf(10, in/20.0f);
 790   return out;
 791 }
 792
 793 // convolution
 794
 795 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
 796
 797
 798   double  sum = 0.0;
 799   for (int i = 0; i < sigin_dim; i++) {
 800
 801     int m = i;
 802     for (int l = impresp_dim - 1; l >=0; l--,m++) {
 803       if (m >= sigin_dim) {
 804         m -= sigin_dim;
 805       }
 806       sum += sigin[m]*impresp[l];
 807     }
 808     sigout[i] = sum;
 809     sum=0.0;
 810     }
 811   return 0;
 812
 813 }
 814
 815
 816 void  inversefft2(double * eqfreqresp, double * ir, int npoints) {
 817   for (int n = 0; n < npoints; n++) {
 818     double parsum = 0.0;
 819     double partial = 0.0;
 820
 821     for (int m = 1; m <= npoints -1; m++) {
 822       partial = cos(2.0*M_PI*((double) m)*( ( ((double) n) - ( ((double) npoints) * 2.0 -1 ) / 2 ) / ( ((double) npoints) * 2.0) ));
 823       parsum = parsum + eqfreqresp[m]*partial;
 824     }
 825     ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
 826     #ifdef DEBUG
 827     printf("%.4f\n", ir[n]);
 828     #endif
 829   }
 830   for (int n = 0; n < npoints; n++) {
 831     ir[npoints+n] = ir[npoints-(n + 1)];
 832     #ifdef DEBUG
 833     printf("%.4f\n", ir[npoints+n]);
 834     #endif
 835   }
 836
 837
 838 }
 839
 840 // scale input according to required calibration
 841 // this could be different for certain digital cinema formats
 842 double inputcalib (double dbdiffch) {
 843
 844     double coeff = pow(10, dbdiffch/20);
 845     return coeff;
 846
 847 }
 848
 849 //rectify, square and sum
 850 int rectify(double * squared, double * inputsamples, int nsamples){
 851   for (int i = 0; i < nsamples; i++) {
 852     squared[i] = (double) powf(inputsamples[i], 2);
 853     }
 854     return 0;
 855
 856 }
 857
 858 int initbuffer(double * buffertoinit, int nsamples) {
 859   for (int i = 0; i < nsamples; i++) {
 860     buffertoinit[i] = 0.0;
 861
 862   }
 863   return 0;
 864 }
 865
 866 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
 867   for (int i = 0; i < nsamples; i++) {
 868     chaccumulator[i] += inputchannel[i];
 869   }
 870   return 0;
 871 }
 872
 873 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples) {
 874   ts->nsamples += nsamples;
 875   for (int i=0; i < nsamples; i++) {
 876     ts->sum  += inputsamples[i];
 877     ts->csum += cinputsamples[i];
 878   }
 879   return 0;
 880
 881 }
 882
 883 int meanoverduration(struct Sum * oldsum) {
 884   oldsum->mean = pow(oldsum->sum / ((double) oldsum->nsamples), 0.500);
 885    oldsum->cmean = pow(oldsum->csum / ((double) oldsum->nsamples), 0.500);
 886    oldsum->rms = 20*log10(oldsum->mean);
 887    oldsum->leqm = 20*log10(oldsum->cmean) + 110.600;//
 888      // 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.
 889    //this value is obtained calibrating with a -20 dBFS Dolby Tone (RMS) I think this is correct
 890    //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
 891    //This is only approximate as you should use a separate calibration according to the Dolby Format. Also for SW the tone should be 100Hz (?)
 892
 893 return 0;
 894 }
 895
 896 double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
 897   double stsum = 0.0;
 898   for (int i=0; i < nsamples; i++) {
 899     stsum += channelaccumulator[i];
 900
 901   }
 902   return stsum / (double) nsamples;
 903 }
 904
 905 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
 906
 907   fprintf(filehandle, "%.4f", featuretimesec);
 908   fprintf(filehandle, "\t");
 909   fprintf(filehandle, "%.4f\n", oldsum->leqm);
 910
 911
 912 }
 913
 914 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
 915   double leqm10 = 20*log10(pow(longaverage, 0.500)) + 110.600;
 916   fprintf(filehandle, "%.4f", featuretimesec);
 917   fprintf(filehandle, "\t");
 918   fprintf(filehandle, "%.4f\n", leqm10);
 919
 920 }