4 * Copyright (C) 1997-2011 Adam Williams <broadcast at earthling dot net>
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.
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.
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
29 #include "transportque.inc"
31 #define HALF_WINDOW (window_size / 2)
42 void FFT::bit_reverse(unsigned int samples,
43 double *real_in, double *imag_in,
44 double *real_out, double *imag_out)
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
53 for( i=1; i<samples1; ++i ) {
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;
60 imag_out[samples1] = 0.0;
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;
72 else { // Do simultaneous data copy and bit-reversal ordering
74 real_out[0] = real_in[0];
76 for( i=0; i<samples1; ++i ) {
77 j = reverse_bits(i, num_bits);
79 real_out[j] = real_in[i];
80 real_out[i] = real_in[j];
81 imag_out[i] = imag_out[j] = 0.0;
83 real_out[samples1] = real_in[samples1];
84 imag_out[samples1] = 0.0;
87 for( i=0; i<samples; ++i ) {
88 j = reverse_bits(i, num_bits);
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];
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
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;
116 if( samples >= 4 ) { // block 2, delta_angle = pi/2
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;
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;
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);
154 for( int i=0; i<samples; i+=bsz2 ) {
155 ar[2] = cm2; ar[1] = cm1;
156 ai[2] = sm2; ai[1] = sm1;
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];
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;
177 // Normalize if inverse transform
179 double scale = 1./samples;
180 for( int i=0; i<samples; ++i ) {
181 real_out[i] *= scale;
182 imag_out[i] *= scale;
189 int FFT::update_progress(int current_position)
194 unsigned int FFT::samples_to_bits(unsigned int samples)
196 if( !samples ) return ~0;
197 int i = 0; --samples;
198 while( (samples>>i) != 0 ) ++i;
203 unsigned int FFT::reverse_bits(unsigned int index, unsigned int bits)
207 for( i = rev = 0; i < bits; i++ ) {
208 rev = (rev << 1) | (index & 1);
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,
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,
237 unsigned int FFT::reverse_bits(unsigned int index, unsigned int bits)
240 union { unsigned int u; uint8_t b[sizeof(unsigned int)]; } data;
243 index = rev_bytes[data.b[0]] >> (8-bits);
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);
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);
259 int FFT::symmetry(int size, double *freq_real, double *freq_imag)
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];
273 CrossfadeFFT::CrossfadeFFT() : FFT()
280 CrossfadeFFT::~CrossfadeFFT()
285 int CrossfadeFFT::reset()
296 // samples in input_buffer and output_buffer
299 input_allocation = 0;
300 output_allocation = 0;
306 int CrossfadeFFT::delete_fft()
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;
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;
323 int CrossfadeFFT::fix_window_size()
325 // fix the window size
326 // window size must be a power of 2
328 while(new_size < window_size) new_size *= 2;
329 window_size = MIN(131072, window_size);
330 window_size = new_size;
334 int CrossfadeFFT::initialize(int window_size, int bands)
336 this->window_size = window_size;
343 long CrossfadeFFT::get_delay()
345 return window_size + HALF_WINDOW;
348 int CrossfadeFFT::reconfigure()
354 void CrossfadeFFT::allocate_output(int new_allocation)
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);
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];
370 output_buffer[i] = new_output;
372 output_allocation = new_allocation;
376 int CrossfadeFFT::process_buffer(int64_t output_sample, long size,
377 Samples *output_ptr, int direction)
379 Samples *output_temp[1];
380 output_temp[0] = output_ptr;
381 return process_buffer(output_sample, size, output_temp, direction);
384 // int CrossfadeFFT::get_read_offset()
388 int CrossfadeFFT::process_buffer(int64_t output_sample,
389 long size, Samples **output_ptr, int direction)
392 int step = (direction == PLAY_FORWARD) ? 1 : -1;
394 // User seeked so output buffer is invalid
395 if( output_sample != this->output_sample ) {
399 this->output_sample = output_sample;
400 input_sample = output_sample;
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);
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;
415 need_samples += HALF_WINDOW;
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);
422 memcpy(new_input_buffer->get_data(),
423 input_buffer->get_data(),
424 input_size * sizeof(double));
428 input_buffer = new_input_buffer;
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);
436 input_buffer->set_offset(0);
437 input_size = need_samples;
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];
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
464 // Overlay processed window on the output buffers
466 // direct copy the 1st window
467 memcpy(output_buffer[band] + output_size,
469 sizeof(double) * window_size);
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;
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);
490 output_size += HALF_WINDOW;
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;
500 // input_buffer->get_data()[j] = input_buffer->get_data()[i];
503 input_size -= HALF_WINDOW;
508 // Transfer output buffer if the user wants it
510 for( int band = 0; band < bands; band++ ) {
511 memcpy(output_ptr[band]->get_data(), output_buffer[band], sizeof(double) * size);
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));
521 this->output_sample += step * size;
522 this->output_size -= size;
528 int CrossfadeFFT::read_samples(int64_t output_sample,
529 int samples, Samples *buffer)
534 int CrossfadeFFT::signal_process()
541 int CrossfadeFFT::post_process()
547 int CrossfadeFFT::signal_process(int band)
555 int CrossfadeFFT::post_process(int band)