cinelerra-5.1/cinelerra/fourier.C

   1
   2 /*
   3  * CINELERRA
   4  * Copyright (C) 1997-2011 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 <math.h>
  23 #include <stdio.h>
  24 #include <string.h>
  25
  26 #include "clip.h"
  27 #include "fourier.h"
  28 #include "samples.h"
  29 #include "transportque.inc"
  30
  31 #define HALF_WINDOW (window_size / 2)
  32
  33
  34 FFT::FFT()
  35 {
  36 }
  37
  38 FFT::~FFT()
  39 {
  40 }
  41
  42 void FFT::bit_reverse(unsigned int samples,
  43                 double *real_in, double *imag_in,
  44                 double *real_out, double *imag_out)
  45 {
  46         unsigned int i, j;
  47         unsigned int num_bits = samples_to_bits(samples);
  48         unsigned int samples1 = samples - 1;
  49         if( !real_out ) real_out = real_in;
  50         if( real_out == real_in ) {  // Do in-place bit-reversal ordering
  51                 if( !imag_in ) {
  52                         imag_out[0] = 0.0;
  53                         for( i=1; i<samples1; ++i ) {
  54                                 imag_out[i] = 0.0;
  55                                 j = reverse_bits(i, num_bits);
  56                                 if( i >= j ) continue;
  57                                 double tr = real_out[j];
  58                                 real_out[j] = real_out[i];  real_out[i] = tr;
  59                         }
  60                         imag_out[samples1] = 0.0;
  61                 }
  62                 else {
  63                         for( i=1; i<samples1; ++i ) {
  64                                 j = reverse_bits(i, num_bits);
  65                                 if( i >= j ) continue;
  66                                 double tr = real_out[j], ti = imag_out[j];
  67                                 real_out[j] = real_out[i];  real_out[i] = tr;
  68                                 imag_out[j] = imag_out[i];  imag_out[i] = ti;
  69                         }
  70                 }
  71         }
  72         else {  // Do simultaneous data copy and bit-reversal ordering
  73                 if( !imag_in ) {
  74                         real_out[0] = real_in[0];
  75                         imag_out[0] = 0.0;
  76                         for( i=0; i<samples1; ++i ) {
  77                                 j = reverse_bits(i, num_bits);
  78                                 if( i > j ) continue;
  79                                 real_out[j] = real_in[i];
  80                                 real_out[i] = real_in[j];
  81                                 imag_out[i] = imag_out[j] = 0.0;
  82                         }
  83                         real_out[samples1] = real_in[samples1];
  84                         imag_out[samples1] = 0.0;
  85                 }
  86                 else {
  87                         for( i=0; i<samples; ++i ) {
  88                                 j = reverse_bits(i, num_bits);
  89                                 if( i > j ) continue;
  90                                 real_out[i] = real_in[j];  imag_out[i] = imag_in[j];
  91                                 real_out[j] = real_in[i];  imag_out[j] = imag_in[i];
  92                         }
  93                 }
  94         }
  95 }
  96
  97 int FFT::do_fft(int samples,    // must be a power of 2
  98                 int inverse,    // 0 = forward FFT, 1 = inverse
  99                 double *real_in, double *imag_in, // complex input
 100                 double *real_out, double *imag_out) // complex output
 101 {
 102         bit_reverse(samples, real_in, imag_in, real_out, imag_out);
 103         double angle_numerator = !inverse ? 2*M_PI : -2*M_PI ;
 104         double ar[3], ai[3], rj, ij, rk, ik;
 105 // In the first two passes, all of the multiplys, and half the adds
 106 //  are saved by unrolling and reducing the x*1,x*0,x+0 operations.
 107         if( samples >= 2 ) { // block 1, delta_angle = pi
 108                 for( int i=0; i<samples; i+=2 ) {
 109                         int j = i, k = i+1; // cos(0)=1, sin(0)=0
 110                         rj = real_out[j];  ij = imag_out[j];
 111                         rk = real_out[k];  ik = imag_out[k];
 112                         real_out[j] += rk;        imag_out[j] += ik;
 113                         real_out[k] = rj - rk;  imag_out[k] = ij - ik;
 114                 }
 115         }
 116         if( samples >= 4 ) { // block 2, delta_angle = pi/2
 117                 if( !inverse ) {
 118                         for( int i=0; i<samples; i+=4 ) {
 119                                 int j = i, k = j+2; // cos(0)=1,sin(0)=0
 120                                 rj = real_out[j];  ij = imag_out[j];
 121                                 rk = real_out[k];  ik = imag_out[k];
 122                                 real_out[j] += rk;        imag_out[j] += ik;
 123                                 real_out[k] = rj - rk;  imag_out[k] = ij - ik;
 124                                 j = i+1;  k = j+2; // cos(-pi/2)=0, sin(-pi/2)=-1
 125                                 rj = real_out[j];  ij = imag_out[j];
 126                                 rk = real_out[k];  ik = imag_out[k];
 127                                 real_out[j] += ik;        imag_out[j] -= rk;
 128                                 real_out[k] = rj - ik;  imag_out[k] = ij + rk;
 129                         }
 130                 }
 131                 else {
 132                         for( int i=0; i<samples; i+=4 ) {
 133                                 int j = i, k = j+2; // cos(0)=1,sin(0)=0
 134                                 rj = real_out[j];  ij = imag_out[j];
 135                                 rk = real_out[k];  ik = imag_out[k];
 136                                 real_out[j] += rk;        imag_out[j] += ik;
 137                                 real_out[k] = rj - rk;  imag_out[k] = ij - ik;
 138                                 j = i+1;  k = j+2; // cos(pi/2)=0, sin(pi/2)=1
 139                                 rj = real_out[j];  ij = imag_out[j];
 140                                 rk = real_out[k];  ik = imag_out[k];
 141                                 real_out[j] -= ik;        imag_out[j] += rk;
 142                                 real_out[k] = rj + ik;  imag_out[k] = ij - rk;
 143                         }
 144                 }
 145         }
 146 // Do the rest of the FFT
 147         for( int bsz=4,bsz2=2*bsz; bsz2<=samples; bsz2<<=1 ) { // block 3,..
 148                 double delta_angle = angle_numerator / bsz2;
 149                 double sm2 = sin(2 * delta_angle);
 150                 double sm1 = sin(delta_angle);
 151                 double cm2 = cos(2 * delta_angle);
 152                 double cm1 = cos(delta_angle);
 153                 double w = 2 * cm1;
 154                 for( int i=0; i<samples; i+=bsz2 ) {
 155                         ar[2] = cm2;  ar[1] = cm1;
 156                         ai[2] = sm2;  ai[1] = sm1;
 157
 158                         double *rjp = real_out+i, *rkp = rjp + bsz;
 159                         double *ijp = imag_out+i, *ikp = ijp + bsz;
 160                         for( int n=bsz; --n>=0; ) {
 161                                 ar[0] = w * ar[1] - ar[2];
 162                                 ar[2] = ar[1];  ar[1] = ar[0];
 163                                 ai[0] = w * ai[1] - ai[2];
 164                                 ai[2] = ai[1];  ai[1] = ai[0];
 165
 166                                 double rj = *rjp, ij = *ijp;
 167                                 double rk = *rkp, ik = *ikp;
 168                                 double tr = ar[0] * rk - ai[0] * ik;
 169                                 double ti = ar[0] * ik + ai[0] * rk;
 170                                 *rjp++ += tr;     *ijp++ += ti;
 171                                 *rkp++ = rj - tr;  *ikp++ = ij - ti;
 172                         }
 173                 }
 174
 175                 bsz = bsz2;
 176         }
 177 // Normalize if inverse transform
 178         if( inverse ) {
 179                 double scale = 1./samples;
 180                 for( int i=0; i<samples; ++i ) {
 181                         real_out[i] *= scale;
 182                         imag_out[i] *= scale;
 183                 }
 184         }
 185
 186         return 0;
 187 }
 188
 189 int FFT::update_progress(int current_position)
 190 {
 191         return 0;
 192 }
 193
 194 unsigned int FFT::samples_to_bits(unsigned int samples)
 195 {
 196         if( !samples ) return ~0;
 197         int i = 0;  --samples;
 198         while( (samples>>i) != 0 ) ++i;
 199         return i;
 200 }
 201
 202 #if 0
 203 unsigned int FFT::reverse_bits(unsigned int index, unsigned int bits)
 204 {
 205         unsigned int i, rev;
 206
 207         for( i = rev = 0; i < bits; i++ ) {
 208                 rev = (rev << 1) | (index & 1);
 209                 index >>= 1;
 210         }
 211
 212         return rev;
 213 }
 214
 215 #else
 216
 217 uint8_t FFT::rev_bytes[256] = {
 218   0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0, 0x10, 0x90, 0x50, 0xd0, 0x30, 0xb0, 0x70, 0xf0,
 219   0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8, 0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8,
 220   0x04, 0x84, 0x44, 0xc4, 0x24, 0xa4, 0x64, 0xe4, 0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4,
 221   0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec, 0x1c, 0x9c, 0x5c, 0xdc, 0x3c, 0xbc, 0x7c, 0xfc,
 222   0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2, 0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2,
 223   0x0a, 0x8a, 0x4a, 0xca, 0x2a, 0xaa, 0x6a, 0xea, 0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa,
 224   0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6, 0x16, 0x96, 0x56, 0xd6, 0x36, 0xb6, 0x76, 0xf6,
 225   0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee, 0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe,
 226
 227   0x01, 0x81, 0x41, 0xc1, 0x21, 0xa1, 0x61, 0xe1, 0x11, 0x91, 0x51, 0xd1, 0x31, 0xb1, 0x71, 0xf1,
 228   0x09, 0x89, 0x49, 0xc9, 0x29, 0xa9, 0x69, 0xe9, 0x19, 0x99, 0x59, 0xd9, 0x39, 0xb9, 0x79, 0xf9,
 229   0x05, 0x85, 0x45, 0xc5, 0x25, 0xa5, 0x65, 0xe5, 0x15, 0x95, 0x55, 0xd5, 0x35, 0xb5, 0x75, 0xf5,
 230   0x0d, 0x8d, 0x4d, 0xcd, 0x2d, 0xad, 0x6d, 0xed, 0x1d, 0x9d, 0x5d, 0xdd, 0x3d, 0xbd, 0x7d, 0xfd,
 231   0x03, 0x83, 0x43, 0xc3, 0x23, 0xa3, 0x63, 0xe3, 0x13, 0x93, 0x53, 0xd3, 0x33, 0xb3, 0x73, 0xf3,
 232   0x0b, 0x8b, 0x4b, 0xcb, 0x2b, 0xab, 0x6b, 0xeb, 0x1b, 0x9b, 0x5b, 0xdb, 0x3b, 0xbb, 0x7b, 0xfb,
 233   0x07, 0x87, 0x47, 0xc7, 0x27, 0xa7, 0x67, 0xe7, 0x17, 0x97, 0x57, 0xd7, 0x37, 0xb7, 0x77, 0xf7,
 234   0x0f, 0x8f, 0x4f, 0xcf, 0x2f, 0xaf, 0x6f, 0xef, 0x1f, 0x9f, 0x5f, 0xdf, 0x3f, 0xbf, 0x7f, 0xff,
 235 };
 236
 237 unsigned int FFT::reverse_bits(unsigned int index, unsigned int bits)
 238 {
 239   unsigned char b;
 240   union { unsigned int u; uint8_t b[sizeof(unsigned int)]; } data;
 241   data.u = index;
 242   if( bits <= 8 ) {
 243         index = rev_bytes[data.b[0]] >> (8-bits);
 244   }
 245   else if( bits <= 16 ) {
 246         b = data.b[0];  data.b[0] = rev_bytes[data.b[1]];  data.b[1] = rev_bytes[b];
 247         index = data.u >> (16-bits);
 248   }
 249   else {
 250         b = data.b[0];  data.b[0] = rev_bytes[data.b[3]];  data.b[3] = rev_bytes[b];
 251         b = data.b[1];  data.b[1] = rev_bytes[data.b[2]];  data.b[2] = rev_bytes[b];
 252         index = data.u >> (32-bits);
 253   }
 254   return index;
 255 }
 256
 257 #endif
 258
 259 int FFT::symmetry(int size, double *freq_real, double *freq_imag)
 260 {
 261         int h = size / 2;
 262         for( int i = h + 1; i < size; i++ ) {
 263                 freq_real[i] = freq_real[size - i];
 264                 freq_imag[i] = -freq_imag[size - i];
 265         }
 266         return 0;
 267 }
 268
 269
 270
 271
 272
 273 CrossfadeFFT::CrossfadeFFT() : FFT()
 274 {
 275         bands = 1;
 276         reset();
 277         window_size = 4096;
 278 }
 279
 280 CrossfadeFFT::~CrossfadeFFT()
 281 {
 282         delete_fft();
 283 }
 284
 285 int CrossfadeFFT::reset()
 286 {
 287         input_buffer = 0;
 288         output_buffer = 0;
 289         freq_real = 0;
 290         freq_imag = 0;
 291         freq_real2 = 0;
 292         freq_imag2 = 0;
 293         output_real = 0;
 294         output_imag = 0;
 295         first_window = 1;
 296 // samples in input_buffer and output_buffer
 297         input_size = 0;
 298         output_size = 0;
 299         input_allocation = 0;
 300         output_allocation = 0;
 301         output_sample = 0;
 302         input_sample = 0;
 303         return 0;
 304 }
 305
 306 int CrossfadeFFT::delete_fft()
 307 {
 308         delete input_buffer;      input_buffer = 0;
 309         if( output_buffer ) {
 310                 for( int i=0; i<bands; ++i ) delete [] output_buffer[i];
 311                 delete [] output_buffer;  output_buffer = 0;
 312         }
 313         delete [] freq_real;      freq_real = 0;
 314         delete [] freq_imag;      freq_imag = 0;
 315         delete [] freq_real2;    freq_real2 = 0;
 316         delete [] freq_imag2;    freq_imag2 = 0;
 317         delete [] output_real;  output_real = 0;
 318         delete [] output_imag;  output_imag = 0;
 319         reset();
 320         return 0;
 321 }
 322
 323 int CrossfadeFFT::fix_window_size()
 324 {
 325 // fix the window size
 326 // window size must be a power of 2
 327         int new_size = 16;
 328         while(new_size < window_size) new_size *= 2;
 329         window_size = MIN(131072, window_size);
 330         window_size = new_size;
 331         return 0;
 332 }
 333
 334 int CrossfadeFFT::initialize(int window_size, int bands)
 335 {
 336         this->window_size = window_size;
 337         this->bands = bands;
 338         first_window = 1;
 339         reconfigure();
 340         return 0;
 341 }
 342
 343 long CrossfadeFFT::get_delay()
 344 {
 345         return window_size + HALF_WINDOW;
 346 }
 347
 348 int CrossfadeFFT::reconfigure()
 349 {
 350         delete_fft();
 351         fix_window_size();
 352         return 0;
 353 }
 354 void CrossfadeFFT::allocate_output(int new_allocation)
 355 {
 356 // Allocate output buffer
 357         if( new_allocation > output_allocation ) {
 358                 if( !output_buffer ) {
 359                         output_buffer = new double*[bands];
 360                         bzero(output_buffer, sizeof(double) * bands);
 361                 }
 362
 363                 for( int i = 0; i < bands; i++ ) {
 364                         double *new_output = new double[new_allocation];
 365                         if( output_buffer[i] ) {
 366                                 memcpy(new_output, output_buffer[i],
 367                                         sizeof(double) * (output_size + HALF_WINDOW));
 368                                 delete [] output_buffer[i];
 369                         }
 370                         output_buffer[i] = new_output;
 371                 }
 372                 output_allocation = new_allocation;
 373         }
 374 }
 375
 376 int CrossfadeFFT::process_buffer(int64_t output_sample, long size,
 377                 Samples *output_ptr, int direction)
 378 {
 379         Samples *output_temp[1];
 380         output_temp[0] = output_ptr;
 381         return process_buffer(output_sample, size, output_temp, direction);
 382 }
 383
 384 // int CrossfadeFFT::get_read_offset()
 385 // {
 386 // }
 387
 388 int CrossfadeFFT::process_buffer(int64_t output_sample,
 389                 long size, Samples **output_ptr, int direction)
 390 {
 391         int result = 0;
 392         int step = (direction == PLAY_FORWARD) ? 1 : -1;
 393
 394 // User seeked so output buffer is invalid
 395         if( output_sample != this->output_sample ) {
 396                 output_size = 0;
 397                 input_size = 0;
 398                 first_window = 1;
 399                 this->output_sample = output_sample;
 400                 input_sample = output_sample;
 401         }
 402
 403
 404 // printf("CrossfadeFFT::process_buffer %d size=%ld input_size=%ld output_size=%ld window_size=%ld\n",
 405 // __LINE__, size, input_size, output_size, window_size);
 406
 407 // must call read_samples once so the upstream plugins don't have to seek
 408 // must be a multiple of 1/2 window
 409         int need_samples = (size - output_size) / HALF_WINDOW * HALF_WINDOW;
 410 // round up a half window
 411         if( need_samples + output_size < size ) {
 412                 need_samples += HALF_WINDOW;
 413         }
 414 // half window tail
 415         need_samples += HALF_WINDOW;
 416
 417 // extend the buffer to need_samples
 418         if( !input_buffer || input_buffer->get_allocated() < need_samples ) {
 419                 Samples *new_input_buffer = new Samples(need_samples);
 420
 421                 if( input_buffer ) {
 422                         memcpy(new_input_buffer->get_data(),
 423                                 input_buffer->get_data(),
 424                                 input_size * sizeof(double));
 425                         delete input_buffer;
 426                 }
 427
 428                 input_buffer = new_input_buffer;
 429         }
 430
 431         input_buffer->set_offset(input_size);
 432         if( need_samples > input_size ) {
 433                 result = read_samples(input_sample, need_samples-input_size, input_buffer);
 434                 input_sample += step * (need_samples - input_size);
 435         }
 436         input_buffer->set_offset(0);
 437         input_size = need_samples;
 438
 439
 440         if( !freq_real ) freq_real = new double[window_size];
 441         if( !freq_imag ) freq_imag = new double[window_size];
 442         if( !freq_real2 ) freq_real2 = new double[window_size];
 443         if( !freq_imag2 ) freq_imag2 = new double[window_size];
 444         if( !output_real ) output_real = new double[window_size];
 445         if( !output_imag ) output_imag = new double[window_size];
 446
 447 // Fill output buffer half a window at a time until size samples are available
 448         while( !result && output_size < size ) {
 449                 do_fft(window_size, 0,  // forward
 450                         input_buffer->get_data(), 0, // input, real only
 451                         freq_real2, freq_imag2); // output
 452                 allocate_output(output_size + window_size);
 453 // process & overlay each band separately
 454                 for( int band=0; band<bands; ++band ) {
 455 // restore from the backup
 456                         memcpy(freq_real, freq_real2, sizeof(double) * window_size);
 457                         memcpy(freq_imag, freq_imag2, sizeof(double) * window_size);
 458                         signal_process(band);
 459                         do_fft(window_size, 1,  // inverse
 460                                 freq_real, freq_imag, // input
 461                                 output_real, output_imag); // output
 462                         post_process(band);
 463
 464 // Overlay processed window on the output buffers
 465                         if( first_window ) {
 466 // direct copy the 1st window
 467                                 memcpy(output_buffer[band] + output_size,
 468                                         output_real,
 469                                         sizeof(double) * window_size);
 470                         }
 471                         else {
 472 // dissolve 1st half of later windows
 473                                 for( int i = 0, j = output_size; i < HALF_WINDOW; i++, j++ ) {
 474                                         double src_level = (double)i / HALF_WINDOW;
 475                                         double dst_level = (double)(HALF_WINDOW - i) / HALF_WINDOW;
 476                                         output_buffer[band][j] = output_buffer[band][j] * dst_level +
 477                                                 output_real[i] * src_level;
 478                                 }
 479
 480 //output_buffer[output_size] = 100.0;
 481 //output_buffer[output_size + HALF_WINDOW] = -100.0;
 482 // copy 2nd half of window
 483                                 memcpy(output_buffer[band] + output_size + HALF_WINDOW,
 484                                         output_real + HALF_WINDOW,
 485                                         sizeof(double) * HALF_WINDOW);
 486                         }
 487                 }
 488
 489                 first_window = 0;
 490                 output_size += HALF_WINDOW;
 491
 492 // Shift input buffer half a window forward
 493                 memcpy(input_buffer->get_data(),
 494                         input_buffer->get_data() + HALF_WINDOW,
 495                         (input_size - HALF_WINDOW) * sizeof(double));
 496 //              for( int i = HALF_WINDOW, j = 0;
 497 //                      i < input_size;
 498 //                      i++, j++ )
 499 //              {
 500 //                      input_buffer->get_data()[j] = input_buffer->get_data()[i];
 501 //              }
 502
 503                 input_size -= HALF_WINDOW;
 504         }
 505
 506
 507
 508 // Transfer output buffer if the user wants it
 509         if( output_ptr ) {
 510                 for( int band = 0; band < bands; band++ ) {
 511                         memcpy(output_ptr[band]->get_data(), output_buffer[band], sizeof(double) * size);
 512                 }
 513         }
 514
 515 // shift output buffers forward
 516         for( int band = 0; band < bands; band++ ) {
 517                 memcpy(output_buffer[band], output_buffer[band] + size,
 518                         sizeof(double) * (output_size + HALF_WINDOW - size));
 519         }
 520
 521         this->output_sample += step * size;
 522         this->output_size -= size;
 523
 524         return 0;
 525 }
 526
 527
 528 int CrossfadeFFT::read_samples(int64_t output_sample,
 529                 int samples, Samples *buffer)
 530 {
 531         return 1;
 532 }
 533
 534 int CrossfadeFFT::signal_process()
 535 {
 536         return 0;
 537 }
 538
 539
 540
 541 int CrossfadeFFT::post_process()
 542 {
 543         return 0;
 544 }
 545
 546
 547 int CrossfadeFFT::signal_process(int band)
 548 {
 549         signal_process();
 550         return 0;
 551 }
 552
 553
 554
 555 int CrossfadeFFT::post_process(int band)
 556 {
 557         post_process();
 558         return 0;
 559 }
 560
 561