leqm-nrt.c

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