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++)
264 freq_real[i] = freq_real[size - i];
265 freq_imag[i] = -freq_imag[size - i];
274 CrossfadeFFT::CrossfadeFFT() : FFT()
280 CrossfadeFFT::~CrossfadeFFT()
285 int CrossfadeFFT::reset()
294 // samples in input_buffer and output_buffer
297 input_allocation = 0;
298 output_allocation = 0;
304 int CrossfadeFFT::delete_fft()
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;
316 int CrossfadeFFT::fix_window_size()
318 // fix the window size
319 // window size must be a power of 2
321 while(new_size < window_size) new_size *= 2;
322 window_size = MIN(131072, window_size);
323 window_size = new_size;
327 int CrossfadeFFT::initialize(int window_size)
329 this->window_size = window_size;
335 long CrossfadeFFT::get_delay()
337 return window_size + HALF_WINDOW;
340 int CrossfadeFFT::reconfigure()
350 // int CrossfadeFFT::process_fifo(long size,
351 // double *input_ptr,
352 // double *output_ptr)
354 // // Load next input buffer
355 // if(input_size + size > input_allocation)
357 // double *new_input = new double[input_size + size];
360 // memcpy(new_input, input_buffer, sizeof(double) * input_size);
361 // delete [] input_buffer;
363 // input_buffer = new_input;
364 // input_allocation = input_size + size;
367 // memcpy(input_buffer + input_size,
369 // size * sizeof(double));
370 // input_size += size;
378 // // Have enough to do some windows
379 // while(input_size >= window_size)
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];
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
395 // int result = signal_process();
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
408 // // Crossfade into the output buffer
409 // long new_allocation = output_size + window_size;
410 // if(new_allocation > output_allocation)
412 // double *new_output = new double[new_allocation];
416 // memcpy(new_output, output_buffer, sizeof(double) * output_size);
417 // delete [] output_buffer;
419 // output_buffer = new_output;
420 // output_allocation = new_allocation;
423 // if(output_size >= HALF_WINDOW)
425 // for(int i = 0, j = output_size - HALF_WINDOW;
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;
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;
441 // // First buffer has no crossfade
442 // memcpy(output_buffer + output_size,
444 // sizeof(double) * window_size);
445 // output_size += window_size;
449 // // Shift input buffer forward
450 // for(int i = window_size - HALF_WINDOW, j = 0;
453 // input_buffer[j] = input_buffer[i];
454 // input_size -= window_size - HALF_WINDOW;
460 // // Have enough to send to output
461 // int samples_rendered = 0;
462 // if(output_size - HALF_WINDOW >= size)
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;
472 // bzero(output_ptr, sizeof(double) * size);
475 // return samples_rendered;
480 int CrossfadeFFT::process_buffer(int64_t output_sample,
486 int step = (direction == PLAY_FORWARD) ? 1 : -1;
488 // User seeked so output buffer is invalid
489 if(output_sample != this->output_sample)
494 this->output_sample = output_sample;
495 this->input_sample = output_sample;
498 // Fill output buffer half a window at a time until size samples are available
499 while(output_size < size)
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];
507 // Fill enough input to make a window starting at output_sample
509 result = read_samples(this->input_sample,
514 input_buffer->set_offset(HALF_WINDOW);
515 // printf("CrossfadeFFT::process_buffer %d %lld %lld\n",
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,
522 input_buffer->set_offset(0);
525 input_size = window_size;
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
535 result = signal_process();
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
546 result = post_process();
548 // Allocate output buffer
549 int new_allocation = output_size + window_size;
550 if(new_allocation > output_allocation)
552 double *new_output = new double[new_allocation];
557 sizeof(double) * (output_size + HALF_WINDOW));
558 delete [] output_buffer;
560 output_buffer = new_output;
561 output_allocation = new_allocation;
564 // Overlay processed buffer
567 memcpy(output_buffer + output_size,
569 sizeof(double) * window_size);
574 for(int i = 0, j = output_size; i < HALF_WINDOW; i++, j++)
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;
582 //output_buffer[output_size] = 100.0;
583 //output_buffer[output_size + HALF_WINDOW] = -100.0;
585 memcpy(output_buffer + output_size + HALF_WINDOW,
586 output_real + HALF_WINDOW,
587 sizeof(double) * HALF_WINDOW);
591 output_size += HALF_WINDOW;
593 // Shift input buffer
594 double *input_samples = input_buffer->get_data();
595 for(int i = window_size - HALF_WINDOW, j = 0;
599 input_samples[j] = input_samples[i];
602 input_size = HALF_WINDOW;
603 this->input_sample += step * HALF_WINDOW;
608 // Transfer output buffer
611 memcpy(output_ptr->get_data(), output_buffer, sizeof(double) * size);
613 for(int i = 0, j = size; j < output_size + HALF_WINDOW; i++, j++)
614 output_buffer[i] = output_buffer[j];
616 this->output_sample += step * size;
617 this->output_size -= size;
623 int CrossfadeFFT::read_samples(int64_t output_sample,
630 int CrossfadeFFT::signal_process()
637 int CrossfadeFFT::post_process()