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