Reported by Fedora team for gcc-13 and Andrew created patch here
[goodguy/cinelerra.git] / 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