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]->argbuffer);
 473       WorkerArgsArray[idxcpu]->argbuffer = NULL;
 474        free(WorkerArgsArray[idxcpu]);
 475        WorkerArgsArray[idxcpu] = NULL;
 476        }
 477               //simply log here your measurement it will be a multiple of your threads and your buffer
 478        if (leqmlog) {
 479          meanoverduration(totsum); //update leq(m) until now and log it
 480        logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
 481                } //endlog
 482    }
 483
 484
 485
 486     //end while worker_id
 487  /// End looping cores
 488   } // main loop through file
 489
 490  //here I should wait for rest worker (< numcpu)
 491  //but I need to dispose of thread id.
 492  if (worker_id != 0) { // worker_id = 0 means the number of samples was divisible through the number of cpus
 493    for (int idxcpu = 0; idxcpu < worker_id; idxcpu++) { //worker_id is at this point one unit more than threads launched
 494      pthread_join(tid[idxcpu], NULL);
 495      free(WorkerArgsArray[idxcpu]->argbuffer);
 496      WorkerArgsArray[idxcpu]->argbuffer = NULL;
 497      free(WorkerArgsArray[idxcpu]);
 498      WorkerArgsArray[idxcpu] = NULL;
 499    }
 500         //also log here for a last value
 501        if (leqmlog) {
 502          meanoverduration(totsum); //update leq(m) until now and log it
 503        logleqm(leqmlogfile, ((double) totsum->nsamples)/((double) sfinfo.samplerate), totsum );
 504                } //endlog
 505         }
 506  // mean of scalar sum over duration
 507
 508  meanoverduration(totsum);
 509  if (leqnw) {
 510  printf("Leq(noW): %.4f\n", totsum->rms); // Leq(no Weighting)
 511  }
 512  printf("Leq(M): %.4f\n", totsum->leqm);
 513
 514   if(timing) {
 515    struct timespec stoptime;
 516    long stoptimenanoseconds;
 517    long executionnanoseconds;
 518    clock_gettime(CLOCK_MONOTONIC, &stoptime);
 519
 520    if (stoptime.tv_nsec < starttime.tv_nsec) {
 521      stoptimenanoseconds = 1000000000 + stoptime.tv_nsec;
 522    } else {
 523      stoptimenanoseconds = stoptime.tv_nsec;
 524    }
 525    executionnanoseconds = stoptimenanoseconds - starttime.tv_nsec;
 526    printf("Total execution time is %.6f seconds\n", ((double) stoptime.tv_sec) - ((double) starttime.tv_sec) + ((double) executionnanoseconds / 1000000000.00));
 527  }
 528
 529
 530  if (leqm10) {
 531
 532    //Take the array with the short term accumulators
 533    double interval = 10.0;
 534    //create a rolling average according to rolling interval
 535    int rollint; // in short 10*60 = 600 sec 600/0.850
 536
 537
 538    //how many element of the array to consider for the rollint?
 539    //that is how many buffersizems in the interval - interval could be parameterized(?)
 540    double tempint = 60.0 * interval / (((double) buffersizems) /1000.0);
 541    rollint = (int) tempint;
 542    //dispose of the rest
 543    if (tempint - ((double) rollint) > 0) {
 544      rollint += 1;
 545    }
 546    //two loops
 547    //external loop
 548    int indexlong = 0;
 549    while(indexlong < (numbershortperiods - rollint)) {
 550
 551      double accumulator = 0;
 552      //internal loop
 553      double averagedaccumulator = 0;
 554      for (int indexshort = 0; indexshort < rollint; indexshort++) {
 555
 556        accumulator += shorttermaveragedarray[indexshort+indexlong];
 557      } //end internal loop
 558      averagedaccumulator = accumulator/((double) rollint);
 559      logleqm10(leqm10logfile, ((double) (indexlong+rollint)) * ((double) buffersizems / 1000.0), averagedaccumulator);
 560      indexlong++;
 561    } //end external loop
 562
 563      fclose(leqm10logfile);
 564    free(shorttermaveragedarray);
 565    shorttermaveragedarray = NULL;
 566  }
 567
 568
 569  if (leqmlog) {
 570
 571    fclose(leqmlogfile);
 572  }
 573
 574  sf_close(file);
 575
 576  free(eqfreqsamples);
 577  eqfreqsamples = NULL;
 578   free(eqfreqresp_db);
 579  eqfreqresp_db=NULL;
 580  free(eqfreqresp);
 581  eqfreqresp = NULL;
 582  free(ir);
 583  ir = NULL;
 584  free(channelconfcalvector);
 585  channelconfcalvector = NULL;
 586  free(WorkerArgsArray);
 587  WorkerArgsArray = NULL;
 588
 589  free(totsum);
 590  totsum = NULL;
 591  free(buffer);
 592  buffer=NULL;
 593    return 0;
 594 }
 595
 596
 597
 598 void * worker_function(void * argstruct) {
 599
 600   struct WorkerArgs * thisWorkerArgs = (struct WorkerArgs *) argstruct;
 601
 602    double * sumandsquarebuffer;
 603    double * csumandsquarebuffer;
 604   double * chsumaccumulator_norm;
 605   double * chsumaccumulator_conv;
 606
 607
 608   sumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 609
 610
 611   csumandsquarebuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 612
 613   chsumaccumulator_norm = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 614
 615   chsumaccumulator_conv = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 616
 617
 618
 619   for (int i = 0; i < thisWorkerArgs->nsamples / thisWorkerArgs->nch; i++) {
 620     sumandsquarebuffer[i] = 0.0;
 621     csumandsquarebuffer[i] = 0.0;
 622   chsumaccumulator_norm[i] = 0.0;
 623   chsumaccumulator_conv[i] = 0.0;
 624   }
 625
 626
 627
 628   for (int ch = 0; ch < thisWorkerArgs->nch; ch++) {
 629
 630     double * normalizedbuffer;
 631     double * convolvedbuffer;
 632
 633
 634     normalizedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 635
 636     convolvedbuffer = malloc(sizeof(double)*(thisWorkerArgs->nsamples / thisWorkerArgs->nch));
 637
 638
 639     for (int n=ch, m= 0; n < thisWorkerArgs->nsamples; n += thisWorkerArgs->nch, m++) {
 640      // 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
 641       //so not normalized but calibrated
 642    normalizedbuffer[m] = thisWorkerArgs->argbuffer[n]*thisWorkerArgs->chconf[ch]; //this scale amplitude according to specified calibration
 643
 644
 645  }
 646
 647  //convolution
 648  convolv_buff(normalizedbuffer, convolvedbuffer, thisWorkerArgs->ir, thisWorkerArgs->nsamples / thisWorkerArgs->nch, thisWorkerArgs->npoints * 2);
 649  //rectify, square und sum
 650  rectify(csumandsquarebuffer,convolvedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 651  rectify(sumandsquarebuffer,normalizedbuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 652
 653  accumulatech(chsumaccumulator_norm, sumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 654  accumulatech(chsumaccumulator_conv, csumandsquarebuffer, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 655
 656
 657  free(normalizedbuffer);
 658  normalizedbuffer= NULL;
 659
 660  free(convolvedbuffer);
 661  convolvedbuffer=NULL;
 662
 663  } // loop through channels
 664
 665     //Create a function for this also a tag so that the worker know if he has to do this or not
 666
 667   if (thisWorkerArgs->leqm10flag) {
 668     thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex] = sumandshorttermavrg(chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 669     #ifdef DEBUG
 670     printf("%d: %.6f\n", thisWorkerArgs->shorttermindex, thisWorkerArgs->shorttermarray[thisWorkerArgs->shorttermindex]);
 671     #endif
 672   }
 673   pthread_mutex_lock(&mutex);
 674   // this should be done under mutex conditions -> shared resources!
 675   sumsamples(thisWorkerArgs->ptrtotsum, chsumaccumulator_norm, chsumaccumulator_conv, thisWorkerArgs->nsamples / thisWorkerArgs->nch);
 676   pthread_mutex_unlock(&mutex);
 677
 678
 679   free(sumandsquarebuffer);
 680   sumandsquarebuffer=NULL;
 681
 682   free(csumandsquarebuffer);
 683   csumandsquarebuffer=NULL;
 684
 685   free(chsumaccumulator_norm);
 686   chsumaccumulator_norm=NULL;
 687
 688   free(chsumaccumulator_conv);
 689   chsumaccumulator_conv=NULL;
 690
 691   free(thisWorkerArgs->argbuffer);
 692   thisWorkerArgs->argbuffer = NULL;
 693   // the memory pointed to by this pointer is freed in main
 694   // it is the same memory for all worker
 695   // but it is necessary to set pointer to NULL otherwise free will not work later (really?)
 696   thisWorkerArgs->chconf = NULL;
 697  pthread_exit(0);
 698
 699 }
 700
 701
 702 //to get impulse response frequency response at equally spaced intervals is needed
 703
 704 int equalinterval( double * freqsamples, double  * freqresp, double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints) {
 705     double freq;
 706     // int findex = 0;
 707     // int rindex = 0;
 708     double pass = ((double) (samplingfreq >> 1)) / ((double) points);
 709     for (int ieq = 0, i = 0; ieq < points; ieq++) {
 710         freq = ieq*pass;
 711         eqfreqsamples[ieq] = freq;
 712
 713         if ((freq == 0.0) || (freq < freqsamples[1])) {
 714           eqfreqresp[ieq] = freqresp[0];
 715             continue;
 716     } else {
 717
 718         if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
 719           eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp[i];
 720         } else if (freq >=freqsamples[i+1]) {
 721             while(freq >= freqsamples[i+1]) {
 722                 i++;
 723                 if ((i + 1) >= origpoints) {
 724                   break;
 725                 }
 726             }
 727             if ((i+1) < origpoints) {
 728             eqfreqresp[ieq] = ((freqresp[i+1] - freqresp[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
 729             } else {
 730               eqfreqresp[ieq] = ((1 - freqresp[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp[i];
 731             }
 732         }
 733         }
 734     }
 735     return 0;
 736 }
 737
 738
 739
 740
 741
 742 //the following is different from version 1 because interpolate between db and not linear. Conversion from db to lin must be done after.
 743 //it is also different for the way it interpolates between DC and 31 Hz
 744 // Pay attention that also arguments to the functions are changed
 745 int equalinterval2( double freqsamples[], double  freqresp_db[], double * eqfreqsamples, double * eqfreqresp, int points, int samplingfreq, int origpoints, int bitdepthsoundfile) {
 746     double freq;
 747
 748
 749     //calculate miminum attenuation depending on the bitdeph (minus one), that is −6.020599913 dB per bit in eccess to sign
 750     double dcatt = ((double) (bitdepthsoundfile - 1))*(-6.020599913) + 20.00; //in dB
 751     //double dcatt = -90.3;
 752     double pass = ((double) (samplingfreq >> 1)) / ((double) points);
 753     for (int ieq = 0, i = 0; ieq < points; ieq++) {
 754         freq = ieq*pass;
 755         eqfreqsamples[ieq] = freq;
 756
 757         if (freq == 0.0) {
 758           eqfreqresp[ieq] = dcatt;
 759         } else if (freq < freqsamples[0]) { // this has a lot of influence on final Leq(M) value
 760           eqfreqresp[ieq] = ((freqresp_db[0] - dcatt) / (freqsamples[0] - 0)) * freq + dcatt;
 761           //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
 762             continue;
 763     } else {
 764
 765         if ((freq >= freqsamples[i]) && (freq < freqsamples[i+1])) {
 766           eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq - freqsamples[i]) + freqresp_db[i];
 767         } else if (freq >=freqsamples[i+1]) {
 768             while(freq >= freqsamples[i+1]) {
 769                 i++;
 770                 if ((i + 1) >= origpoints) {
 771                   break;
 772                 }
 773             }
 774             if ((i+1) < origpoints) {
 775             eqfreqresp[ieq] = ((freqresp_db[i+1] - freqresp_db[i])/(freqsamples[i+1] - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
 776             } else {
 777               eqfreqresp[ieq] = ((1 - freqresp_db[i])/(((double) (samplingfreq >> 1)) - freqsamples[i]))*(freq- freqsamples[i]) + freqresp_db[i];
 778             }
 779         }
 780         }
 781     }
 782     return 0;
 783 }
 784
 785
 786
 787
 788
 789
 790 int convloglin(double * in, double * out, int points) {
 791   for (int i = 0; i < points; i++) {
 792     out[i] = powf(10, (in[i]/20.0));
 793   }
 794
 795   return 0;
 796 }
 797
 798
 799 double convlinlog_single(double in) {
 800   double out;
 801     out = log(in)*20.0f;
 802   return out;
 803 }
 804
 805
 806 double convloglin_single(double in) {
 807   double out;
 808   out = powf(10, in/20.0f);
 809   return out;
 810 }
 811
 812 // convolution
 813
 814 int convolv_buff(double * sigin, double * sigout, double * impresp, int sigin_dim, int impresp_dim) {
 815
 816
 817   double  sum = 0.0;
 818   for (int i = 0; i < sigin_dim; i++) {
 819
 820     int m = i;
 821     for (int l = impresp_dim - 1; l >=0; l--,m++) {
 822       if (m >= sigin_dim) {
 823         m -= sigin_dim;
 824       }
 825       sum += sigin[m]*impresp[l];
 826     }
 827     sigout[i] = sum;
 828     sum=0.0;
 829     }
 830   return 0;
 831
 832 }
 833
 834
 835 void  inversefft2(double * eqfreqresp, double * ir, int npoints) {
 836   for (int n = 0; n < npoints; n++) {
 837     double parsum = 0.0;
 838     double partial = 0.0;
 839
 840     for (int m = 1; m <= npoints -1; m++) {
 841       partial = cos(2.0*M_PI*((double) m)*( ( ((double) n) - ( ((double) npoints) * 2.0 -1 ) / 2 ) / ( ((double) npoints) * 2.0) ));
 842       parsum = parsum + eqfreqresp[m]*partial;
 843     }
 844     ir[n] = (eqfreqresp[0] + 2.0 * parsum)/((double) npoints * 2.0);
 845     #ifdef DEBUG
 846     printf("%.4f\n", ir[n]);
 847     #endif
 848   }
 849   for (int n = 0; n < npoints; n++) {
 850     ir[npoints+n] = ir[npoints-(n + 1)];
 851     #ifdef DEBUG
 852     printf("%.4f\n", ir[npoints+n]);
 853     #endif
 854   }
 855
 856
 857 }
 858
 859 // scale input according to required calibration
 860 // this could be different for certain digital cinema formats
 861 double inputcalib (double dbdiffch) {
 862
 863     double coeff = pow(10, dbdiffch/20);
 864     return coeff;
 865
 866 }
 867
 868 //rectify, square and sum
 869 int rectify(double * squared, double * inputsamples, int nsamples){
 870   for (int i = 0; i < nsamples; i++) {
 871     squared[i] = (double) powf(inputsamples[i], 2);
 872     }
 873     return 0;
 874
 875 }
 876
 877 int initbuffer(double * buffertoinit, int nsamples) {
 878   for (int i = 0; i < nsamples; i++) {
 879     buffertoinit[i] = 0.0;
 880
 881   }
 882   return 0;
 883 }
 884
 885 int accumulatech(double * chaccumulator, double * inputchannel, int nsamples) {
 886   for (int i = 0; i < nsamples; i++) {
 887     chaccumulator[i] += inputchannel[i];
 888   }
 889   return 0;
 890 }
 891
 892 int sumsamples(struct Sum * ts, double * inputsamples, double * cinputsamples, int nsamples) {
 893   ts->nsamples += nsamples;
 894   for (int i=0; i < nsamples; i++) {
 895     ts->sum  += inputsamples[i];
 896     ts->csum += cinputsamples[i];
 897   }
 898   return 0;
 899
 900 }
 901
 902 int meanoverduration(struct Sum * oldsum) {
 903   oldsum->mean = pow(oldsum->sum / ((double) oldsum->nsamples), 0.500);
 904    oldsum->cmean = pow(oldsum->csum / ((double) oldsum->nsamples), 0.500);
 905    oldsum->rms = 20*log10(oldsum->mean);
 906    oldsum->leqm = 20*log10(oldsum->cmean) +  108.0851;//
 907      //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.
 908    //this value is obtained calibrating with a -20 dBFS.
 909
 910 return 0;
 911 }
 912
 913 double sumandshorttermavrg(double * channelaccumulator, int nsamples) {
 914   double stsum = 0.0;
 915   for (int i=0; i < nsamples; i++) {
 916     stsum += channelaccumulator[i];
 917
 918   }
 919   return stsum / (double) nsamples;
 920 }
 921
 922 void logleqm(FILE * filehandle, double featuretimesec, struct Sum * oldsum) {
 923
 924   fprintf(filehandle, "%.4f", featuretimesec);
 925   fprintf(filehandle, "\t");
 926   fprintf(filehandle, "%.4f\n", oldsum->leqm);
 927
 928
 929 }
 930
 931 void logleqm10(FILE * filehandle, double featuretimesec, double longaverage) {
 932   double leqm10 = 20*log10(pow(longaverage, 0.500)) +  108.0851;
 933   fprintf(filehandle, "%.4f", featuretimesec);
 934   fprintf(filehandle, "\t");
 935   fprintf(filehandle, "%.4f\n", leqm10);
 936
 937 }