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