cinelerra/resample.C

   1 #include "clip.h"
   2 #include "file.h"
   3 #include "resample.h"
   4
   5 #include <math.h>
   6 #include <stdio.h>
   7 #include <string.h>
   8
   9 // Resampling from Lame
  10
  11 Resample::Resample(File *file, int channels)
  12 {
  13 //printf("Resample::Resample 1 %d\n", channels);
  14         this->file = file;
  15         this->channels = channels;
  16
  17         old = new double*[channels];
  18         for(int i = 0; i < channels; i++)
  19         {
  20                 old[i] = new double[BLACKSIZE];
  21         }
  22         itime = new double[channels];
  23         output_temp_start = new long[channels];
  24         bzero(output_temp_start, sizeof(long) * channels);
  25         resample_init = new int[channels];
  26         bzero(resample_init, sizeof(int) * channels);
  27         last_ratio = 0;
  28         output_temp = 0;
  29         output_size = new long[channels];
  30         bzero(output_size, sizeof(long) * channels);
  31         output_allocation = 0;
  32         input_size = RESAMPLE_CHUNKSIZE;
  33         input_chunk_end = new long[channels];
  34         bzero(input_chunk_end, sizeof(long) * channels);
  35         input = new double[input_size];
  36         last_out_end = new long[channels];
  37         bzero(last_out_end, sizeof(long) * channels);
  38 //printf("Resample::Resample 2 %d\n", channels);
  39 }
  40
  41
  42 Resample::~Resample()
  43 {
  44         for(int i = 0; i < channels; i++)
  45         {
  46                 delete [] old[i];
  47         }
  48         if(output_temp)
  49         {
  50                 for(int i = 0; i < channels; i++)
  51                 {
  52                         delete [] output_temp[i];
  53                 }
  54                 delete [] output_temp;
  55         }
  56
  57         delete [] input_chunk_end;
  58         delete [] input;
  59         delete [] old;
  60         delete [] itime;
  61         delete [] output_temp_start;
  62         delete [] output_size;
  63         delete [] last_out_end;
  64 }
  65
  66 void Resample::reset(int channel)
  67 {
  68 //printf("Resample::reset 1 channel=%d normalized_sample_rate=%d\n", channel, file->normalized_sample_rate);
  69         if(channel < 0)
  70         {
  71                 bzero(resample_init, sizeof(int) * channels);
  72                 bzero(output_size, sizeof(long) * channels);
  73                 bzero(last_out_end, sizeof(long) * channels);
  74                 bzero(input_chunk_end, sizeof(long) * channels);
  75         }
  76         else
  77         {
  78                 resample_init[channel] = 0;
  79                 output_size[channel] = 0;
  80                 last_out_end[channel] = 0;
  81                 input_chunk_end[channel] = 0;
  82         }
  83 }
  84
  85 double Resample::blackman(int i, double offset, double fcn, int l)
  86 {
  87   /* This algorithm from:
  88 SIGNAL PROCESSING ALGORITHMS IN FORTRAN AND C
  89 S.D. Stearns and R.A. David, Prentice-Hall, 1992
  90   */
  91
  92         double bkwn;
  93         double wcn = (M_PI * fcn);
  94         double dly = l / 2.0;
  95         double x = i-offset;
  96         if(x < 0) x = 0;
  97         else
  98         if(x > l) x = l;
  99
 100         bkwn = 0.42 - 0.5 * cos((x * 2) * M_PI /l) + 0.08 * cos((x * 4) * M_PI /l);
 101         if(fabs(x - dly) < 1e-9)
 102                 return wcn / M_PI;
 103     else
 104         return (sin((wcn * (x - dly))) / (M_PI * (x - dly)) * bkwn);
 105 }
 106
 107
 108 int Resample::get_output_size(int channel)
 109 {
 110         return output_size[channel];
 111 }
 112
 113 void Resample::read_output(double *output, int channel, int size)
 114 {
 115         memcpy(output, output_temp[channel], size * sizeof(double));
 116 // Shift leftover forward
 117         for(int i = size; i < output_size[channel]; i++)
 118                 output_temp[channel][i - size] = output_temp[channel][i];
 119         output_size[channel] -= size;
 120 }
 121
 122
 123
 124 void Resample::resample_chunk(double *input,
 125         long in_len,
 126         int in_rate,
 127         int out_rate,
 128         int channel)
 129 {
 130         double resample_ratio = (double)in_rate / out_rate;
 131         int filter_l;
 132         double fcn, intratio;
 133         double offset, xvalue;
 134         int num_used;
 135         int i, j, k;
 136
 137         intratio = (fabs(resample_ratio - floor(.5 + resample_ratio)) < .0001);
 138         fcn = .90 / resample_ratio;
 139         if(fcn > .90) fcn = .90;
 140         filter_l = BLACKSIZE - 6;
 141 /* must be odd */
 142         if(0 == filter_l % 2 ) --filter_l;
 143
 144 /* if resample_ratio = int, filter_l should be even */
 145         filter_l += (int)intratio;
 146
 147 // Blackman filter initialization must be called whenever there is a
 148 // sampling ratio change
 149         if(!resample_init[channel] || last_ratio != resample_ratio)
 150         {
 151                 resample_init[channel] = 1;
 152                 itime[channel] = 0;
 153                 bzero(old[channel], sizeof(double) * BLACKSIZE);
 154
 155 // precompute blackman filter coefficients
 156         for (j = 0; j <= 2 * BPC; ++j)
 157                 {
 158                         for(j = 0; j <= 2 * BPC; j++)
 159                         {
 160                                 offset = (double)(j - BPC) / (2 * BPC);
 161                                 for(i = 0; i <= filter_l; i++)
 162                                 {
 163                                         blackfilt[j][i] = blackman(i, offset, fcn, filter_l);
 164                                 }
 165                         }
 166                 }
 167         }
 168
 169 // Main loop
 170         double *inbuf_old = old[channel];
 171         for(k = 0; 1; k++)
 172         {
 173                 double time0;
 174                 int joff;
 175
 176                 time0 = k * resample_ratio;
 177                 j = (int)floor(time0 - itime[channel]);
 178
 179 //              if(j + filter_l / 2 >= input_size) break;
 180                 if(j + (filter_l + 1) / 2 >= in_len) break;
 181
 182 /* blackman filter.  by default, window centered at j+.5(filter_l%2) */
 183 /* but we want a window centered at time0.   */
 184                 offset = (time0 - itime[channel] - (j + .5 * (filter_l % 2)));
 185                 joff = (int)floor((offset * 2 * BPC) + BPC + .5);
 186                 xvalue = 0;
 187
 188                 for(i = 0; i <= filter_l; i++)
 189                 {
 190                         int j2 = i + j - filter_l / 2;
 191                         double y = ((j2 < 0) ? inbuf_old[BLACKSIZE + j2] : input[j2]);
 192
 193                         xvalue += y * blackfilt[joff][i];
 194                 }
 195
 196                 if(output_allocation <= output_size[channel])
 197                 {
 198                         double **new_output = new double*[channels];
 199                         long new_allocation = output_allocation ? (output_allocation * 2) : 16384;
 200                         for(int l = 0; l < channels; l++)
 201                         {
 202                                 new_output[l] = new double[new_allocation];
 203                                 if(output_temp)
 204                                 {
 205                                         bcopy(output_temp[l], new_output[l], output_allocation * sizeof(double));
 206                                         delete [] output_temp[l];
 207                                 }
 208                         }
 209
 210                         if(output_temp) delete [] output_temp;
 211                         output_temp = new_output;
 212                         output_allocation = new_allocation;
 213                 }
 214
 215                 output_temp[channel][output_size[channel]++] = xvalue;
 216         }
 217
 218         num_used = MIN(in_len, j + filter_l / 2);
 219         itime[channel] += num_used - k * resample_ratio;
 220 //      for(i = 0; i < BLACKSIZE; i++)
 221 //              inbuf_old[i] = input[num_used + i - BLACKSIZE];
 222         bcopy(input + num_used - BLACKSIZE, inbuf_old, BLACKSIZE * sizeof(double));
 223
 224         last_ratio = resample_ratio;
 225 }
 226
 227 void Resample::read_chunk(double *input, long len, int &reseek, int iteration)
 228 {
 229 //printf("Resample::read_chunk 1\n");
 230         if(reseek)
 231         {
 232                 file->set_audio_position(file->current_sample, 0);
 233                 reseek= 0;
 234         }
 235         else
 236         if(iteration == 0)
 237         {
 238 // Resume at the end of the last resample call
 239                 file->set_audio_position(input_chunk_end[file->current_channel], 0);
 240         }
 241
 242         file->read_samples(input, len, 0);
 243         input_chunk_end[file->current_channel] = file->current_sample;
 244
 245 //printf("Resample::read_chunk 2\n");
 246 }
 247
 248 int Resample::resample(double *output,
 249         long out_len,
 250         int in_rate,
 251         int out_rate,
 252         int channel,
 253         long in_position,
 254         long out_position)
 255 {
 256         int total_input = 0;
 257         int reseek = 0;
 258
 259 #define REPOSITION(x, y) \
 260         (labs((x) - (y)) > 1)
 261
 262
 263
 264 //printf("Resample::resample 1 last_out_end=%d out_position=%d\n", last_out_end[channel], out_position);
 265
 266         if(REPOSITION(last_out_end[channel], out_position))
 267         {
 268                 reseek = 1;
 269                 reset(channel);
 270         }
 271
 272
 273
 274
 275
 276
 277         output_temp_start[channel] = file->get_audio_position(out_rate) + out_len;
 278         last_out_end[channel] = out_position + out_len;
 279
 280         int i = 0;
 281         while(out_len > 0)
 282         {
 283 // Drain output buffer
 284                 if(output_size[channel])
 285                 {
 286                         int fragment_len = output_size[channel];
 287                         if(fragment_len > out_len) fragment_len = out_len;
 288
 289 //printf("Resample::resample 1 %d %d %d\n", out_len, output_size[channel], channel);
 290                         bcopy(output_temp[channel], output, fragment_len * sizeof(double));
 291
 292 // Shift leftover forward
 293 //                      for(int i = fragment_len; i < output_size[channel]; i++)
 294 //                              output_temp[channel][i - fragment_len] = output_temp[channel][i];
 295                         bcopy(output_temp[channel] + fragment_len, output_temp[channel], (output_size[channel] - fragment_len) * sizeof(double));
 296
 297                         output_size[channel] -= fragment_len;
 298                         out_len -= fragment_len;
 299                         output += fragment_len;
 300                 }
 301
 302 // Import new samples
 303 //printf("Resample::resample 2 %d %d\n", out_len, channel);
 304                 if(out_len > 0)
 305                 {
 306 //printf("Resample::resample 3 input_size=%d reseek=%d out_position=%d channel=%d\n", input_size, reseek, out_position, channel);
 307                         read_chunk(input, input_size, reseek, i);
 308                         resample_chunk(input,
 309                                 input_size,
 310                                 in_rate,
 311                                 out_rate,
 312                                 channel);
 313                         total_input += input_size;
 314                 }
 315
 316                 i++;
 317         }
 318 //printf("Resample::resample 2 %d %d\n", last_out_end[channel], out_position);
 319 //printf("Resample::resample 2 %d %d %d\n", out_len, output_size[channel], channel);
 320
 321 //printf("Resample::resample 2 %d %d\n", channel, output_size[channel]);
 322
 323         return total_input;
 324 }
 325
 326
 327 Resample_float::Resample_float(File *file, int channels)
 328 {
 329 //printf("Resample_float::Resample_float 1 %d\n", channels);
 330         this->file = file;
 331         this->channels = channels;
 332
 333         old = new float*[channels];
 334         for(int i = 0; i < channels; i++)
 335         {
 336                 old[i] = new float[BLACKSIZE];
 337         }
 338         itime = new float[channels];
 339         output_temp_start = new long[channels];
 340         bzero(output_temp_start, sizeof(long) * channels);
 341         resample_init = new int[channels];
 342         bzero(resample_init, sizeof(int) * channels);
 343         last_ratio = 0;
 344         output_temp = 0;
 345         output_size = new long[channels];
 346         bzero(output_size, sizeof(long) * channels);
 347         output_allocation = 0;
 348         input_size = RESAMPLE_CHUNKSIZE;
 349         input_chunk_end = new long[channels];
 350         bzero(input_chunk_end, sizeof(long) * channels);
 351         input = new float[input_size];
 352         last_out_end = new long[channels];
 353         bzero(last_out_end, sizeof(long) * channels);
 354 }
 355
 356
 357 Resample_float::~Resample_float()
 358 {
 359         for(int i = 0; i < channels; i++)
 360         {
 361                 delete [] old[i];
 362         }
 363         if(output_temp)
 364         {
 365                 for(int i = 0; i < channels; i++)
 366                 {
 367                         delete [] output_temp[i];
 368                 }
 369                 delete [] output_temp;
 370         }
 371
 372         delete [] input_chunk_end;
 373         delete [] input;
 374         delete [] old;
 375         delete [] itime;
 376         delete [] output_temp_start;
 377         delete [] output_size;
 378         delete [] last_out_end;
 379 }
 380
 381 void Resample_float::reset(int channel)
 382 {
 383 //printf("Resample_float::reset 1 channel=%d normalized_sample_rate=%d\n", channel, file->normalized_sample_rate);
 384         if(channel < 0)
 385         {
 386                 bzero(resample_init, sizeof(int) * channels);
 387                 bzero(output_size, sizeof(long) * channels);
 388                 bzero(last_out_end, sizeof(long) * channels);
 389                 bzero(input_chunk_end, sizeof(long) * channels);
 390         }
 391         else
 392         {
 393                 resample_init[channel] = 0;
 394                 output_size[channel] = 0;
 395                 last_out_end[channel] = 0;
 396                 input_chunk_end[channel] = 0;
 397         }
 398 }
 399
 400 float Resample_float::blackman(int i, float offset, float fcn, int l)
 401 {
 402   /* This algorithm from:
 403 SIGNAL PROCESSING ALGORITHMS IN FORTRAN AND C
 404 S.D. Stearns and R.A. David, Prentice-Hall, 1992
 405   */
 406
 407         float bkwn;
 408         float wcn = (M_PI * fcn);
 409         float dly = l / 2.0;
 410         float x = i-offset;
 411         if(x < 0) x = 0;
 412         else
 413         if(x > l) x = l;
 414
 415         bkwn = 0.42 - 0.5 * cos((x * 2) * M_PI /l) + 0.08 * cos((x * 4) * M_PI /l);
 416         if(fabs(x - dly) < 1e-9)
 417                 return wcn / M_PI;
 418     else
 419         return (sin((wcn * (x - dly))) / (M_PI * (x - dly)) * bkwn);
 420 }
 421
 422
 423 int Resample_float::get_output_size(int channel)
 424 {
 425         return output_size[channel];
 426 }
 427
 428 void Resample_float::read_output(double *output, int channel, int size)
 429 {
 430         memcpy(output, output_temp[channel], size * sizeof(double));
 431 // Shift leftover forward
 432         for(int i = size; i < output_size[channel]; i++)
 433                 output_temp[channel][i - size] = output_temp[channel][i];
 434         output_size[channel] -= size;
 435 }
 436
 437
 438
 439 void Resample_float::resample_chunk(float *input,
 440         long in_len,
 441         int in_rate,
 442         int out_rate,
 443         int channel)
 444 {
 445         float resample_ratio = (float)in_rate / out_rate;
 446         int filter_l;
 447         float fcn, intratio;
 448         float offset, xvalue;
 449         int num_used;
 450         int i, j, k;
 451
 452 //printf("Resample_float::resample_chunk 1\n");
 453         intratio = (fabs(resample_ratio - floor(.5 + resample_ratio)) < .0001);
 454         fcn = .90 / resample_ratio;
 455         if(fcn > .90) fcn = .90;
 456         filter_l = BLACKSIZE - 6;
 457 /* must be odd */
 458         if(0 == filter_l % 2 ) --filter_l;
 459
 460 //printf("Resample_float::resample_chunk 2\n");
 461 /* if resample_ratio = int, filter_l should be even */
 462         filter_l += (int)intratio;
 463
 464 //printf("Resample_float::resample_chunk 3\n");
 465 // Blackman filter initialization must be called whenever there is a
 466 // sampling ratio change
 467         if(!resample_init[channel] || last_ratio != resample_ratio)
 468         {
 469                 resample_init[channel] = 1;
 470                 itime[channel] = 0;
 471                 bzero(old[channel], sizeof(float) * BLACKSIZE);
 472
 473 //printf("Resample_float::resample_chunk 4\n");
 474 // precompute blackman filter coefficients
 475         for (j = 0; j <= 2 * BPC; ++j)
 476                 {
 477                         for(j = 0; j <= 2 * BPC; j++)
 478                         {
 479                                 offset = (float)(j - BPC) / (2 * BPC);
 480                                 for(i = 0; i <= filter_l; i++)
 481                                 {
 482                                         blackfilt[j][i] = blackman(i, offset, fcn, filter_l);
 483                                 }
 484                         }
 485                 }
 486         }
 487
 488 //printf("Resample_float::resample_chunk 5\n");
 489 // Main loop
 490         float *inbuf_old = old[channel];
 491         for(k = 0; 1; k++)
 492         {
 493                 float time0;
 494                 int joff;
 495
 496 //printf("Resample_float::resample_chunk 6\n");
 497                 time0 = k * resample_ratio;
 498                 j = (int)floor(time0 - itime[channel]);
 499
 500 //              if(j + filter_l / 2 >= input_size) break;
 501                 if(j + (filter_l + 1) / 2 >= in_len) break;
 502
 503 //printf("Resample_float::resample_chunk 7\n");
 504 /* blackman filter.  by default, window centered at j+.5(filter_l%2) */
 505 /* but we want a window centered at time0.   */
 506                 offset = (time0 - itime[channel] - (j + .5 * (filter_l % 2)));
 507                 joff = (int)floor((offset * 2 * BPC) + BPC + .5);
 508                 xvalue = 0;
 509
 510 //printf("Resample_float::resample_chunk 8\n");
 511                 for(i = 0; i <= filter_l; i++)
 512                 {
 513                         int j2 = i + j - filter_l / 2;
 514                         float y = ((j2 < 0) ? inbuf_old[BLACKSIZE + j2] : input[j2]);
 515
 516 //printf("Resample_float::resample_chunk 9\n");
 517                         xvalue += (double)y * blackfilt[joff][i];
 518                 }
 519
 520 //printf("Resample_float::resample_chunk 10\n");
 521                 if(output_allocation <= output_size[channel])
 522                 {
 523                         double **new_output = new double*[channels];
 524                         long new_allocation = output_allocation ? (output_allocation * 2) : 16384;
 525                         for(int l = 0; l < channels; l++)
 526                         {
 527                                 new_output[l] = new double[new_allocation];
 528                                 if(output_temp)
 529                                 {
 530                                         bcopy(output_temp[l], new_output[l], output_allocation * sizeof(double));
 531                                         delete [] output_temp[l];
 532                                 }
 533                         }
 534
 535                         if(output_temp) delete [] output_temp;
 536                         output_temp = new_output;
 537                         output_allocation = new_allocation;
 538                 }
 539
 540 //printf("Resample_float::resample_chunk 11 %d %d\n", output_size[channel], output_allocation);
 541                 output_temp[channel][output_size[channel]++] = xvalue;
 542         }
 543
 544 //printf("Resample_float::resample_chunk 12\n");
 545         num_used = MIN(in_len, j + filter_l / 2);
 546         itime[channel] += num_used - k * resample_ratio;
 547 //      for(i = 0; i < BLACKSIZE; i++)
 548 //              inbuf_old[i] = input[num_used + i - BLACKSIZE];
 549         bcopy(input + num_used - BLACKSIZE, inbuf_old, BLACKSIZE * sizeof(float));
 550
 551 //printf("Resample_float::resample_chunk 13\n");
 552         last_ratio = resample_ratio;
 553 }
 554
 555 void Resample_float::read_chunk(float *input, long len, int &reseek, int iteration)
 556 {
 557 //printf("Resample_float::read_chunk 1\n");
 558         if(reseek)
 559         {
 560                 file->set_audio_position(file->current_sample, 0);
 561                 reseek= 0;
 562         }
 563         else
 564         if(iteration == 0)
 565         {
 566 // Resume at the end of the last resample call
 567                 file->set_audio_position(input_chunk_end[file->current_channel], 0);
 568         }
 569
 570         file->read_samples(0, len, 0, input);
 571         input_chunk_end[file->current_channel] = file->current_sample;
 572
 573 //printf("Resample_float::read_chunk 2\n");
 574 }
 575
 576 int Resample_float::resample(double *output,
 577         long out_len,
 578         int in_rate,
 579         int out_rate,
 580         int channel,
 581         long in_position,
 582         long out_position)
 583 {
 584         int total_input = 0;
 585         int reseek = 0;
 586
 587 #define REPOSITION(x, y) \
 588         (labs((x) - (y)) > 1)
 589
 590
 591
 592 //printf("Resample_float::resample 1 last_out_end=%d out_position=%d\n", last_out_end[channel], out_position);
 593
 594         if(REPOSITION(last_out_end[channel], out_position))
 595         {
 596                 reseek = 1;
 597                 reset(channel);
 598         }
 599
 600         output_temp_start[channel] = file->get_audio_position(out_rate) + out_len;
 601         last_out_end[channel] = out_position + out_len;
 602
 603         int i = 0;
 604         while(out_len > 0)
 605         {
 606 // Drain output buffer
 607                 if(output_size[channel])
 608                 {
 609                         int fragment_len = output_size[channel];
 610                         if(fragment_len > out_len) fragment_len = out_len;
 611
 612 //printf("Resample_float::resample 1 %d %d %d\n", out_len, output_size[channel], channel);
 613                         bcopy(output_temp[channel], output, fragment_len * sizeof(double));
 614
 615
 616 // Shift leftover forward
 617                         //for(int i = fragment_len; i < output_size[channel]; i++)
 618                         //      output_temp[channel][i - fragment_len] = output_temp[channel][i];
 619                         bcopy(output_temp[channel] + fragment_len, output_temp[channel], (output_size[channel] - fragment_len) * sizeof(double));
 620
 621                         output_size[channel] -= fragment_len;
 622                         out_len -= fragment_len;
 623                         output += fragment_len;
 624                 }
 625
 626 // Import new samples
 627 //printf("Resample_float::resample 2 %d %d\n", out_len, channel);
 628                 if(out_len > 0)
 629                 {
 630 //printf("Resample_float::resample 3 input_size=%d reseek=%d out_position=%d channel=%d\n", input_size, reseek, out_position, channel);
 631                         read_chunk(input, input_size, reseek, i);
 632                         resample_chunk(input,
 633                                 input_size,
 634                                 in_rate,
 635                                 out_rate,
 636                                 channel);
 637                         total_input += input_size;
 638                 }
 639
 640                 i++;
 641         }
 642 //printf("Resample_float::resample 2 %d %d\n", last_out_end[channel], out_position);
 643 //printf("Resample_float::resample 2 %d %d %d\n", out_len, output_size[channel], channel);
 644
 645 //printf("Resample_float::resample 2 %d %d\n", channel, output_size[channel]);
 646
 647         return total_input;
 648 }