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