7 #include "transportque.inc"
9 #define HALF_WINDOW (window_size / 2)
12 // we need to do some trickery to get around of fftw thread unsafetyness
13 fftw_plan_desc *FFT::fftw_plans = 0;
14 Mutex FFT::plans_lock = Mutex();
24 int FFT::do_fft(unsigned int samples, // must be a power of 2
25 int inverse, // 0 = forward FFT, 1 = inverse
26 double *real_in, // array of input's real samples
27 double *imag_in, // array of input's imag samples
28 double *real_out, // array of output's reals
31 unsigned int num_bits; // Number of bits needed to store indices
32 unsigned int i, j, k, n;
33 unsigned int block_size, block_end;
35 double angle_numerator = 2.0 * M_PI;
36 double tr, ti; // temp real, temp imaginary
39 angle_numerator = -angle_numerator;
41 num_bits = samples_to_bits(samples);
43 // Do simultaneous data copy and bit-reversal ordering into outputs
45 for(i = 0; i < samples; i++)
47 j = reverse_bits(i, num_bits);
48 real_out[j] = real_in[i];
49 imag_out[j] = (imag_in == 0) ? 0.0 : imag_in[i];
63 for(block_size = 2; block_size <= samples; block_size <<= 1)
65 delta_angle = angle_numerator / (double)block_size;
66 sm2 = sin(-2 * delta_angle);
67 sm1 = sin(-delta_angle);
68 cm2 = cos(-2 * delta_angle);
69 cm1 = cos(-delta_angle);
72 for(i = 0; i < samples; i += block_size)
80 for(j = i, n = 0; n < block_end; j++, n++)
82 ar[0] = w * ar[1] - ar[2];
86 ai[0] = w * ai[1] - ai[2];
91 tr = ar[0] * real_out[k] - ai[0] * imag_out[k];
92 ti = ar[0] * imag_out[k] + ai[0] * real_out[k];
94 real_out[k] = real_out[j] - tr;
95 imag_out[k] = imag_out[j] - ti;
102 block_end = block_size;
105 // Normalize if inverse transform
109 double denom = (double)samples;
111 for (i = 0; i < samples; i++)
113 real_out[i] /= denom;
114 imag_out[i] /= denom;
120 int FFT::update_progress(int current_position)
125 unsigned int FFT::samples_to_bits(unsigned int samples)
131 if(samples & (1 << i))
137 unsigned int FFT::reverse_bits(unsigned int index, unsigned int bits)
141 for(i = rev = 0; i < bits; i++)
143 rev = (rev << 1) | (index & 1);
150 int FFT::symmetry(int size, double *freq_real, double *freq_imag)
153 for(int i = h + 1; i < size; i++)
155 freq_real[i] = freq_real[size - i];
156 freq_imag[i] = -freq_imag[size - i];
161 // Create a proper fftw plan to be used later
162 int FFT::ready_fftw(unsigned int samples)
164 // FFTW plan generation is not thread safe, so we have to take precausions
165 FFT::plans_lock.lock();
166 fftw_plan_desc *plan;
170 for (plan = fftw_plans; plan; plan = plan->next)
171 if (plan->samples == samples)
179 fftw_complex *temp_data = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * samples);
180 my_fftw_plan = new fftw_plan_desc; // we never discard this, since they are static
181 my_fftw_plan->samples = samples;
182 my_fftw_plan->plan_forward = fftw_plan_dft_1d(samples, temp_data, temp_data, FFTW_FORWARD, FFTW_ESTIMATE);
183 my_fftw_plan->plan_backward = fftw_plan_dft_1d(samples, temp_data, temp_data, FFTW_BACKWARD, FFTW_ESTIMATE);
184 // We will use this plan only in guru mode so we can now discard the temp_data
185 fftw_free(temp_data);
187 // Put the plan into the linked list
188 my_fftw_plan->next = fftw_plans;
189 fftw_plans = my_fftw_plan;
192 FFT::plans_lock.unlock();
196 int FFT::do_fftw_inplace(unsigned int samples,
201 fftw_execute_dft(my_fftw_plan->plan_forward, data, data);
203 fftw_execute_dft(my_fftw_plan->plan_backward, data, data);
209 CrossfadeFFT::CrossfadeFFT() : FFT()
215 CrossfadeFFT::~CrossfadeFFT()
220 int CrossfadeFFT::reset()
229 // samples in input_buffer and output_buffer
232 input_allocation = 0;
233 output_allocation = 0;
244 int CrossfadeFFT::delete_fft()
246 if(input_buffer) delete [] input_buffer;
247 if(output_buffer) delete [] output_buffer;
248 if(freq_real) delete [] freq_real;
249 if(freq_imag) delete [] freq_imag;
250 if(temp_real) delete [] temp_real;
251 if(temp_imag) delete [] temp_imag;
252 if(pre_window) delete [] pre_window;
253 if(post_window) delete [] post_window;
254 if(fftw_data) fftw_free(fftw_data);
259 int CrossfadeFFT::fix_window_size()
261 // fix the window size
262 // window size must be a power of 2
264 while(new_size < window_size) new_size *= 2;
265 window_size = MIN(131072, window_size);
266 window_size = new_size;
270 int CrossfadeFFT::initialize(int window_size)
272 this->window_size = window_size;
278 long CrossfadeFFT::get_delay()
280 return window_size + HALF_WINDOW;
283 int CrossfadeFFT::reconfigure()
295 int CrossfadeFFT::process_buffer(int64_t output_sample,
303 int step = (direction == PLAY_FORWARD) ? 1 : -1;
305 int overlap_size = window_size / oversample;
311 printf("ERROR, no output pointer!\n");
314 if(output_sample != this->output_sample || first_window)
318 this->output_sample = output_sample;
320 start_skip = window_size - overlap_size;
321 total_size = size + start_skip;
322 // signal_process() will always be able to know which window it has by looking at input_sample
323 this->input_sample = output_sample - step * start_skip;
324 if (step == -1) this->input_sample += overlap_size;
332 // Find out how big output buffer we will need, take overlapping into account
333 int new_allocation = total_size + window_size;
334 if(new_allocation > output_allocation)
336 double *new_output = new double[new_allocation];
341 sizeof(double) * (samples_ready + window_size - overlap_size));
342 delete [] output_buffer;
345 output_buffer = new_output;
346 output_allocation = new_allocation;
348 // Fill output buffer by overlap_size at a time until size samples are available
349 while(samples_ready < total_size)
351 if(!input_buffer) input_buffer = new double[window_size];
352 if(!fftw_data) fftw_data = (fftw_complex *)fftw_malloc(window_size * sizeof(fftw_complex));
354 // Fill enough input to make a window starting at output_sample
362 read_start = this->input_sample;
364 read_start = this->input_sample - window_size;
366 read_len = window_size;
371 read_start = this->input_sample + window_size - overlap_size;
372 write_pos = window_size - overlap_size;
375 read_start = this->input_sample - window_size;
378 read_len = overlap_size;
383 // completely outside the track
384 memset (input_buffer + write_pos, 0, read_len * sizeof(double));
387 if (read_start + read_len *step < 0)
389 // special case for reading before the track - in general it would be sensible that this behaviour is done by read_samples()
390 memset (input_buffer, 0, (read_len - read_start) * sizeof(double));
391 result = read_samples(read_start,
393 input_buffer + read_len *step + read_start + write_pos);
396 //printf("Readstart: %lli, read len: %i, write pos: %i\n", read_start, read_len, write_pos);
397 result = read_samples(read_start,
399 input_buffer + write_pos);
403 // apply Hanning window to input samples
404 for (int i = 0; i< window_size; i++)
406 fftw_data[i][0] = input_buffer[i] * pre_window[i];
412 do_fftw_inplace(window_size, 0, fftw_data);
414 result = signal_process();
415 // result = signal_process_oversample(first_window);
417 do_fftw_inplace(window_size, 1, fftw_data);
419 // Overlay over existing output - overlap processing
422 for (int i = 0; i < window_size - overlap_size; i++)
423 output_buffer[i + samples_ready] += fftw_data[i][0] * post_window[i];
424 for (int i = window_size - overlap_size; i < window_size; i++)
425 output_buffer[i + samples_ready] = fftw_data[i][0] * post_window[i];
428 int offset = output_allocation - samples_ready - window_size;
429 for (int i = 0; i < overlap_size; i++)
430 output_buffer[i + offset] = fftw_data[i][0] * post_window[i];
431 for (int i = overlap_size; i < window_size; i++)
432 output_buffer[i + offset] += fftw_data[i][0] * post_window[i];
436 // Shift input buffer
438 memmove(input_buffer, input_buffer + overlap_size, (window_size - overlap_size) * sizeof(double));
440 memmove(input_buffer + overlap_size, input_buffer, (window_size - overlap_size) * sizeof(double));
442 this->input_sample += step * overlap_size;
444 samples_ready += overlap_size;
450 memcpy(output_ptr, output_buffer + start_skip , size * sizeof(double));
451 samples_ready -= total_size;
453 memmove(output_buffer,
454 output_buffer + total_size,
455 (samples_ready + window_size - overlap_size) * sizeof(double));
456 this->output_sample += size;
460 memcpy(output_ptr, output_buffer + output_allocation - total_size , size * sizeof(double));
461 samples_ready -= total_size;
463 memmove(output_buffer + output_allocation - (samples_ready + window_size - overlap_size),
464 output_buffer + output_allocation - (samples_ready + window_size - overlap_size) - total_size,
465 (samples_ready + window_size - overlap_size) * sizeof(double));
467 this->output_sample -= size;
474 int step = (direction == PLAY_FORWARD) ? 1 : -1;
476 if(output_sample != this->output_sample || first_window)
481 this->output_sample = output_sample;
482 this->input_sample = output_sample;
485 // Fill output buffer half a window at a time until size samples are available
486 while(output_size < size)
488 if(!input_buffer) input_buffer = new double[window_size];
489 if(!freq_real) freq_real = new double[window_size];
490 if(!freq_imag) freq_imag = new double[window_size];
491 if(!temp_real) temp_real = new double[window_size];
492 if(!temp_imag) temp_imag = new double[window_size];
494 // Fill enough input to make a window starting at output_sample
496 result = read_samples(this->input_sample,
500 result = read_samples(this->input_sample + step * HALF_WINDOW,
502 input_buffer + HALF_WINDOW);
504 input_size = window_size;
507 do_fft(window_size, // must be a power of 2
508 0, // 0 = forward FFT, 1 = inverse
509 input_buffer, // array of input's real samples
510 0, // array of input's imag samples
511 freq_real, // array of output's reals
514 result = signal_process();
516 do_fft(window_size, // must be a power of 2
517 1, // 0 = forward FFT, 1 = inverse
518 freq_real, // array of input's real samples
519 freq_imag, // array of input's imag samples
520 temp_real, // array of output's reals
521 temp_imag); // array of output's imaginaries
523 // Allocate output buffer
524 int new_allocation = output_size + window_size;
525 if(new_allocation > output_allocation)
527 double *new_output = new double[new_allocation];
532 sizeof(double) * (output_size + HALF_WINDOW));
533 delete [] output_buffer;
535 output_buffer = new_output;
536 output_allocation = new_allocation;
539 // Overlay processed buffer
542 memcpy(output_buffer + output_size,
544 sizeof(double) * window_size);
549 for(int i = 0, j = output_size; i < HALF_WINDOW; i++, j++)
551 double src_level = (double)i / HALF_WINDOW;
552 double dst_level = (double)(HALF_WINDOW - i) / HALF_WINDOW;
553 output_buffer[j] = output_buffer[j] * dst_level +
554 temp_real[i] * src_level;
557 memcpy(output_buffer + output_size + HALF_WINDOW,
558 temp_real + HALF_WINDOW,
559 sizeof(double) * HALF_WINDOW);
562 output_size += HALF_WINDOW;
564 // Shift input buffer
565 for(int i = window_size - HALF_WINDOW, j = 0;
569 input_buffer[j] = input_buffer[i];
571 input_size = HALF_WINDOW;
572 this->input_sample += step * HALF_WINDOW;
577 // Transfer output buffer
580 memcpy(output_ptr, output_buffer, sizeof(double) * size);
582 for(int i = 0, j = size; j < output_size + HALF_WINDOW; i++, j++)
583 output_buffer[i] = output_buffer[j];
584 this->output_sample += step * size;
585 this->output_size -= size;
591 void CrossfadeFFT::set_oversample(int oversample)
593 // Only powers of two can be used for oversample
594 int oversample_fix = 2;
595 while(oversample_fix < oversample) oversample_fix *= 2;
596 this->oversample = oversample = oversample_fix;
598 // Precalculate the pre-envelope hanning window
599 pre_window = new double[window_size];
600 for (int i = 0; i< window_size; i++)
601 pre_window[i] = 0.5 - 0.5 *cos(2 * M_PI * i / window_size);
603 // Precalculate the post-envelope hanning window also, we could have triangle here also
604 post_window = new double[window_size];
605 /* for (int i = 0; i< window_size/2; i++)
606 post_window[i] = 1.0 * i / (window_size/2) / oversample * 2;
607 for (int i = window_size/2; i< window_size; i++)
608 post_window[i] = 1.0 * (window_size - i) / (window_size/2) / oversample * 2;
610 for (int i = 0; i< window_size; i++)
611 post_window[i] = (0.5 - 0.5 *cos(2 * M_PI * i / window_size)) * 6/ oversample / window_size;
613 ready_fftw(window_size);
617 void smbFft(double *fftBuffer, long fftFrameSize, long sign);
624 int CrossfadeFFT::read_samples(int64_t output_sample,
631 int CrossfadeFFT::signal_process()