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