HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CL_FreqFilter.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: CL_FreqFilter.h ( Clip Library, C++)
7  *
8  * COMMENTS:
9  * Does continuous digital filtering in the frequency domain
10  */
11 
12 
13 #ifndef __CL_FreqFilter__
14 #define __CL_FreqFilter__
15 
16 #include "CL_API.h"
17 #include <UT/UT_Array.h>
18 #include <UT/UT_Complex.h>
19 #include <UT/UT_FFT.h>
20 #include <UT/UT_Interrupt.h>
21 #include <UT/UT_UniquePtr.h>
22 #include <SYS/SYS_Types.h>
23 
24 #define SQRT_TWO 1.4142136
25 
26 class CL_API CL_FreqFilter final
27 {
28 public:
29 
30  // Filter Type, used for Pass filters
31  enum class Type
32  {
33  LOW_PASS,
34  HIGH_PASS,
35  BAND_PASS,
36  BAND_STOP
37  };
38 
39  // Filter Design, used for Pass filters
40  enum class Design
41  {
42  STANDARD,
43  BUTTERWORTH
44  };
45 
46  // Filter shape, used for Equalization Filters
47  enum class Shape
48  {
49  GAUSSIAN,
50  BOX
51  };
52 
53  CL_FreqFilter(int size,
54  fpreal overlap = 0.1,
55  fpreal discard = 0.1,
56  bool blend_chunks = true,
57  bool filter_phase = false);
58 
59  ~CL_FreqFilter();
60 
61  void resize(int size);
62 
63  void setFilter(int size, fpreal *filter, fpreal *phase = nullptr);
64 
65  /// Ensures the filter has been set and the window size is valid
66  bool isReady() const {
67  if (myFilterPhase && myPhaseFilter == nullptr)
68  return false;
69 
70  return myFilter != nullptr && myWindowSize > 0;
71  }
72 
73  void setOverlap(fpreal overlap) { myOverlap = overlap; }
74  void setDiscard(fpreal discard) { myDiscard = discard; }
75  void setBlendChunks(bool blend) { myBlendChunks = blend; }
76  void setFilterPhase(bool filter_phase) { myFilterPhase = filter_phase; }
77 
78  /// Designs a pass filter, filling out the filter and phase_filter arrays
79  /// Filter MUST be size floats long.
80  template <class F>
81  static void designPassFilter(Type type, Design design,
82  fpreal freq1, fpreal freq2,
83  fpreal attenuation, fpreal pass_gain, int order,
84  int size, fpreal freqstep, int index, int num,
85  const F& override_parameters,
86  fpreal *filter, fpreal *phase_filter = nullptr);
87 
88  /// Variant of the designPassFilter method which does not have a parameter
89  /// override method
90  static void designPassFilter(Type type, Design design,
91  fpreal freq1, fpreal freq2,
92  fpreal attenuation, fpreal pass_gain, int order,
93  int size, fpreal freqstep,
94  fpreal *mag_filter, fpreal *phase_filter = nullptr)
95  {
96  // Override Parameters does nothing
97  auto override_parameters
98  = [](fpreal &, fpreal &, fpreal &, fpreal &, int, int) {};
99 
101  type, design, freq1, freq2, attenuation, pass_gain, order,
102  size, freqstep, 0, 0, override_parameters,
103  mag_filter, phase_filter);
104  }
105 
106  // Designs a eq filter, filling out the filter array
107  // Filter MUST be size floats long.
108  template <class F>
109  static void designBandEQFilter(Shape shape, fpreal basefreq, fpreal *bands, int num_bands,
110  int size, fpreal freqstep, int index, int num,
111  const F& override_parameters,
112  fpreal *filter);
113 
114  // Designs a eq filter, filling out the filter array
115  // Filter MUST be size floats long.
116  template <class F>
117  static void designParaEQFilter(fpreal shape, fpreal freq, fpreal fgain, fpreal pgain,
118  fpreal width, fpreal attenuation,
119  int size, fpreal freqstep, int index, int num,
120  const F& override_parameters,
121  fpreal *filter);
122 
123  // Designs an average sideband filter, filling out the filter and phase_filter arrays
124  // Filter MUST be size floats long.
125  template <class F>
126  static void designAverageSideband(fpreal fgain, fpreal gain, fpreal effect, int length, int passfreq,
127  int size, const F& evaluate_data,
128  fpreal *filter);
129 
130  // Designs an average sideband filter, filling out the filter and phase_filter arrays
131  // Filter MUST be size floats long.
132  template <class F>
133  static void designContinuousSideband(fpreal fgain, fpreal gain, fpreal effect, int passfreq,
134  int size, int index, const F& evaluate_data,
135  fpreal *filter);
136 
137  /// Filter input data by applying the frequency domain filter
138  /// @param evaluate_data lambda function which returns the data value at the
139  /// given index
140  /// @param design_filter lambda function which redesigns the filter at the
141  /// currently processed chunk.
142  template <class F, class FD>
143  void filter(fpreal* dest, int dest_capacity, fpreal freqstep,
144  const F& evaluate_data, const FD& design_filter);
145 
146  /// Variant of the filter method which does not need to redesign the filter depending on the chunk
147  template <class F>
148  void filter(
149  fpreal *dest,
150  int dest_capacity,
151  fpreal freqstep,
152  const F &evaluate_data)
153  {
154  // Don't redesign filter depending on the chunk
155  auto design_filter = [](fpreal, int, int, fpreal *, fpreal *) {};
156  filter(dest, dest_capacity, freqstep, evaluate_data, design_filter);
157  }
158 
159  /// Filter method which constrains the endpoints by "bending" the smoothed
160  /// curve to match the start and end point of the original curve.
161  ///
162  /// @param constrain_region - Number of samples over which to bend the curve
163  /// at the beginning and end of the curve. Setting this equal to the sample
164  /// rate, attenuated by the high frequency cutoff tends to give good
165  /// results for smoothing animation curves.
166  /// Example: constrain_region = SYSmax((1.0 - high_cutoff) * sample_rate), 1)
167  /// Setting this to a value <= 0 will skip the constraining step.
168  ///
169  /// @param shift - If true, will apply a y-offset to the input data so the
170  /// starting value is zero before smoothing. The smoothed curve will then be
171  /// shifted back after smoothing. This helps avoid the step response which
172  /// occurs when smoothing data with a non-zero starting y-value. This can be
173  /// especially seen when applying the filter to a constant function
174  template <class F>
175  void constrainedFilter(
176  fpreal *dest,
177  int dest_capacity,
178  fpreal freqstep,
179  const F &evaluate_data,
180  int constrain_region,
181  bool shift = true);
182 
183 private:
184  void window(fpreal *data1, fpreal *data2, fpreal *windower);
185  void computeWindow();
186 
187  static fpreal attenuate(fpreal value, fpreal power);
188 
189  int myWindowSize;
190 
191  fpreal *myFilter;
192  fpreal *myPhaseFilter;
193 
194  fpreal myLastDataSample;
195  fpreal myLastSlope;
196 
197  // Overlap is the region (specified as a fraction of the window size) of
198  // which neightboring chunks should overlap
199  fpreal myOverlap;
200  // Discard is the region (specified as a fraction of the window size) of
201  // data at the beginning and ending of each chunk which should be ignored
202  fpreal myDiscard;
203  // Whether or not the smoothed curves for neighboring chunks should be
204  // blended in the overlap region
205  bool myBlendChunks;
206  // Whether or not the phase filter should be applied, or just the magnitude
207  // filter. (phase filter must be non-null)
208  bool myFilterPhase;
209 
210  UT_Array<fpreal> myWindow;
211  UT_Array<fpreal> myUnWindow;
212 
213  UT_Array<fpreal> myReal;
214  UT_Array<fpreal> myImag;
215 
216  UT_FFT<fpreal> myTransform;
217 };
218 
219 template <class F, class FD>
220 void
222  fpreal *dest,
223  int dest_capacity,
224  fpreal freqstep,
225  const F &evaluate_data,
226  const FD &design_filter)
227 {
228  int j, index, n;
229  int garbage = int(myWindowSize*(myDiscard / 2.0));
230  int datachunk = int(myWindowSize - garbage*2);
231  int blendreg = int(datachunk * myOverlap);
232  int chunk;
233  UT_Interrupt *boss = UTgetInterrupt();
234  fpreal cph, sph;
235  fpreal v1, v2;
236  short int intcnt = 0;
237 
238  fpreal *data1 = new fpreal[myWindowSize];
239  fpreal *data2 = new fpreal[myWindowSize];
240  fpreal *readdata;
241  fpreal *blend = nullptr;
242 
243  if(blendreg >= datachunk)
244  blendreg = (int)SYSfloor(datachunk*0.95);
245 
246  readdata = new fpreal[blendreg];
247 
248  if(datachunk-blendreg < dest_capacity)
249  {
250  blend = new fpreal[blendreg];
251  cph = M_PI / fpreal(blendreg+1);
252  sph = -0.5 * (fpreal)M_PI;
253  for(j=0; j<blendreg; j++)
254  blend[j] = 0.5 + 0.5*SYSsin(sph + cph * fpreal(j+1));
255  }
256 
257  if (boss->opStart("Filtering Data"))
258  {
259  for(index=0, chunk=0; index<dest_capacity;
260  index+=(datachunk-blendreg)*2, chunk+=2)
261  {
262  if(!intcnt-- && boss->opInterrupt(100*index/dest_capacity))
263  break;
264 
265  // get data
266  if(myBlendChunks)
267  {
268  evaluate_data(fpreal(index), fpreal(index+datachunk-1),
269  data1+garbage, datachunk);
270 
271  evaluate_data(fpreal(index+datachunk-blendreg),
272  fpreal(index-blendreg+datachunk*2-1),
273  data2+garbage, datachunk);
274  }
275  else
276  {
277  evaluate_data(fpreal(index), fpreal(index+datachunk-1-blendreg),
278  data1+garbage, datachunk-blendreg);
279  evaluate_data(fpreal(index+datachunk-blendreg),
280  fpreal(index-blendreg*2+datachunk*2-1),
281  data2+garbage, datachunk-blendreg);
282 
283  for(j=blendreg; j--;)
284  {
285  data1[j+datachunk-blendreg] = 0;
286  data2[j+datachunk-blendreg] = 0;
287  }
288  }
289 
290  if(index>0)
291  {
292  for(j=blendreg; j--;)
293  data1[j+garbage] = readdata[j];
294  }
295 
296  for(j=blendreg; j--;)
297  readdata[j] = data2[garbage+datachunk-blendreg + j];
298 
299  for(j=garbage; j--;)
300  {
301  data1[j] =0;
302  data2[j] =0;
303  data1[myWindowSize-garbage+j] =0;
304  data2[myWindowSize-garbage+j] =0;
305  }
306 
307  window(data1, data2, myWindow.data());
308 
309  // do freq transform; multiply; do time transform
310  myTransform.toFrequency(data1, data2, myReal.data(), myImag.data(), myWindowSize);
311 
312  // Update the filter depending on the current index & chunk
313  design_filter(freqstep, index, chunk, myFilter, myPhaseFilter);
314 
315  if(!myFilterPhase)
316  {
317  // Magnitude scale only
318  for(j=0; j<myWindowSize; j++)
319  {
320  myReal[j] *= myFilter[j];
321  myImag[j] *= myFilter[j];
322  }
323  }
324  else
325  {
326  // Mag & phase filter
327  for(j=myWindowSize; j--;)
328  {
329  myReal[j] *= myFilter[j];
330  myImag[j] *= myPhaseFilter[j];
331  }
332  }
333 
334  myTransform.toTime(myReal.data(), myImag.data(), data1, data2, myWindowSize);
335 
336  window(data1, data2, myUnWindow.data());
337 
338  if(index>0)
339  {
340  // blend data part 1
341  n = blendreg;
342  if(index+blendreg >= dest_capacity)
343  {
344  if(index >= dest_capacity)
345  n = 0;
346  else
347  n = dest_capacity - index;
348  }
349 
350  if(myBlendChunks)
351  {
352  for(j=n; j--;)
353  {
354  v1 = dest[index+j]*(1-blend[j]);
355  v2 = data1[j+garbage] * blend[j];
356 
357  dest[index+j] = v2 + v1;
358  }
359  }
360  else
361  {
362  for(j=n; j--;)
363  dest[index+j] += data1[j+garbage];
364  }
365 
366  // store data part1 & 2
367  for(j=blendreg; j<datachunk && index+j<dest_capacity; j++)
368  dest[index+j] = data1[j+garbage];
369 
370  for(j=blendreg; j<datachunk && index+datachunk-blendreg+j<dest_capacity; j++)
371  dest[index+datachunk-blendreg+j] = data2[j+garbage];
372 
373  // blend data part 2
374  if(index+datachunk-blendreg < dest_capacity)
375  {
376  n = blendreg;
377  if(index+datachunk >= dest_capacity)
378  n = dest_capacity - (index+datachunk-blendreg);
379 
380  if(myBlendChunks)
381  {
382  for(j=n; j--;)
383  {
384  v1 = dest[index+datachunk-blendreg+j]*(1-blend[j]);
385  v2 = data2[j+garbage] * blend[j];
386 
387  dest[index+datachunk-blendreg+j] = v2 + v1;
388  }
389  }
390  else
391  {
392  for(j=n; j--;)
393  dest[index+datachunk-blendreg+j]+=data2[j+garbage];
394  }
395  }
396  }
397  else
398  {
399  // store data parts 1
400  n = datachunk;
401  if(index+datachunk >= dest_capacity)
402  {
403  if(index >= dest_capacity)
404  n = 0;
405  else
406  n = dest_capacity - index;
407  }
408 
409  for(j=n; j--;)
410  dest[index+j] = data1[j+garbage];
411 
412  // blend data part 2
413  if(index+datachunk-blendreg < dest_capacity)
414  {
415  n = blendreg;
416  if(index+datachunk >= dest_capacity)
417  n = dest_capacity - (index+datachunk-blendreg);
418 
419  if(myBlendChunks)
420  {
421  for(j=n; j--;)
422  {
423  v1 = dest[index+datachunk-blendreg+j]*(1-blend[j]);
424  v2 = data2[j+garbage] * blend[j];
425 
426  dest[index-blendreg+datachunk+j] = v2 + v1;
427  }
428  }
429  else
430  {
431  for(j=n; j--;)
432  dest[index-blendreg+datachunk+j] += data2[j+garbage];
433  }
434  }
435 
436  // store data part 2
437  for(j=blendreg; j<datachunk && index+datachunk-blendreg+j < dest_capacity; j++)
438  dest[index+datachunk-blendreg+j] = data2[j+garbage];
439 
440  }
441  }
442  }
443  boss->opEnd();
444 
445  if (blend)
446  delete [] blend;
447 
448  delete [] data1;
449  delete [] data2;
450  delete [] readdata;
451 }
452 
453 template <class F>
454 void
456  fpreal *dest,
457  int dest_capacity,
458  fpreal freqstep,
459  const F &evaluate_data,
460  int constrain_region,
461  bool shift)
462 {
463  auto shifted = UTmakeUnique<fpreal[]>(dest_capacity);
464  evaluate_data(0, dest_capacity - 1, shifted.get(), dest_capacity);
465 
466  fpreal original_start = shifted[0];
467 
468  // shift so that start = 0
469  // this avoids the initial step response, especially when applied to constant functions
470  if (shift)
471  {
472  for (exint i = 0; i < dest_capacity; i++)
473  shifted[i] -= original_start;
474  }
475 
476  fpreal start = shifted[0];
477  fpreal end = shifted[dest_capacity - 1];
478 
479  auto evaluate_shifted = [&](int start, int, fpreal *data, int size)
480  {
481  for (exint i = 0; i < size; i++)
482  data[i] = (start + i) < dest_capacity ? shifted[start + i] : shifted[dest_capacity - 1];
483  };
484 
485  filter(dest, dest_capacity, freqstep, evaluate_shifted);
486 
487  if (constrain_region > dest_capacity / 2)
488  constrain_region = dest_capacity / 2;
489 
490  fpreal delta_start = start - dest[0];
491 
492  fpreal delta_end = end - dest[dest_capacity - 1];
493 
494  for (exint i = 0; i < dest_capacity / 2; i++)
495  {
496  fpreal t = static_cast<fpreal>(i) / constrain_region;
497  fpreal a = SYSpow(2.0, -10.0 * t);
498 
499  dest[i] += delta_start * a;
500 
501  dest[dest_capacity - 1 - i] += delta_end * a;
502  }
503 
504  // shift back to original range
505  if (shift)
506  {
507  for (exint i = 0; i < dest_capacity; i++)
508  dest[i] += original_start;
509  }
510 }
511 
512 // Designs a filter, populating the filter and phase_filter arrays
513 template <class F>
514 void
516  fpreal freq1, fpreal freq2,
517  fpreal attenuation, fpreal pass_gain, int order,
518  int chunk_size, fpreal freqstep, int index, int num,
519  const F& override_parameters,
520  fpreal *filter, fpreal *phase_filter)
521 {
522  fpreal f,k;
523  UT_Complex s,w1,w2,w3, w1_inv, response;
524  fpreal dc_gain,gain, exp1, filter_order;
525  int i;
526  int half = chunk_size/2;
527  fpreal max;
528 
529  exp1 = 2.0;
530  filter_order = static_cast<fpreal>(order); // Order of the filter. Used for Butterworth filters
531 
532  // Potentially overrides the frequency, gain, and myAttenuationuation parameters
533  // depending on the index and num
534  override_parameters(freq1, freq2, gain, attenuation, index, num);
535 
536  if(freq1<=SYS_FTOLERANCE)
537  freq1 = SYS_FTOLERANCE;
538  if(freq2<=SYS_FTOLERANCE)
539  freq2 = SYS_FTOLERANCE;
540 
541  switch(type)
542  {
543  case Type::LOW_PASS:
544  w2.set(1,0);
545  w1.set(freq2,0); // Cutoff frequency (high)
546  w1_inv = w2 / w1;
547 
548  s.set(0,freq2/100);
549 
550  k = SYSpow(freq2,exp1);
551 
552  if (design == Design::BUTTERWORTH)
553  // Frequency response of the Butterworth Filter
554  // See: https://en.wikipedia.org/wiki/Butterworth_filter#Transfer_function
555  response = (w2) / (w2 + (-(s * w1_inv).pow(exp1)).pow(filter_order));
556  else
557  response = (w2*k)/((s+w1).pow(exp1)); // default type
558 
559  dc_gain = response.magnitude();
560  if(SYSequalZero(dc_gain))
561  dc_gain=1.0;
562  gain = 1.0/dc_gain;
563 
564  for(i=0,f=0; i<=half; i++,f+=freqstep)
565  {
566  s.set(0,f);
567 
568  if (design == Design::BUTTERWORTH)
569  {
570  if (f < SYS_FTOLERANCE)
571  s.set(0, SYS_FTOLERANCE);
572  response = (w2) / (w2 + (-(s * w1_inv).pow(exp1)).pow(filter_order));
573  }
574  else
575  response = (w2*k)/((s+w1).pow(exp1));
576 
577  if(phase_filter)
578  {
579  filter[i] = response.real();
580  phase_filter[i] = response.imaginary();
581  }
582  else
583  filter[i] = response.magnitude()*gain;
584  }
585  break;
586 
587  case Type::HIGH_PASS:
588  w1.set(freq1,0);
589  w2.set(1, 0);
590  s.set(0,freq1*100);
591 
592  k = SYSpow(freq1,exp1);
593  if (design == Design::BUTTERWORTH)
594  response = (w2) / (w2 + (-(w1 / s).pow(exp1)).pow(filter_order));
595  else
596  response = (s.pow(exp1)*k) / (s+w1).pow(exp1);
597 
598  dc_gain = response.magnitude();
599  if(SYSequalZero(dc_gain))
600  dc_gain=1.0;
601 
602  gain = 1.0/dc_gain;
603 
604  for(i=0,f=0; i<=half; i++,f+=freqstep)
605  {
606  s.set(0,f);
607 
608  if (design == Design::BUTTERWORTH)
609  {
610  if (f < SYS_FTOLERANCE)
611  s.set(0, SYS_FTOLERANCE);
612 
613  response = (w2) / (w2 + (-(w1 / s).pow(exp1)).pow(filter_order));
614  }
615  else
616  response = (s.pow(exp1)*k)/(s+w1).pow(exp1);
617 
618  if(phase_filter)
619  {
620  filter[i] = response.real();
621  phase_filter[i] = response.imaginary();
622  }
623  else
624  filter[i] = response.magnitude()*gain;
625 
626  }
627  break;
628 
629  case Type::BAND_STOP:
630  case Type::BAND_PASS:
631 
632  if(freq1>freq2)
633  {
634  f=freq1;
635  freq1=freq2;
636  freq2=f;
637  }
638 
639  w1.set(freq1,0);
640  w2.set(freq2,0);
641  s.set(0,freq1);
642 
643  k = freq1 * freq2;
644  response = (s*k) /((s+w1) * (s+w2));
645 
646  dc_gain = response.magnitude();
647  if(SYSequalZero(dc_gain))
648  dc_gain=1.0;
649 
650  gain = 1.0/dc_gain;
651 
652  for(i=0,f=0; i<=half; i++,f+=freqstep)
653  {
654  s.set(0,f);
655 
656  response = (s*k)/((s+w1)*(s+w2));
657 
658  if(phase_filter)
659  {
660  filter[i] = response.real();
661  phase_filter[i] = response.imaginary();
662  }
663  else
664  filter[i] = response.magnitude()*gain;
665  }
666  break;
667  }
668 
669  if(!SYSisEqual(attenuation,1.0))
670  {
671  for(i=0; i<=half; i++)
672  filter[i] = SYSpow(filter[i],attenuation);
673  }
674 
675  gain = SYSpow(CONST_FPREAL(10), pass_gain/20);
676  max = 0;
677  for(i=0; i<=half; i++)
678  if(filter[i] > max)
679  max = filter[i];
680 
681  max = (max!=0)?(gain/max):gain;
682 
683  for(i=0; i<=half; i++)
684  filter[i] *=max;
685 
686  // band stop is just 1-bandpass
687  if(type == Type::BAND_STOP)
688  {
689  if(phase_filter)
690  {
691  for(i=0; i<=half; i++)
692  {
693  max = SYSsqrt(filter[i]*filter[i] + phase_filter[i]*phase_filter[i]);
694  max = (gain-max)/max;
695  filter[i]*= max;
696  phase_filter[i]*=max;
697  }
698  }
699  else
700  {
701  for(i=0; i<=half; i++)
702  filter[i] = gain - filter[i];
703  }
704  }
705 
706 
707  for(i=1; i<half; i++)
708  {
709  filter[chunk_size-i] = filter[i];
710  if(phase_filter)
711  phase_filter[chunk_size-i] = -phase_filter[i];
712  }
713 }
714 
715 template <class F>
716 void
717 CL_FreqFilter::designBandEQFilter(Shape shape, fpreal basefreq, fpreal *bands, int num_bands,
718  int chunk_size, fpreal freqstep, int index, int num,
719  const F& override_parameters,
720  fpreal *filter)
721 {
722  int i, j;
723  fpreal band;
724  int starti,endi;
725  fpreal freq,amp;
726  fpreal lstart, lend, lfreq, lwidth, lcenter;
727 
728  override_parameters(basefreq, bands, index, num);
729 
730  band = CONST_FPREAL(10) * SYSpow(CONST_FPREAL(2), basefreq);
731 
733 
734  for(i=0; i< chunk_size; i++)
735  filter[i] = 1.0;
736 
737  endi = 0;
738  for(i=0; i<num_bands; i++)
739  {
740  if(endi>=chunk_size/2)
741  continue;
742 
743  if(shape == Shape::BOX)
744  starti = endi;
745  else
746  starti = int((band/width)/freqstep + 0.5);
747  endi = int((band*width)/freqstep + 1.5);
748 
749  if(bands[i] != 1.0)
750  {
751  lstart = SYSlog(band/width);
752  lend = SYSlog(band*width);
753  lwidth = (lend-lstart)/2.0;
754  lcenter = SYSlog(band);
755 
756  // filter contraints.
757  if(starti>chunk_size/2)
758  continue;
759  if(endi>chunk_size/2)
760  endi = chunk_size/2;
761  if(starti<0)
762  starti=0;
763 
764 
765  for(j=starti; j<=endi; j++)
766  {
767  lfreq = SYSlog(j*freqstep);
768  if(shape == Shape::GAUSSIAN)
769  {
770  freq = 3.0*(lfreq-lcenter)/lwidth;
771 
772  amp = SYSexp( -0.5 *freq*freq);
773  amp = (bands[i] - 1.0) * amp;
774 
775  filter[j] +=amp;
776  }
777  else
778  {
779  filter[j] = bands[i];
780  }
781  }
782  }
783 
784  band *=2;
785  }
786 
787  for(i=1; i<chunk_size/2; i++)
788  filter[chunk_size-i] = filter[i];
789 
790 }
791 
792 template <class F>
793 void
795  fpreal width, fpreal atten,
796  int chunk_size, fpreal freqstep, int index, int num,
797  const F& override_parameters,
798  fpreal *filter)
799 {
800  int i, j;
801  int starti,endi;
802  fpreal amp,ind;
803  fpreal lstart, lend, lfreq, lwidth, lcenter;
804 
805  override_parameters(freq, width, fgain, pgain, shape, atten, index, num);
806 
807  for(i=0; i<chunk_size; i++)
808  filter[i] = pgain;
809 
810  width = SYSpow(2.0,width/2.0);
811 
812  starti = int((freq/width)/freqstep + 0.5);
813  endi = int((freq*width)/freqstep + 1.5);
814 
815  // filter contraints.
816  if(starti>chunk_size/2)
817  return;
818  if(endi>chunk_size/2)
819  endi = chunk_size/2;
820  if(starti<0)
821  starti=0;
822 
823  lstart = SYSlog(freq/width);
824  lend = SYSlog(freq*width);
825  lcenter = SYSlog(freq);
826  lwidth = (lend - lstart)/2.0;
827 
828  ind = freq/width;
829  for(j=starti; j<=endi; j++, ind+=freqstep)
830  {
831  lfreq = SYSlog(ind);
832 
833  if(shape <= 0.5)
834  {
835  // blend between box (0) and triangle (0.5)
836  amp = (1.0-SYSabs(lfreq-lcenter)/lwidth);
837 
838  amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
839  amp = fgain*(0.5-shape)*2 + amp*(shape*2);
840 
841  filter[j] =amp;
842  }
843  else
844  {
845  // blend between triangle (0.5) and gaussian (1.0)
846  freq = 3.0*(lfreq-lcenter)/lwidth;
847 
848  amp = SYSexp( -0.5 *freq*freq);
849  amp = (fgain - pgain)*attenuate(amp,atten) + pgain;
850  filter[j] = amp*(shape-0.5)*2;
851 
852  amp = 1.0-SYSabs(lfreq-lcenter)/lwidth;
853  amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
854  filter[j] += amp*(1.0-shape)*2;
855  }
856  }
857 
858  for(i=1; i<chunk_size/2; i++)
859  filter[chunk_size-i] = filter[i];
860 
861 }
862 
863 template <class F>
864 void
865 CL_FreqFilter::designAverageSideband(fpreal fgain, fpreal gain, fpreal effect, int length, int passfreq,
866  int chunk_size, const F& evaluate_data,
867  fpreal *filter)
868 {
869  UT_FFT<fpreal> freqtransform(chunk_size);
870 
871  fpreal *temp;
872  int i,j,count;
873  fpreal *real,*imag;
874  fpreal factor;
875  int start,end;
876  double max;
877 
878  temp = new fpreal[chunk_size];
879  real = new fpreal[chunk_size];
880  imag = new fpreal[chunk_size];
881 
882  for(j=1; j<chunk_size; j++)
883  filter[j]=0;
884 
885 
886  for(i=0,count=0; i<length; i+=chunk_size,count++)
887  {
888  evaluate_data((fpreal)i, (fpreal)i+chunk_size-1,temp,chunk_size);
889 
890  // do freq transform; multiply; do time transform
891  freqtransform.toFrequency(temp,real,imag,chunk_size);
892 
893  // calculate power spectrum
894  for(j=1; j<chunk_size; j++)
895  temp[j] = (real[j]*real[j] + imag[j]*imag[j]);
896 
897  for(j=1; j<chunk_size; j++)
898  filter[j] += temp[j];
899  }
900 
901  if(count>1)
902  for(j=1; j<chunk_size; j++)
903  filter[j] /= count;
904 
905  // 'effect' widens the sideband filter on a sample basis
906  if(!SYSisEqual(effect,1.0))
907  {
908  for(j=(chunk_size>>1); j--;)
909  {
910  temp[j] = filter[j];
911  filter[j] = 0;
912  }
913 
914  for(j=(chunk_size>>1); j--;)
915  {
916  start = int(j/effect);
917  end = int(j*effect);
918 
919  if(start < 0)
920  start=0;
921  if(end > chunk_size/2)
922  end = chunk_size/2;
923 
924  for(i=start; i<=end; i++)
925  {
926  // geometric interp
927  factor = SYSlog(fpreal(i)/fpreal(j))/SYSlog(effect);
928 
929  // gaussian distribution
930  factor = SYSexp(-1.5*factor*factor);
931  filter[i] += temp[j] * factor;
932  }
933  }
934 
935  for(j=1; j<(chunk_size>>1); j++)
936  filter[chunk_size-j] = filter[j];
937  }
938  // don't filter DC
939 // filter[0] = 1.0;
940 
941  for(j=0; j<chunk_size; j++)
942  if(j==0 || filter[j]>max)
943  max = filter[j];
944 
945  if(!SYSequalZero(max))
946  for(j=0; j<chunk_size; j++)
947  filter[j] /=max;
948 
949  if(!SYSisEqual(fgain,1.0))
950  for(j=0; j<chunk_size; j++)
951  filter[j] *= fgain;
952 
953  if(passfreq)
954  {
955  for(j=0; j<chunk_size; j++)
956  filter[j] = (1.0+gain) - filter[j];
957  }
958  else
959  {
960  for(j=0; j<chunk_size; j++)
961  filter[j] += gain;
962  }
963 
964  delete [] temp;
965  delete [] real;
966  delete [] imag;
967 }
968 
969 template <class F>
970 void
972  int chunk_size, int index,
973  const F& evaluate_data,
974  fpreal *filter)
975 {
976  UT_FFT<fpreal> freqtransform(chunk_size);
977 
978  int i,j;
979  int start,end;
980  fpreal factor;
981  double max;
982 
983  fpreal *real = new fpreal[chunk_size];
984  fpreal *imag = new fpreal[chunk_size];
985 
986  evaluate_data((fpreal)index, (fpreal)index+chunk_size-1, filter, chunk_size);
987 
988  // do freq transform; multiply; do time transform
989  freqtransform.toFrequency(filter,real,imag,chunk_size);
990 
991  // calculate power spectrum
992  for(j=1; j<chunk_size; j++)
993  filter[j] = (real[j]*real[j] + imag[j]*imag[j]);
994 
995  // 'effect' widens the sideband filter on a sample basis
996  if(!SYSisEqual(effect,1.0))
997  {
998  for(j=0; j<=chunk_size/2; j++)
999  {
1000  real[j] = filter[j];
1001  filter[j] = 0;
1002  }
1003 
1004  for(j=0; j<chunk_size/2; j++)
1005  {
1006  start = int(j/effect);
1007  end = int(j*effect);
1008 
1009  if(start < 0)
1010  start=0;
1011  if(end > chunk_size/2)
1012  end = chunk_size/2;
1013 
1014  for(i=start; i<=end; i++)
1015  {
1016  // geometric interp
1017  factor = SYSlog(fpreal(i)/fpreal(j))/SYSlog(effect);
1018 
1019  // gaussian distribution
1020  factor = SYSexp(-1.5*factor*factor);
1021  filter[i] += real[j] * factor;
1022  }
1023  }
1024 
1025  for(j=1; j<chunk_size/2; j++)
1026  filter[chunk_size-j] = filter[j];
1027  }
1028 
1029  // don't filter DC
1030 // filter[0]=1;
1031 
1032  for(j=0; j<chunk_size; j++)
1033  if(j==0 || filter[j]>max)
1034  max = filter[j];
1035 
1036  if(!SYSequalZero(max))
1037  for(j=0; j<chunk_size; j++)
1038  filter[j] /=max;
1039 
1040  if(!SYSisEqual(fgain,1.0))
1041  for(j=0; j<chunk_size; j++)
1042  filter[j] *= fgain;
1043 
1044  if(passfreq)
1045  {
1046  for(j=0; j<chunk_size; j++)
1047  filter[j] = (gain+1.0) - filter[j];
1048  }
1049  else
1050  {
1051  for(j=0; j<chunk_size; j++)
1052  filter[j] += gain;
1053  }
1054 
1055  delete [] real;
1056  delete [] imag;
1057 }
1058 
1059 #endif
bool isReady() const
Ensures the filter has been set and the window size is valid.
Definition: CL_FreqFilter.h:66
void toTime(T *sourcer, T *sourcei, T *dest, exint npnts)
T real() const
Definition: UT_Complex.h:68
void constrainedFilter(fpreal *dest, int dest_capacity, fpreal freqstep, const F &evaluate_data, int constrain_region, bool shift=true)
GLuint GLdouble GLdouble GLint GLint order
Definition: glew.h:3460
imath_half_bits_t half
if we're in a C-only context, alias the half bits type to half
Definition: half.h:266
#define SQRT_TWO
Definition: CL_FreqFilter.h:24
#define M_PI
Definition: fmath.h:90
GLuint start
Definition: glcorearb.h:475
static void designAverageSideband(fpreal fgain, fpreal gain, fpreal effect, int length, int passfreq, int size, const F &evaluate_data, fpreal *filter)
void filter(fpreal *dest, int dest_capacity, fpreal freqstep, const F &evaluate_data, const FD &design_filter)
int64 exint
Definition: SYS_Types.h:125
#define SYSabs(a)
Definition: SYS_Math.h:1515
void setBlendChunks(bool blend)
Definition: CL_FreqFilter.h:75
GLdouble GLdouble t
Definition: glew.h:1403
UT_ComplexT pow(T exp) const
Definition: UT_Complex.h:195
static void designContinuousSideband(fpreal fgain, fpreal gain, fpreal effect, int passfreq, int size, int index, const F &evaluate_data, fpreal *filter)
static void designParaEQFilter(fpreal shape, fpreal freq, fpreal fgain, fpreal pgain, fpreal width, fpreal attenuation, int size, fpreal freqstep, int index, int num, const F &override_parameters, fpreal *filter)
GLsizeiptr size
Definition: glcorearb.h:664
ImageBuf OIIO_API pow(const ImageBuf &A, cspan< float > B, ROI roi={}, int nthreads=0)
#define CL_API
Definition: CL_API.h:10
vint4 blend(const vint4 &a, const vint4 &b, const vbool4 &mask)
Definition: simd.h:4784
int opInterrupt(int percent=-1)
T magnitude() const
Definition: UT_Complex.h:241
static void designPassFilter(Type type, Design design, fpreal freq1, fpreal freq2, fpreal attenuation, fpreal pass_gain, int order, int size, fpreal freqstep, int index, int num, const F &override_parameters, fpreal *filter, fpreal *phase_filter=nullptr)
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:818
void setOverlap(fpreal overlap)
Definition: CL_FreqFilter.h:73
static void designBandEQFilter(Shape shape, fpreal basefreq, fpreal *bands, int num_bands, int size, fpreal freqstep, int index, int num, const F &override_parameters, fpreal *filter)
#define CONST_FPREAL(c)
Definition: SYS_Types.h:327
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
GLuint GLuint end
Definition: glcorearb.h:475
GLuint num
Definition: glew.h:2695
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
GLint GLint GLint GLint GLint GLint GLint GLbitfield GLenum filter
Definition: glcorearb.h:1297
void setDiscard(fpreal discard)
Definition: CL_FreqFilter.h:74
GLint GLsizei count
Definition: glcorearb.h:405
SYS_API fpreal32 SYSfloor(fpreal32 val)
GLint GLsizei width
Definition: glcorearb.h:103
void toFrequency(T *source, T *destr, T *desti, exint npnts)
GLdouble n
Definition: glcorearb.h:2008
T imaginary() const
Definition: UT_Complex.h:70
bool SYSequalZero(const UT_Vector3T< T > &v)
Definition: UT_Vector3.h:1052
GLboolean * data
Definition: glcorearb.h:131
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
ImageBuf OIIO_API resize(const ImageBuf &src, string_view filtername="", float filterwidth=0.0f, ROI roi={}, int nthreads=0)
fpreal64 fpreal
Definition: SYS_Types.h:277
T * data()
Definition: UT_Array.h:785
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
GLuint index
Definition: glcorearb.h:786
void opEnd(int opid=-1)
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GLfloat f
Definition: glcorearb.h:1926
void setFilterPhase(bool filter_phase)
Definition: CL_FreqFilter.h:76
int opStart(const char *opname=0, int npoll=0, int immediate_dialog=0, int *opid=0)
#define SYS_FTOLERANCE
Definition: SYS_Types.h:208
void set(T r, T i)
Definition: UT_Complex.h:112
Definition: core.h:1131
GLfloat GLfloat v1
Definition: glcorearb.h:817
type
Definition: core.h:1059
static void designPassFilter(Type type, Design design, fpreal freq1, fpreal freq2, fpreal attenuation, fpreal pass_gain, int order, int size, fpreal freqstep, fpreal *mag_filter, fpreal *phase_filter=nullptr)
Definition: CL_FreqFilter.h:90
GLdouble s
Definition: glew.h:1395
void filter(fpreal *dest, int dest_capacity, fpreal freqstep, const F &evaluate_data)
Variant of the filter method which does not need to redesign the filter depending on the chunk...