Credit Andrew - improve in-tree documentation
[goodguy/cinelerra.git] / 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         {
264                 freq_real[i] = freq_real[size - i];
265                 freq_imag[i] = -freq_imag[size - i];
266         }
267         return 0;
268 }
269
270
271
272
273
274 CrossfadeFFT::CrossfadeFFT() : FFT()
275 {
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         output_real = 0;
292         output_imag = 0;
293         first_window = 1;
294 // samples in input_buffer and output_buffer
295         input_size = 0;
296         output_size = 0;
297         input_allocation = 0;
298         output_allocation = 0;
299         output_sample = 0;
300         input_sample = 0;
301         return 0;
302 }
303
304 int CrossfadeFFT::delete_fft()
305 {
306         if(input_buffer) delete input_buffer;
307         if(output_buffer) delete [] output_buffer;
308         if(freq_real) delete [] freq_real;
309         if(freq_imag) delete [] freq_imag;
310         if(output_real) delete [] output_real;
311         if(output_imag) delete [] output_imag;
312         reset();
313         return 0;
314 }
315
316 int CrossfadeFFT::fix_window_size()
317 {
318 // fix the window size
319 // window size must be a power of 2
320         int new_size = 16;
321         while(new_size < window_size) new_size *= 2;
322         window_size = MIN(131072, window_size);
323         window_size = new_size;
324         return 0;
325 }
326
327 int CrossfadeFFT::initialize(int window_size)
328 {
329         this->window_size = window_size;
330         first_window = 1;
331         reconfigure();
332         return 0;
333 }
334
335 long CrossfadeFFT::get_delay()
336 {
337         return window_size + HALF_WINDOW;
338 }
339
340 int CrossfadeFFT::reconfigure()
341 {
342         delete_fft();
343         fix_window_size();
344
345
346
347         return 0;
348 }
349
350 // int CrossfadeFFT::process_fifo(long size,
351 //      double *input_ptr,
352 //      double *output_ptr)
353 // {
354 // // Load next input buffer
355 //      if(input_size + size > input_allocation)
356 //      {
357 //              double *new_input = new double[input_size + size];
358 //              if(input_buffer)
359 //              {
360 //                      memcpy(new_input, input_buffer, sizeof(double) * input_size);
361 //                      delete [] input_buffer;
362 //              }
363 //              input_buffer = new_input;
364 //              input_allocation = input_size + size;
365 //      }
366 //
367 //      memcpy(input_buffer + input_size,
368 //              input_ptr,
369 //              size * sizeof(double));
370 //      input_size += size;
371 //
372 //
373 //
374 //
375 //
376 //
377 //
378 // // Have enough to do some windows
379 //      while(input_size >= window_size)
380 //      {
381 //              if(!freq_real) freq_real = new double[window_size];
382 //              if(!freq_imag) freq_imag = new double[window_size];
383 //              if(!output_real) output_real = new double[window_size];
384 //              if(!output_imag) output_imag = new double[window_size];
385 //
386 //
387 //
388 //              do_fft(window_size,  // must be a power of 2
389 //                      0,               // 0 = forward FFT, 1 = inverse
390 //                      input_buffer,    // array of input's real samples
391 //                      0,       // array of input's imag samples
392 //                      freq_real,      // array of output's reals
393 //                      freq_imag);
394 //
395 //              int result = signal_process();
396 //
397 //              if(!result)
398 //              {
399 //                      do_fft(window_size,  // must be a power of 2
400 //                              1,                         // 0 = forward FFT, 1 = inverse
401 //                              freq_real,       // array of input's real samples
402 //                              freq_imag,       // array of input's imag samples
403 //                              output_real,     // array of output's reals
404 //                              output_imag);
405 //              }
406 //
407 //
408 // // Crossfade into the output buffer
409 //              long new_allocation = output_size + window_size;
410 //              if(new_allocation > output_allocation)
411 //              {
412 //                      double *new_output = new double[new_allocation];
413 //
414 //                      if(output_buffer)
415 //                      {
416 //                              memcpy(new_output, output_buffer, sizeof(double) * output_size);
417 //                              delete [] output_buffer;
418 //                      }
419 //                      output_buffer = new_output;
420 //                      output_allocation = new_allocation;
421 //              }
422 //
423 //              if(output_size >= HALF_WINDOW)
424 //              {
425 //                      for(int i = 0, j = output_size - HALF_WINDOW;
426 //                              i < HALF_WINDOW;
427 //                              i++, j++)
428 //                      {
429 //                              double src_level = (double)i / HALF_WINDOW;
430 //                              double dst_level = (double)(HALF_WINDOW - i) / HALF_WINDOW;
431 //                              output_buffer[j] = output_buffer[j] * dst_level + output_real[i] * src_level;
432 //                      }
433 //
434 //                      memcpy(output_buffer + output_size,
435 //                              output_real + HALF_WINDOW,
436 //                              sizeof(double) * (window_size - HALF_WINDOW));
437 //                      output_size += window_size - HALF_WINDOW;
438 //              }
439 //              else
440 //              {
441 // // First buffer has no crossfade
442 //                      memcpy(output_buffer + output_size,
443 //                              output_real,
444 //                              sizeof(double) * window_size);
445 //                      output_size += window_size;
446 //              }
447 //
448 //
449 // // Shift input buffer forward
450 //              for(int i = window_size - HALF_WINDOW, j = 0;
451 //                      i < input_size;
452 //                      i++, j++)
453 //                      input_buffer[j] = input_buffer[i];
454 //              input_size -= window_size - HALF_WINDOW;
455 //      }
456 //
457 //
458 //
459 //
460 // // Have enough to send to output
461 //      int samples_rendered = 0;
462 //      if(output_size - HALF_WINDOW >= size)
463 //      {
464 //              memcpy(output_ptr, output_buffer, sizeof(double) * size);
465 //              for(int i = size, j = 0; i < output_size; i++, j++)
466 //                      output_buffer[j] = output_buffer[i];
467 //              output_size -= size;
468 //              samples_rendered = size;
469 //      }
470 //      else
471 //      {
472 //              bzero(output_ptr, sizeof(double) * size);
473 //      }
474 //
475 //      return samples_rendered;
476 // }
477
478
479
480 int CrossfadeFFT::process_buffer(int64_t output_sample,
481         long size,
482         Samples *output_ptr,
483         int direction)
484 {
485         int result = 0;
486         int step = (direction == PLAY_FORWARD) ? 1 : -1;
487
488 // User seeked so output buffer is invalid
489         if(output_sample != this->output_sample)
490         {
491                 output_size = 0;
492                 input_size = 0;
493                 first_window = 1;
494                 this->output_sample = output_sample;
495                 this->input_sample = output_sample;
496         }
497
498 // Fill output buffer half a window at a time until size samples are available
499         while(output_size < size)
500         {
501                 if(!input_buffer) input_buffer = new Samples(window_size);
502                 if(!freq_real) freq_real = new double[window_size];
503                 if(!freq_imag) freq_imag = new double[window_size];
504                 if(!output_real) output_real = new double[window_size];
505                 if(!output_imag) output_imag = new double[window_size];
506
507 // Fill enough input to make a window starting at output_sample
508                 if(first_window)
509                         result = read_samples(this->input_sample,
510                                 window_size,
511                                 input_buffer);
512                 else
513                 {
514                         input_buffer->set_offset(HALF_WINDOW);
515 // printf("CrossfadeFFT::process_buffer %d %lld %lld\n",
516 // __LINE__,
517 // this->input_sample + step * HALF_WINDOW,
518 // this->input_sample + step * HALF_WINDOW + HALF_WINDOW);
519                         result = read_samples(this->input_sample + step * HALF_WINDOW,
520                                 HALF_WINDOW,
521                                 input_buffer);
522                         input_buffer->set_offset(0);
523                 }
524
525                 input_size = window_size;
526
527                 if(!result)
528                         do_fft(window_size,   // must be a power of 2
529                                 0,                              // 0 = forward FFT, 1 = inverse
530                                 input_buffer->get_data(),        // array of input's real samples
531                                 0,                              // array of input's imag samples
532                                 freq_real,              // array of output's reals
533                                 freq_imag);
534                 if(!result)
535                         result = signal_process();
536
537                 if(!result)
538                         do_fft(window_size,  // must be a power of 2
539                                 1,                         // 0 = forward FFT, 1 = inverse
540                                 freq_real,         // array of input's real samples
541                                 freq_imag,         // array of input's imag samples
542                                 output_real,       // array of output's reals
543                                 output_imag);     // array of output's imaginaries
544
545                 if(!result)
546                         result = post_process();
547
548 // Allocate output buffer
549                 int new_allocation = output_size + window_size;
550                 if(new_allocation > output_allocation)
551                 {
552                         double *new_output = new double[new_allocation];
553                         if(output_buffer)
554                         {
555                                 memcpy(new_output,
556                                         output_buffer,
557                                         sizeof(double) * (output_size + HALF_WINDOW));
558                                 delete [] output_buffer;
559                         }
560                         output_buffer = new_output;
561                         output_allocation = new_allocation;
562                 }
563
564 // Overlay processed buffer
565                 if(first_window)
566                 {
567                         memcpy(output_buffer + output_size,
568                                 output_real,
569                                 sizeof(double) * window_size);
570                         first_window = 0;
571                 }
572                 else
573                 {
574                         for(int i = 0, j = output_size; i < HALF_WINDOW; i++, j++)
575                         {
576                                 double src_level = (double)i / HALF_WINDOW;
577                                 double dst_level = (double)(HALF_WINDOW - i) / HALF_WINDOW;
578                                 output_buffer[j] = output_buffer[j] * dst_level +
579                                         output_real[i] * src_level;
580                         }
581
582 //output_buffer[output_size] = 100.0;
583 //output_buffer[output_size + HALF_WINDOW] = -100.0;
584
585                         memcpy(output_buffer + output_size + HALF_WINDOW,
586                                 output_real + HALF_WINDOW,
587                                 sizeof(double) * HALF_WINDOW);
588                 }
589
590
591                 output_size += HALF_WINDOW;
592
593 // Shift input buffer
594                 double *input_samples = input_buffer->get_data();
595                 for(int i = window_size - HALF_WINDOW, j = 0;
596                         i < input_size;
597                         i++, j++)
598                 {
599                         input_samples[j] = input_samples[i];
600                 }
601
602                 input_size = HALF_WINDOW;
603                 this->input_sample += step * HALF_WINDOW;
604         }
605
606
607
608 // Transfer output buffer
609         if(output_ptr)
610         {
611                 memcpy(output_ptr->get_data(), output_buffer, sizeof(double) * size);
612         }
613         for(int i = 0, j = size; j < output_size + HALF_WINDOW; i++, j++)
614                 output_buffer[i] = output_buffer[j];
615
616         this->output_sample += step * size;
617         this->output_size -= size;
618
619         return 0;
620 }
621
622
623 int CrossfadeFFT::read_samples(int64_t output_sample,
624                 int samples,
625                 Samples *buffer)
626 {
627         return 1;
628 }
629
630 int CrossfadeFFT::signal_process()
631 {
632         return 0;
633 }
634
635
636
637 int CrossfadeFFT::post_process()
638 {
639         return 0;
640 }
641
642
643
644
645
646
647
648