4 * Copyright (C) 2008 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
22 #include "bcdisplayinfo.h"
38 #define WINDOW_BORDER (window_size / 2)
39 #define SGN(x) (x<0 ? -1: 1)
42 REGISTER_PLUGIN(DenoiseEffect)
48 DenoiseEffect::DenoiseEffect(PluginServer *server)
49 : PluginAClient(server)
55 DenoiseEffect::~DenoiseEffect()
62 LOAD_CONFIGURATION_MACRO(DenoiseEffect, DenoiseConfig)
64 NEW_WINDOW_MACRO(DenoiseEffect, DenoiseWindow)
66 void DenoiseEffect::delete_dsp()
68 if(ex_coeff_d) delete ex_coeff_d;
69 if(ex_coeff_r) delete ex_coeff_r;
70 if(ex_coeff_rn) delete ex_coeff_rn;
71 if(wave_coeff_d) delete wave_coeff_d;
72 if(wave_coeff_r) delete wave_coeff_r;
73 if(decomp_filter) delete decomp_filter;
74 if(recon_filter) delete recon_filter;
75 if(input_buffer) delete [] input_buffer;
76 if(output_buffer) delete [] output_buffer;
77 if(dsp_in) delete [] dsp_in;
78 if(dsp_out) delete [] dsp_out;
79 if(dsp_iteration) delete [] dsp_iteration;
96 void DenoiseEffect::reset()
111 input_allocation = 0;
112 output_allocation = 0;
129 const char* DenoiseEffect::plugin_title() { return N_("Denoise"); }
130 int DenoiseEffect::is_realtime() { return 1; }
134 void DenoiseEffect::read_data(KeyFrame *keyframe)
137 input.set_shared_input(keyframe->xbuf);
142 result = input.read_tag();
146 if(input.tag.title_is("DENOISE"))
148 config.level = input.tag.get_property("LEVEL", config.level);
154 void DenoiseEffect::save_data(KeyFrame *keyframe)
157 output.set_shared_output(keyframe->xbuf);
159 output.tag.set_title("DENOISE");
160 output.tag.set_property("LEVEL", config.level);
162 output.tag.set_title("/DENOISE");
164 output.append_newline();
165 output.terminate_string();
169 void DenoiseEffect::update_gui()
173 ((DenoiseWindow*)thread->window)->lock_window();
174 ((DenoiseWindow*)thread->window)->update();
175 ((DenoiseWindow*)thread->window)->unlock_window();
181 double DenoiseEffect::dot_product(double *data, double *filter, char filtlen)
187 for(i = 0; i < filtlen; i++) sum += *data-- * *filter++;
191 int DenoiseEffect::convolve_dec_2(double *input_sequence,
195 double *output_sequence)
197 // convolve the input sequence with the filter and decimate by two
198 int i, shortlen, offset;
199 int64_t lengthp4 = length + 4;
200 int64_t lengthm4 = length - 4;
201 int64_t lengthp5 = length + 5;
202 int64_t lengthp8 = length + 8;
204 for(i = 0; (i <= lengthp8) && ((i - filtlen) <= lengthp8); i += 2)
207 *output_sequence++ = dot_product(input_sequence + i, filter, i + 1);
211 offset = i - lengthm4;
212 shortlen = filtlen - offset;
213 *output_sequence++ = dot_product(input_sequence + lengthp4,
214 filter + offset, shortlen);
217 *output_sequence++ = dot_product(input_sequence + i, filter, filtlen);
222 int64_t DenoiseEffect::decompose_branches(double *in_data,
224 WaveletFilters *decomp_filter,
228 // Take input data and filters and form two branches of half the
229 // original length. Length of branches is returned.
230 convolve_dec_2(in_data, length, decomp_filter->h, decomp_filter->length, out_low);
231 convolve_dec_2(in_data, length, decomp_filter->g, decomp_filter->length, out_high);
235 int DenoiseEffect::wavelet_decomposition(double *in_data,
239 for(int i = 0; i < levels; i++)
241 in_length = decompose_branches(in_data,
245 out_data[(2 * i) + 1]);
247 in_data = out_data[2 * i];
252 int DenoiseEffect::tree_copy(double **output,
259 for(i = 0, k = 1; k < levels; i++, k++)
265 for(j = 0; j < length + 5; j++)
268 output[m][j] = input[m][j];
276 for(j = 0; j < length + 5; j++)
278 output[l][j] = input[l][j];
279 output[m][j] = input[m][j];
284 int DenoiseEffect::threshold(int window_size, double gammas, int levels)
287 double threshold, cv, cvb, abs_coeff_r;
288 double *coeff_r, *coeff_l;
291 for(i = 0; i < levels; i++)
293 length = (window_size >> (i + 1)) + 5;
294 threshold = sqrt(2 * log(length) / log(2)) * gammas / sqrt(length);
296 for(j = 0; j < length; j++)
298 coeff_r = &(ex_coeff_r->values[(2 * i) + 1][j]);
299 coeff_l = &(ex_coeff_rn->values[(2 * i) + 1][j]);
302 abs_coeff_r = fabs(*coeff_r);
303 cvb = abs_coeff_r - threshold;
306 if(abs_coeff_r > threshold)
322 double DenoiseEffect::dot_product_even(double *data, double *filter, int filtlen)
328 for(i = 0; i < filtlen; i += 2) sum += *data-- * filter[i];
333 double DenoiseEffect::dot_product_odd(double *data, double *filter, int filtlen)
339 for(i = 1; i < filtlen; i += 2) sum += *data-- * filter[i];
343 int DenoiseEffect::convolve_int_2(double *input_sequence,
348 double *output_sequence)
349 // insert zeros between each element of the input sequence and
350 // convolve with the filter to interpolate the data
353 int endpoint = length + filtlen - 2;
357 // summation with previous convolution
358 // every other dot product interpolates the data
359 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
361 *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
362 *output_sequence++ += dot_product_even(input_sequence + j, filter, filtlen);
365 *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
369 // first convolution of pair
370 // every other dot product interpolates the data
371 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
373 *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
374 *output_sequence++ = dot_product_even(input_sequence + j, filter, filtlen);
377 *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
383 int64_t DenoiseEffect::reconstruct_branches(double *in_low,
386 WaveletFilters *recon_filter,
389 // take input data and filters and form two branches of half the
390 // original length. length of branches is returned
391 convolve_int_2(in_low, in_length, recon_filter->h,
392 recon_filter->length, 0, output);
393 convolve_int_2(in_high, in_length, recon_filter->g,
394 recon_filter->length, 1, output);
395 return in_length * 2;
398 int DenoiseEffect::wavelet_reconstruction(double **in_data,
405 in_length = in_length >> levels;
406 // destination of all but last branch reconstruction is the next
407 // higher intermediate approximation
408 for(i = levels - 1; i > 0; i--)
410 output = in_data[2 * (i - 1)];
411 in_length = reconstruct_branches(in_data[2 * i],
412 in_data[(2 * i) + 1],
418 // destination of the last branch reconstruction is the output data
419 reconstruct_branches(in_data[0],
428 void DenoiseEffect::process_window()
431 for(j = 0; j < iterations; j++)
433 wavelet_decomposition(dsp_in, window_size, ex_coeff_d->values);
435 tree_copy(ex_coeff_r->values, ex_coeff_d->values, window_size, levels);
436 tree_copy(ex_coeff_rn->values, ex_coeff_d->values, window_size, levels);
439 //printf("DenoiseEffect::process_window %f\n", config.level);
440 threshold(window_size, config.level * 10.0, levels);
442 wavelet_reconstruction(ex_coeff_r->values, window_size, dsp_iteration);
443 wavelet_reconstruction(ex_coeff_rn->values, window_size, dsp_in);
445 for(i = 0; i < window_size; i++)
446 dsp_out[i] += dsp_iteration[i];
453 int DenoiseEffect::process_realtime(int64_t size, Samples *input_ptr, Samples *output_ptr)
455 load_configuration();
459 int64_t size_factor = (int)(pow(2, levels));
460 dsp_in = new double[window_size * size_factor];
461 dsp_out = new double[window_size * 2];
462 dsp_iteration = new double[window_size * 2];
465 ex_coeff_d = new Tree(window_size, levels);
466 ex_coeff_r = new Tree(window_size, levels);
467 ex_coeff_rn = new Tree(window_size, levels);
468 wave_coeff_d = new WaveletCoeffs(alpha, beta);
469 wave_coeff_r = new WaveletCoeffs(alpha, beta);
470 decomp_filter = new WaveletFilters(wave_coeff_d, DECOMP);
471 recon_filter = new WaveletFilters(wave_coeff_r, RECON);
472 in_scale = 65535 / sqrt(window_size) / iterations;
473 out_scale = output_level / 65535 * sqrt(window_size);
477 // Append input buffer
478 if(input_size + size > input_allocation)
480 double *new_input = new double[input_size + size];
483 memcpy(new_input, input_buffer, sizeof(double) * input_size);
484 delete [] input_buffer;
486 input_buffer = new_input;
487 input_allocation = input_size + size;
489 memcpy(input_buffer + input_size,
490 input_ptr->get_data(),
491 size * sizeof(double));
495 // Have enough to do some windows
496 while(input_size >= window_size)
499 for(int i = 0; i < window_size; i++)
501 dsp_in[i] = input_buffer[i] * in_scale;
503 bzero(dsp_out, sizeof(double) * window_size);
510 // First window produces garbage
520 // Crossfade into the output buffer
521 int64_t new_allocation = output_size + window_size;
522 if(new_allocation > output_allocation)
524 double *new_output = new double[new_allocation];
528 memcpy(new_output, output_buffer, sizeof(double) * output_size);
529 //printf("CrossfadeFFT::process_fifo 1 %p\n", output_buffer);
530 delete [] output_buffer;
531 //printf("CrossfadeFFT::process_fifo 2\n");
533 output_buffer = new_output;
534 output_allocation = new_allocation;
537 if(output_size >= WINDOW_BORDER)
539 for(int i = 0, j = output_size - WINDOW_BORDER;
543 double src_level = (double)i / WINDOW_BORDER;
544 double dst_level = (double)(WINDOW_BORDER - i) / WINDOW_BORDER;
545 output_buffer[j] = output_buffer[j] * dst_level + out_scale * dsp_out[i] * src_level;
548 for(int i = 0; i < window_size - WINDOW_BORDER; i++)
549 output_buffer[output_size + i] = dsp_out[WINDOW_BORDER + i] * out_scale;
550 output_size += window_size - WINDOW_BORDER;
554 // First buffer has no crossfade
555 memcpy(output_buffer + output_size,
557 sizeof(double) * window_size);
558 output_size += window_size;
562 // Shift input buffer forward
563 for(int i = window_size - WINDOW_BORDER, j = 0;
566 input_buffer[j] = input_buffer[i];
567 input_size -= window_size - WINDOW_BORDER;
571 // Have enough to send to output
572 if(output_size - WINDOW_BORDER >= size)
574 memcpy(output_ptr->get_data(), output_buffer, sizeof(double) * size);
575 for(int i = size, j = 0; i < output_size; i++, j++)
576 output_buffer[j] = output_buffer[i];
581 //printf("DenoiseEffect::process_realtime 1\n");
582 bzero(output_ptr->get_data(), sizeof(double) * size);
594 Tree::Tree(int input_length, int levels)
596 this->input_length = input_length;
597 this->levels = levels;
600 // create decomposition tree
601 values = new double*[2 * levels];
603 for (i = 0; i < levels; i++)
611 values[2 * i] = new double[j + 5];
612 values[2 * i + 1] = new double[j + 5];
620 for (i = 2 * levels - 1; i >= 0; i--)
626 WaveletCoeffs::WaveletCoeffs(double alpha, double beta)
629 double tcosa = cos(alpha);
630 double tcosb = cos(beta);
631 double tsina = sin(alpha);
632 double tsinb = sin(beta);
634 // calculate first two wavelet coefficients a = a(-2) and b = a(-1)
635 values[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
636 + 2.0 * tsinb * tcosa) / 4.0;
637 values[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
638 - 2.0 * tsinb * tcosa) / 4.0;
640 tcosa = cos(alpha - beta);
641 tsina = sin(alpha - beta);
643 // calculate last four wavelet coefficients c = a(0), d = a(1),
644 // e = a(2), and f = a(3)
645 values[2] = (1.0 + tcosa + tsina) / 2.0;
646 values[3] = (1.0 + tcosa - tsina) / 2.0;
647 values[4] = 1 - values[0] - values[2];
648 values[5] = 1 - values[1] - values[3];
650 // zero out very small coefficient values caused by truncation error
651 for (i = 0; i < 6; i++)
653 if (fabs(values[i]) < 1.0e-15) values[i] = 0.0;
657 WaveletCoeffs::~WaveletCoeffs()
662 WaveletFilters::WaveletFilters(WaveletCoeffs *wave_coeffs, wavetype transform)
666 // find the first non-zero wavelet coefficient
668 while(wave_coeffs->values[i] == 0.0) i++;
670 // find the last non-zero wavelet coefficient
672 while(wave_coeffs->values[j] == 0.0) j--;
674 // Form the decomposition filters h~ and g~ or the reconstruction
675 // filters h and g. The division by 2 in the construction
676 // of the decomposition filters is for normalization.
678 for(k = 0; k < length; ++i, --j, ++k)
680 if (transform == DECOMP)
682 h[k] = wave_coeffs->values[j] / 2.0;
683 g[k] = (double) (((i & 0x01) * 2) - 1) * wave_coeffs->values[i] / 2.0;
687 h[k] = wave_coeffs->values[i];
688 g[k] = (double) (((j & 0x01) * 2) - 1) * wave_coeffs->values[j];
692 // clear out the additional array locations, if any
700 WaveletFilters::~WaveletFilters()
712 DenoiseConfig::DenoiseConfig()
717 void DenoiseConfig::copy_from(DenoiseConfig &that)
722 int DenoiseConfig::equivalent(DenoiseConfig &that)
724 return EQUIV(level, that.level);
727 void DenoiseConfig::interpolate(DenoiseConfig &prev,
731 int64_t current_frame)
733 double next_scale = (double)(current_frame - prev_frame) / (next_frame - prev_frame);
734 double prev_scale = (double)(next_frame - current_frame) / (next_frame - prev_frame);
735 this->level = prev.level * prev_scale + next.level * next_scale;
756 DenoiseWindow::DenoiseWindow(DenoiseEffect *plugin)
757 : PluginClientWindow(plugin, xS(200), yS(60), xS(200), yS(60), 0)
759 this->plugin = plugin;
762 void DenoiseWindow::create_objects()
764 int xs10 = xS(10), xs70 = xS(70);
766 int x = xs10, y = ys10;
768 add_subwindow(new BC_Title(x, y, _("Level:")));
770 add_subwindow(scale = new DenoiseLevel(plugin, x, y));
777 void DenoiseWindow::update()
779 scale->update(plugin->config.level);
793 DenoiseLevel::DenoiseLevel(DenoiseEffect *plugin, int x, int y)
794 : BC_FPot(x, y, (float)plugin->config.level, 0, 1.0)
796 this->plugin = plugin;
800 int DenoiseLevel::handle_event()
802 plugin->config.level = get_value();
803 plugin->send_configure_change();