cinelerra-5.1/cinelerra/resample.C

   1
   2 /*
   3  * CINELERRA
   4  * Copyright (C) 2009 Adam Williams <broadcast at earthling dot net>
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 2 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19  *
  20  */
  21
  22 #include "bcsignals.h"
  23 #include "clip.h"
  24 #include "resample.h"
  25 #include "samples.h"
  26 #include "transportque.inc"
  27
  28 #include <math.h>
  29 #include <stdio.h>
  30 #include <stdlib.h>
  31 #include <string.h>
  32
  33 // Resampling from Lame
  34
  35 Resample::Resample()
  36 {
  37         old.allocate(BLACKSIZE, 0);
  38         double *old_data = old.get_data();
  39         memset(old_data, 0, BLACKSIZE*sizeof(*old_data));
  40         resample_init = 0;
  41         last_ratio = 0;
  42         output_temp = 0;
  43         output_size = 0;
  44         output_allocation = 0;
  45         input_size = RESAMPLE_CHUNKSIZE;
  46         input_position = 0;
  47         input = new Samples(input_size + 1);
  48         output_position = 0;
  49         itime = 0;
  50         direction = PLAY_FORWARD;
  51 }
  52
  53
  54 Resample::~Resample()
  55 {
  56         delete [] output_temp;
  57         delete input;
  58 }
  59
  60 int Resample::read_samples(Samples *buffer, int64_t start, int64_t len, int direction)
  61 {
  62         return 0;
  63 }
  64
  65 void Resample::reset()
  66 {
  67         resample_init = 0;
  68         output_size = 0;
  69         output_position = 0;
  70         input_position = 0;
  71 }
  72
  73 /* This algorithm from:
  74  * SIGNAL PROCESSING ALGORITHMS IN FORTRAN AND C
  75  * S.D. Stearns and R.A. David, Prentice-Hall, 1992 */
  76 void Resample::blackman(double fcn, int filter_l)
  77 {
  78         double wcn = M_PI * fcn;
  79         double ctr = filter_l / 2.0;
  80         double cir = 2*M_PI/filter_l;
  81         for( int j=0; j<=2*BPC; ++j ) {
  82                 double offset = (j-BPC) / (2.*BPC);  // -0.5 ... 0.5
  83                 for( int i=0; i<=filter_l; ++i ) {
  84                         double x = i - offset;
  85                         bclamp(x, 0,filter_l);
  86                         double v, dx = x - ctr;
  87                         if( fabs(dx) >= 1e-9 ) {
  88                                 double curve = sin(wcn * dx) / (M_PI * dx);
  89                                 double th = x * cir;
  90                                 double blkmn = 0.42 - 0.5 * cos(th) + 0.08 * cos(2*th);
  91                                 v = blkmn * curve;
  92                         }
  93                         else
  94                                 v = fcn;
  95                         blackfilt[j][i] = v;
  96                 }
  97         }
  98 }
  99
 100
 101 int Resample::get_output_size()
 102 {
 103         return output_size;
 104 }
 105
 106 // void Resample::read_output(double *output, int size)
 107 // {
 108 //      memcpy(output, output_temp, size * sizeof(double));
 109 // // Shift leftover forward
 110 //      for( int i = size; i < output_size; i++ )
 111 //              output_temp[i - size] = output_temp[i];
 112 //      output_size -= size;
 113 // }
 114
 115
 116 // starts odd = (even-1)
 117 #define FILTER_N (BLACKSIZE-6)
 118 #define FILTER_L (FILTER_N - (~FILTER_N & 1));
 119
 120 void Resample::resample_chunk(Samples *input_buffer, int64_t in_len,
 121         int in_rate, int out_rate)
 122 {
 123 //printf("Resample::resample_chunk %d in_len=%jd input_size=%d\n",
 124 // __LINE__, in_len, input_size);
 125         double *input = input_buffer->get_data();
 126         double resample_ratio = (double)in_rate / out_rate;
 127         double fcn = .90 / resample_ratio;
 128         if( fcn > .90 ) fcn = .90;
 129         int filter_l = FILTER_L;
 130 // if resample_ratio = int, filter_l should include right edge
 131         if( fabs(resample_ratio - floor(.5 + resample_ratio)) < .0001 )
 132                 ++filter_l;
 133 // Blackman filter initialization must be called whenever there is a
 134 // sampling ratio change
 135         if( !resample_init || last_ratio != resample_ratio ) {
 136                 resample_init = 1;
 137                 last_ratio = resample_ratio;
 138                 blackman(fcn, filter_l);
 139                 itime = 0;
 140         }
 141
 142         double filter_l2 = filter_l/2.;
 143         int l2 = filter_l2;
 144         int64_t end_time = itime + in_len + l2;
 145         int64_t out_time = end_time / resample_ratio + 1;
 146         int64_t demand = out_time - output_position;
 147         if( demand >= output_allocation ) {
 148 // demand 2**n buffer
 149                 int64_t new_allocation = output_allocation ? output_allocation : 16384;
 150                 while( new_allocation < demand ) new_allocation <<= 1;
 151                 double *new_output = new double[new_allocation];
 152                 if( output_temp ) {
 153                         memmove(new_output, output_temp, output_allocation*sizeof(double));
 154                         delete [] output_temp;
 155                 }
 156                 output_temp = new_output;
 157                 output_allocation = new_allocation;
 158         }
 159
 160 // Main loop
 161         double *old_data = old.get_data();
 162         double ctr_pos = 0;
 163         int otime = 0, last_used = 0;
 164         while( output_size < output_allocation ) {
 165                 double in_pos = otime * resample_ratio;
 166 // window centered at ctr_pos
 167                 ctr_pos = in_pos + itime;
 168                 double pos = ctr_pos - filter_l2;
 169                 int ipos = floor(pos);
 170                 last_used = ipos + filter_l;
 171                 if( last_used >= in_len ) break;
 172                 double fraction = pos - ipos;
 173                 int phase = floor(fraction * 2*BPC + .5);
 174                 int i = ipos, j = filter_l;  // fir filter
 175                 double xvalue = 0, *filt = blackfilt[phase];
 176                 for( ; j>=0 && i<0; ++i,--j ) xvalue += *filt++ * old_data[BLACKSIZE + i];
 177                 for( ; j>=0;        ++i,--j ) xvalue += *filt++ * input[i];
 178                 output_temp[output_size++] = xvalue;
 179                 ++otime;
 180         }
 181 // move ctr_pos backward by in_len as new itime offset
 182 // the next read will be in the history, itime is negative
 183         itime = ctr_pos - in_len;
 184         memmove(old_data, input+in_len-BLACKSIZE, BLACKSIZE*sizeof(double));
 185 }
 186
 187 void Resample::reverse_buffer(double *buffer, int64_t len)
 188 {
 189         double *ap = buffer;
 190         double *bp = ap + len;
 191         while( ap < --bp ) {
 192                 double t = *ap;
 193                 *ap++ = *bp;
 194                 *bp = t;
 195         }
 196 }
 197
 198 int Resample::set_input_position(int64_t in_pos, int in_dir)
 199 {
 200         reset();
 201         input_position = in_pos;
 202         direction = in_dir;
 203 // update old, just before/after input going fwd/rev;
 204         int dir = direction == PLAY_FORWARD ? -1 : 1;
 205         in_pos += dir * BLACKSIZE;
 206         return read_samples(&old, in_pos, BLACKSIZE, in_dir);
 207 }
 208
 209 int Resample::resample(Samples *output, int64_t out_len,
 210                 int in_rate, int out_rate, int64_t out_position, int direction)
 211 {
 212         int result = 0;
 213         if( this->output_position != out_position ||
 214             this->direction != direction ) {
 215 //printf("missed  %jd!=%jd\n", output_position, out_position);
 216 // starting point in input rate.
 217                 int64_t in_pos = out_position * in_rate / out_rate;
 218                 set_input_position(in_pos, direction);
 219         }
 220 //else
 221 //printf("matched %jd==%jd\n", output_position, out_position);
 222
 223         int dir = direction == PLAY_REVERSE ? -1 : 1;
 224         int remaining_len = out_len;
 225         double *output_ptr = output->get_data();
 226         while( remaining_len > 0 && !result ) {
 227                 if( output_size ) {
 228                         int len = bmin(output_size, remaining_len);
 229                         memmove(output_ptr, output_temp, len*sizeof(double));
 230                         memmove(output_temp, output_temp+len, (output_size-=len)*sizeof(double));
 231                         output_ptr += len;  remaining_len -= len;
 232                 }
 233                 if( remaining_len > 0 ) {
 234                         result = read_samples(input, input_position, input_size, direction);
 235                         if( result ) break;
 236                         resample_chunk(input, input_size, in_rate, out_rate);
 237                         input_position += dir * input_size;
 238                 }
 239         }
 240         if( !result )
 241                 this->output_position = out_position + dir * out_len;
 242         return result;
 243 }
 244
 245
 246