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(
185  UT_UniquePtr<fpreal[]> &data1,
186  UT_UniquePtr<fpreal[]> &data2,
187  const UT_Array<fpreal> &windower);
188  void computeWindow();
189 
190  static fpreal attenuate(fpreal value, fpreal power);
191 
192  int myWindowSize;
193 
194  fpreal *myFilter;
195  fpreal *myPhaseFilter;
196 
197  fpreal myLastDataSample;
198  fpreal myLastSlope;
199 
200  // Overlap is the region (specified as a fraction of the window size) of
201  // which neightboring chunks should overlap
202  fpreal myOverlap;
203  // Discard is the region (specified as a fraction of the window size) of
204  // data at the beginning and ending of each chunk which should be ignored
205  fpreal myDiscard;
206  // Whether or not the smoothed curves for neighboring chunks should be
207  // blended in the overlap region
208  bool myBlendChunks;
209  // Whether or not the phase filter should be applied, or just the magnitude
210  // filter. (phase filter must be non-null)
211  bool myFilterPhase;
212 
213  UT_Array<fpreal> myWindow;
214  UT_Array<fpreal> myUnWindow;
215 
216  UT_Array<fpreal> myReal;
217  UT_Array<fpreal> myImag;
218 
219  UT_FFT<fpreal> myTransform;
220 };
221 
222 template <class F, class FD>
223 void
225  fpreal *dest,
226  int dest_capacity,
227  fpreal freqstep,
228  const F &evaluate_data,
229  const FD &design_filter)
230 {
231  int j, index, n;
232  int garbage = int(myWindowSize*(myDiscard / 2.0));
233  int datachunk = int(myWindowSize - garbage*2);
234  int blendreg = int(datachunk * myOverlap);
235  int chunk;
236  UT_Interrupt *boss = UTgetInterrupt();
237  fpreal cph, sph;
238  fpreal v1, v2;
239  short int intcnt = 0;
240 
241  auto data1 = UTmakeUnique<fpreal[]>(myWindowSize);
242  auto data2 = UTmakeUnique<fpreal[]>(myWindowSize);
243  UT_UniquePtr<fpreal[]> blend = nullptr;
244 
245  if(blendreg >= datachunk)
246  blendreg = (int)SYSfloor(datachunk*0.95);
247 
248  auto readdata = UTmakeUnique<fpreal[]>(blendreg);
249 
250  if(datachunk-blendreg < dest_capacity)
251  {
252  blend = UTmakeUnique<fpreal[]>(blendreg);
253  cph = M_PI / fpreal(blendreg+1);
254  sph = -0.5 * (fpreal)M_PI;
255  for(j=0; j<blendreg; j++)
256  blend[j] = 0.5 + 0.5*SYSsin(sph + cph * fpreal(j+1));
257  }
258 
259  if (boss->opStart("Filtering Data"))
260  {
261  for(index=0, chunk=0; index<dest_capacity;
262  index+=(datachunk-blendreg)*2, chunk+=2)
263  {
264  if(!intcnt-- && boss->opInterrupt(100*index/dest_capacity))
265  break;
266 
267  // get data
268  if(myBlendChunks)
269  {
270  evaluate_data(fpreal(index), fpreal(index+datachunk-1),
271  data1.get() + garbage, datachunk);
272 
273  evaluate_data(fpreal(index+datachunk-blendreg),
274  fpreal(index-blendreg+datachunk*2-1),
275  data2.get() + garbage, datachunk);
276  }
277  else
278  {
279  evaluate_data(fpreal(index), fpreal(index+datachunk-1-blendreg),
280  data1.get() + garbage, datachunk-blendreg);
281  evaluate_data(fpreal(index+datachunk-blendreg),
282  fpreal(index-blendreg*2+datachunk*2-1),
283  data2.get() + garbage, datachunk-blendreg);
284 
285  for(j=blendreg; j--;)
286  {
287  data1[j+datachunk-blendreg] = 0;
288  data2[j+datachunk-blendreg] = 0;
289  }
290  }
291 
292  if(index>0)
293  {
294  for(j=blendreg; j--;)
295  data1[j+garbage] = readdata[j];
296  }
297 
298  for(j=blendreg; j--;)
299  readdata[j] = data2[garbage+datachunk-blendreg + j];
300 
301  for(j=garbage; j--;)
302  {
303  data1[j] =0;
304  data2[j] =0;
305  data1[myWindowSize-garbage+j] =0;
306  data2[myWindowSize-garbage+j] =0;
307  }
308 
309  window(data1, data2, myWindow);
310 
311  // do freq transform; multiply; do time transform
312  myTransform.toFrequency(
313  data1.get(), data2.get(), myReal.data(), myImag.data(),
314  myWindowSize);
315 
316  // Update the filter depending on the current index & chunk
317  design_filter(freqstep, index, chunk, myFilter, myPhaseFilter);
318 
319  if(!myFilterPhase)
320  {
321  // Magnitude scale only
322  for(j=0; j<myWindowSize; j++)
323  {
324  myReal[j] *= myFilter[j];
325  myImag[j] *= myFilter[j];
326  }
327  }
328  else
329  {
330  // Mag & phase filter
331  for(j=myWindowSize; j--;)
332  {
333  myReal[j] *= myFilter[j];
334  myImag[j] *= myPhaseFilter[j];
335  }
336  }
337 
338  myTransform.toTime(
339  myReal.data(), myImag.data(), data1.get(), data2.get(),
340  myWindowSize);
341 
342  window(data1, data2, myUnWindow);
343 
344  if(index>0)
345  {
346  // blend data part 1
347  n = blendreg;
348  if(index+blendreg >= dest_capacity)
349  {
350  if(index >= dest_capacity)
351  n = 0;
352  else
353  n = dest_capacity - index;
354  }
355 
356  if(myBlendChunks)
357  {
358  for(j=n; j--;)
359  {
360  v1 = dest[index+j]*(1-blend[j]);
361  v2 = data1[j+garbage] * blend[j];
362 
363  dest[index+j] = v2 + v1;
364  }
365  }
366  else
367  {
368  for(j=n; j--;)
369  dest[index+j] += data1[j+garbage];
370  }
371 
372  // store data part1 & 2
373  for(j=blendreg; j<datachunk && index+j<dest_capacity; j++)
374  dest[index+j] = data1[j+garbage];
375 
376  for(j=blendreg; j<datachunk && index+datachunk-blendreg+j<dest_capacity; j++)
377  dest[index+datachunk-blendreg+j] = data2[j+garbage];
378 
379  // blend data part 2
380  if(index+datachunk-blendreg < dest_capacity)
381  {
382  n = blendreg;
383  if(index+datachunk >= dest_capacity)
384  n = dest_capacity - (index+datachunk-blendreg);
385 
386  if(myBlendChunks)
387  {
388  for(j=n; j--;)
389  {
390  v1 = dest[index+datachunk-blendreg+j]*(1-blend[j]);
391  v2 = data2[j+garbage] * blend[j];
392 
393  dest[index+datachunk-blendreg+j] = v2 + v1;
394  }
395  }
396  else
397  {
398  for(j=n; j--;)
399  dest[index+datachunk-blendreg+j]+=data2[j+garbage];
400  }
401  }
402  }
403  else
404  {
405  // store data parts 1
406  n = datachunk;
407  if(index+datachunk >= dest_capacity)
408  {
409  if(index >= dest_capacity)
410  n = 0;
411  else
412  n = dest_capacity - index;
413  }
414 
415  for(j=n; j--;)
416  dest[index+j] = data1[j+garbage];
417 
418  // blend data part 2
419  if(index+datachunk-blendreg < dest_capacity)
420  {
421  n = blendreg;
422  if(index+datachunk >= dest_capacity)
423  n = dest_capacity - (index+datachunk-blendreg);
424 
425  if(myBlendChunks)
426  {
427  for(j=n; j--;)
428  {
429  v1 = dest[index+datachunk-blendreg+j]*(1-blend[j]);
430  v2 = data2[j+garbage] * blend[j];
431 
432  dest[index-blendreg+datachunk+j] = v2 + v1;
433  }
434  }
435  else
436  {
437  for(j=n; j--;)
438  dest[index-blendreg+datachunk+j] += data2[j+garbage];
439  }
440  }
441 
442  // store data part 2
443  for(j=blendreg; j<datachunk && index+datachunk-blendreg+j < dest_capacity; j++)
444  dest[index+datachunk-blendreg+j] = data2[j+garbage];
445 
446  }
447  }
448  }
449  boss->opEnd();
450 }
451 
452 template <class F>
453 void
455  fpreal *dest,
456  int dest_capacity,
457  fpreal freqstep,
458  const F &evaluate_data,
459  int constrain_region,
460  bool shift)
461 {
462  auto shifted = UTmakeUnique<fpreal[]>(dest_capacity);
463  evaluate_data(0, dest_capacity - 1, shifted.get(), dest_capacity);
464 
465  fpreal original_start = shifted[0];
466 
467  // shift so that start = 0
468  // this avoids the initial step response, especially when applied to constant functions
469  if (shift)
470  {
471  for (exint i = 0; i < dest_capacity; i++)
472  shifted[i] -= original_start;
473  }
474 
475  fpreal start = shifted[0];
476  fpreal end = shifted[dest_capacity - 1];
477 
478  auto evaluate_shifted = [&](int start, int, fpreal *data, int size)
479  {
480  for (exint i = 0; i < size; i++)
481  data[i] = (start + i) < dest_capacity ? shifted[start + i] : shifted[dest_capacity - 1];
482  };
483 
484  filter(dest, dest_capacity, freqstep, evaluate_shifted);
485 
486  if (constrain_region > dest_capacity / 2)
487  constrain_region = dest_capacity / 2;
488 
489  fpreal delta_start = start - dest[0];
490 
491  fpreal delta_end = end - dest[dest_capacity - 1];
492 
493  for (exint i = 0; i < dest_capacity / 2; i++)
494  {
495  fpreal t = static_cast<fpreal>(i) / constrain_region;
496  fpreal a = SYSpow(2.0, -10.0 * t);
497 
498  dest[i] += delta_start * a;
499 
500  dest[dest_capacity - 1 - i] += delta_end * a;
501  }
502 
503  // shift back to original range
504  if (shift)
505  {
506  for (exint i = 0; i < dest_capacity; i++)
507  dest[i] += original_start;
508  }
509 }
510 
511 // Designs a filter, populating the filter and phase_filter arrays
512 template <class F>
513 void
515  fpreal freq1, fpreal freq2,
516  fpreal attenuation, fpreal pass_gain, int order,
517  int chunk_size, fpreal freqstep, int index, int num,
518  const F& override_parameters,
519  fpreal *filter, fpreal *phase_filter)
520 {
521  fpreal f,k;
522  UT_Complex s,w1,w2,w3, w1_inv, response;
523  fpreal dc_gain,gain, exp1, filter_order;
524  int i;
525  int half = chunk_size/2;
526  fpreal max;
527 
528  exp1 = 2.0;
529  filter_order = static_cast<fpreal>(order); // Order of the filter. Used for Butterworth filters
530 
531  // Potentially overrides the frequency, gain, and myAttenuationuation parameters
532  // depending on the index and num
533  override_parameters(freq1, freq2, gain, attenuation, index, num);
534 
535  if(freq1<=SYS_FTOLERANCE)
536  freq1 = SYS_FTOLERANCE;
537  if(freq2<=SYS_FTOLERANCE)
538  freq2 = SYS_FTOLERANCE;
539 
540  switch(type)
541  {
542  case Type::LOW_PASS:
543  w2.set(1,0);
544  w1.set(freq2,0); // Cutoff frequency (high)
545  w1_inv = w2 / w1;
546 
547  s.set(0,freq2/100);
548 
549  k = SYSpow(freq2,exp1);
550 
551  if (design == Design::BUTTERWORTH)
552  // Frequency response of the Butterworth Filter
553  // See: https://en.wikipedia.org/wiki/Butterworth_filter#Transfer_function
554  response = (w2) / (w2 + (-(s * w1_inv).pow(exp1)).pow(filter_order));
555  else
556  response = (w2*k)/((s+w1).pow(exp1)); // default type
557 
558  dc_gain = response.magnitude();
559  if(SYSequalZero(dc_gain))
560  dc_gain=1.0;
561  gain = 1.0/dc_gain;
562 
563  for(i=0,f=0; i<=half; i++,f+=freqstep)
564  {
565  s.set(0,f);
566 
567  if (design == Design::BUTTERWORTH)
568  {
569  if (f < SYS_FTOLERANCE)
570  s.set(0, SYS_FTOLERANCE);
571  response = (w2) / (w2 + (-(s * w1_inv).pow(exp1)).pow(filter_order));
572  }
573  else
574  response = (w2*k)/((s+w1).pow(exp1));
575 
576  if(phase_filter)
577  {
578  filter[i] = response.real();
579  phase_filter[i] = response.imaginary();
580  }
581  else
582  filter[i] = response.magnitude()*gain;
583  }
584  break;
585 
586  case Type::HIGH_PASS:
587  w1.set(freq1,0);
588  w2.set(1, 0);
589  s.set(0,freq1*100);
590 
591  k = SYSpow(freq1,exp1);
592  if (design == Design::BUTTERWORTH)
593  response = (w2) / (w2 + (-(w1 / s).pow(exp1)).pow(filter_order));
594  else
595  response = (s.pow(exp1)*k) / (s+w1).pow(exp1);
596 
597  dc_gain = response.magnitude();
598  if(SYSequalZero(dc_gain))
599  dc_gain=1.0;
600 
601  gain = 1.0/dc_gain;
602 
603  for(i=0,f=0; i<=half; i++,f+=freqstep)
604  {
605  s.set(0,f);
606 
607  if (design == Design::BUTTERWORTH)
608  {
609  if (f < SYS_FTOLERANCE)
610  s.set(0, SYS_FTOLERANCE);
611 
612  response = (w2) / (w2 + (-(w1 / s).pow(exp1)).pow(filter_order));
613  }
614  else
615  response = (s.pow(exp1)*k)/(s+w1).pow(exp1);
616 
617  if(phase_filter)
618  {
619  filter[i] = response.real();
620  phase_filter[i] = response.imaginary();
621  }
622  else
623  filter[i] = response.magnitude()*gain;
624 
625  }
626  break;
627 
628  case Type::BAND_STOP:
629  case Type::BAND_PASS:
630 
631  if(freq1>freq2)
632  {
633  f=freq1;
634  freq1=freq2;
635  freq2=f;
636  }
637 
638  w1.set(freq1,0);
639  w2.set(freq2,0);
640  s.set(0,freq1);
641 
642  k = freq1 * freq2;
643  response = (s*k) /((s+w1) * (s+w2));
644 
645  dc_gain = response.magnitude();
646  if(SYSequalZero(dc_gain))
647  dc_gain=1.0;
648 
649  gain = 1.0/dc_gain;
650 
651  for(i=0,f=0; i<=half; i++,f+=freqstep)
652  {
653  s.set(0,f);
654 
655  response = (s*k)/((s+w1)*(s+w2));
656 
657  if(phase_filter)
658  {
659  filter[i] = response.real();
660  phase_filter[i] = response.imaginary();
661  }
662  else
663  filter[i] = response.magnitude()*gain;
664  }
665  break;
666  }
667 
668  if(!SYSisEqual(attenuation,1.0))
669  {
670  for(i=0; i<=half; i++)
671  filter[i] = SYSpow(filter[i],attenuation);
672  }
673 
674  gain = SYSpow(CONST_FPREAL(10), pass_gain/20);
675  max = 0;
676  for(i=0; i<=half; i++)
677  if(filter[i] > max)
678  max = filter[i];
679 
680  max = (max!=0)?(gain/max):gain;
681 
682  for(i=0; i<=half; i++)
683  filter[i] *=max;
684 
685  // band stop is just 1-bandpass
686  if(type == Type::BAND_STOP)
687  {
688  if(phase_filter)
689  {
690  for(i=0; i<=half; i++)
691  {
692  max = SYSsqrt(filter[i]*filter[i] + phase_filter[i]*phase_filter[i]);
693  max = (gain-max)/max;
694  filter[i]*= max;
695  phase_filter[i]*=max;
696  }
697  }
698  else
699  {
700  for(i=0; i<=half; i++)
701  filter[i] = gain - filter[i];
702  }
703  }
704 
705 
706  for(i=1; i<half; i++)
707  {
708  filter[chunk_size-i] = filter[i];
709  if(phase_filter)
710  phase_filter[chunk_size-i] = -phase_filter[i];
711  }
712 }
713 
714 template <class F>
715 void
716 CL_FreqFilter::designBandEQFilter(Shape shape, fpreal basefreq, fpreal *bands, int num_bands,
717  int chunk_size, fpreal freqstep, int index, int num,
718  const F& override_parameters,
719  fpreal *filter)
720 {
721  int i, j;
722  fpreal band;
723  int starti,endi;
724  fpreal freq,amp;
725  fpreal lstart, lend, lfreq, lwidth, lcenter;
726 
727  override_parameters(basefreq, bands, index, num);
728 
729  band = CONST_FPREAL(10) * SYSpow(CONST_FPREAL(2), basefreq);
730 
732 
733  for(i=0; i< chunk_size; i++)
734  filter[i] = 1.0;
735 
736  endi = 0;
737  for(i=0; i<num_bands; i++)
738  {
739  if(endi>=chunk_size/2)
740  continue;
741 
742  if(shape == Shape::BOX)
743  starti = endi;
744  else
745  starti = int((band/width)/freqstep + 0.5);
746  endi = int((band*width)/freqstep + 1.5);
747 
748  if(bands[i] != 1.0)
749  {
750  lstart = SYSlog(band/width);
751  lend = SYSlog(band*width);
752  lwidth = (lend-lstart)/2.0;
753  lcenter = SYSlog(band);
754 
755  // filter contraints.
756  if(starti>chunk_size/2)
757  continue;
758  if(endi>chunk_size/2)
759  endi = chunk_size/2;
760  if(starti<0)
761  starti=0;
762 
763 
764  for(j=starti; j<=endi; j++)
765  {
766  lfreq = SYSlog(j*freqstep);
767  if(shape == Shape::GAUSSIAN)
768  {
769  freq = 3.0*(lfreq-lcenter)/lwidth;
770 
771  amp = SYSexp( -0.5 *freq*freq);
772  amp = (bands[i] - 1.0) * amp;
773 
774  filter[j] +=amp;
775  }
776  else
777  {
778  filter[j] = bands[i];
779  }
780  }
781  }
782 
783  band *=2;
784  }
785 
786  for(i=1; i<chunk_size/2; i++)
787  filter[chunk_size-i] = filter[i];
788 
789 }
790 
791 template <class F>
792 void
794  fpreal width, fpreal atten,
795  int chunk_size, fpreal freqstep, int index, int num,
796  const F& override_parameters,
797  fpreal *filter)
798 {
799  int i, j;
800  int starti,endi;
801  fpreal amp,ind;
802  fpreal lstart, lend, lfreq, lwidth, lcenter;
803 
804  override_parameters(freq, width, fgain, pgain, shape, atten, index, num);
805 
806  for(i=0; i<chunk_size; i++)
807  filter[i] = pgain;
808 
809  width = SYSpow(2.0,width/2.0);
810 
811  starti = int((freq/width)/freqstep + 0.5);
812  endi = int((freq*width)/freqstep + 1.5);
813 
814  // filter contraints.
815  if(starti>chunk_size/2)
816  return;
817  if(endi>chunk_size/2)
818  endi = chunk_size/2;
819  if(starti<0)
820  starti=0;
821 
822  lstart = SYSlog(freq/width);
823  lend = SYSlog(freq*width);
824  lcenter = SYSlog(freq);
825  lwidth = (lend - lstart)/2.0;
826 
827  ind = freq/width;
828  for(j=starti; j<=endi; j++, ind+=freqstep)
829  {
830  lfreq = SYSlog(ind);
831 
832  if(shape <= 0.5)
833  {
834  // blend between box (0) and triangle (0.5)
835  amp = (1.0-SYSabs(lfreq-lcenter)/lwidth);
836 
837  amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
838  amp = fgain*(0.5-shape)*2 + amp*(shape*2);
839 
840  filter[j] =amp;
841  }
842  else
843  {
844  // blend between triangle (0.5) and gaussian (1.0)
845  freq = 3.0*(lfreq-lcenter)/lwidth;
846 
847  amp = SYSexp( -0.5 *freq*freq);
848  amp = (fgain - pgain)*attenuate(amp,atten) + pgain;
849  filter[j] = amp*(shape-0.5)*2;
850 
851  amp = 1.0-SYSabs(lfreq-lcenter)/lwidth;
852  amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
853  filter[j] += amp*(1.0-shape)*2;
854  }
855  }
856 
857  for(i=1; i<chunk_size/2; i++)
858  filter[chunk_size-i] = filter[i];
859 
860 }
861 
862 template <class F>
863 void
864 CL_FreqFilter::designAverageSideband(fpreal fgain, fpreal gain, fpreal effect, int length, int passfreq,
865  int chunk_size, const F& evaluate_data,
866  fpreal *filter)
867 {
868  UT_FFT<fpreal> freqtransform(chunk_size);
869 
870  int i,j,count;
871  fpreal factor;
872  int start,end;
873  double max;
874 
875  auto temp = UTmakeUnique<fpreal[]>(chunk_size);
876  auto real = UTmakeUnique<fpreal[]>(chunk_size);
877  auto imag = UTmakeUnique<fpreal[]>(chunk_size);
878 
879  for(j=1; j<chunk_size; j++)
880  filter[j]=0;
881 
882 
883  for(i=0,count=0; i<length; i+=chunk_size,count++)
884  {
885  evaluate_data(
886  (fpreal)i, (fpreal)i + chunk_size - 1, temp.get(), chunk_size);
887 
888  // do freq transform; multiply; do time transform
889  freqtransform.toFrequency(temp.get(),real.get(),imag.get(),chunk_size);
890 
891  // calculate power spectrum
892  for(j=1; j<chunk_size; j++)
893  temp[j] = (real[j]*real[j] + imag[j]*imag[j]);
894 
895  for(j=1; j<chunk_size; j++)
896  filter[j] += temp[j];
897  }
898 
899  if(count>1)
900  for(j=1; j<chunk_size; j++)
901  filter[j] /= count;
902 
903  // 'effect' widens the sideband filter on a sample basis
904  if(!SYSisEqual(effect,1.0))
905  {
906  for(j=(chunk_size>>1); j--;)
907  {
908  temp[j] = filter[j];
909  filter[j] = 0;
910  }
911 
912  for(j=(chunk_size>>1); j--;)
913  {
914  start = int(j/effect);
915  end = int(j*effect);
916 
917  if(start < 0)
918  start=0;
919  if(end > chunk_size/2)
920  end = chunk_size/2;
921 
922  for(i=start; i<=end; i++)
923  {
924  // geometric interp
925  factor = SYSlog(fpreal(i)/fpreal(j))/SYSlog(effect);
926 
927  // gaussian distribution
928  factor = SYSexp(-1.5*factor*factor);
929  filter[i] += temp[j] * factor;
930  }
931  }
932 
933  for(j=1; j<(chunk_size>>1); j++)
934  filter[chunk_size-j] = filter[j];
935  }
936  // don't filter DC
937 // filter[0] = 1.0;
938 
939  for(j=0; j<chunk_size; j++)
940  if(j==0 || filter[j]>max)
941  max = filter[j];
942 
943  if(!SYSequalZero(max))
944  for(j=0; j<chunk_size; j++)
945  filter[j] /=max;
946 
947  if(!SYSisEqual(fgain,1.0))
948  for(j=0; j<chunk_size; j++)
949  filter[j] *= fgain;
950 
951  if(passfreq)
952  {
953  for(j=0; j<chunk_size; j++)
954  filter[j] = (1.0+gain) - filter[j];
955  }
956  else
957  {
958  for(j=0; j<chunk_size; j++)
959  filter[j] += gain;
960  }
961 }
962 
963 template <class F>
964 void
966  int chunk_size, int index,
967  const F& evaluate_data,
968  fpreal *filter)
969 {
970  UT_FFT<fpreal> freqtransform(chunk_size);
971 
972  int i,j;
973  int start,end;
974  fpreal factor;
975  double max;
976 
977  auto real = UTmakeUnique<fpreal[]>(chunk_size);
978  auto imag = UTmakeUnique<fpreal[]>(chunk_size);
979 
980  evaluate_data((fpreal)index, (fpreal)index+chunk_size-1, filter, chunk_size);
981 
982  // do freq transform; multiply; do time transform
983  freqtransform.toFrequency(filter,real.get(),imag.get(),chunk_size);
984 
985  // calculate power spectrum
986  for(j=1; j<chunk_size; j++)
987  filter[j] = (real[j]*real[j] + imag[j]*imag[j]);
988 
989  // 'effect' widens the sideband filter on a sample basis
990  if(!SYSisEqual(effect,1.0))
991  {
992  for(j=0; j<=chunk_size/2; j++)
993  {
994  real[j] = filter[j];
995  filter[j] = 0;
996  }
997 
998  for(j=0; j<chunk_size/2; j++)
999  {
1000  start = int(j/effect);
1001  end = int(j*effect);
1002 
1003  if(start < 0)
1004  start=0;
1005  if(end > chunk_size/2)
1006  end = chunk_size/2;
1007 
1008  for(i=start; i<=end; i++)
1009  {
1010  // geometric interp
1011  factor = SYSlog(fpreal(i)/fpreal(j))/SYSlog(effect);
1012 
1013  // gaussian distribution
1014  factor = SYSexp(-1.5*factor*factor);
1015  filter[i] += real[j] * factor;
1016  }
1017  }
1018 
1019  for(j=1; j<chunk_size/2; j++)
1020  filter[chunk_size-j] = filter[j];
1021  }
1022 
1023  // don't filter DC
1024 // filter[0]=1;
1025 
1026  for(j=0; j<chunk_size; j++)
1027  if(j==0 || filter[j]>max)
1028  max = filter[j];
1029 
1030  if(!SYSequalZero(max))
1031  for(j=0; j<chunk_size; j++)
1032  filter[j] /=max;
1033 
1034  if(!SYSisEqual(fgain,1.0))
1035  for(j=0; j<chunk_size; j++)
1036  filter[j] *= fgain;
1037 
1038  if(passfreq)
1039  {
1040  for(j=0; j<chunk_size; j++)
1041  filter[j] = (gain+1.0) - filter[j];
1042  }
1043  else
1044  {
1045  for(j=0; j<chunk_size; j++)
1046  filter[j] += gain;
1047  }
1048 }
1049 
1050 #endif
bool isReady() const
Ensures the filter has been set and the window size is valid.
Definition: CL_FreqFilter.h:66
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
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)
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
GLboolean * data
Definition: glcorearb.h:131
#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
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
GLdouble s
Definition: glad.h:3009
#define SYSabs(a)
Definition: SYS_Math.h:1572
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
void setBlendChunks(bool blend)
Definition: CL_FreqFilter.h:75
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:818
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)
std::unique_ptr< T, Deleter > UT_UniquePtr
A smart pointer for unique ownership of dynamically allocated objects.
Definition: UT_UniquePtr.h:39
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
GLdouble n
Definition: glcorearb.h:2008
GLfloat f
Definition: glcorearb.h:1926
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)
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
GLuint GLuint end
Definition: glcorearb.h:475
GLdouble GLdouble GLint GLint order
Definition: glad.h:2676
void setDiscard(fpreal discard)
Definition: CL_FreqFilter.h:74
SYS_API fpreal32 SYSfloor(fpreal32 val)
void toFrequency(T *source, T *destr, T *desti, exint npnts)
GLdouble t
Definition: glad.h:2397
T imaginary() const
Definition: UT_Complex.h:70
bool SYSequalZero(const UT_Vector3T< T > &v)
Definition: UT_Vector3.h:1069
GLint j
Definition: glad.h:2733
GLsizeiptr size
Definition: glcorearb.h:664
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:842
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
GLuint index
Definition: glcorearb.h:786
void opEnd(int opid=-1)
GLfloat GLfloat v1
Definition: glcorearb.h:817
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
void setFilterPhase(bool filter_phase)
Definition: CL_FreqFilter.h:76
GLint GLsizei width
Definition: glcorearb.h:103
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
bool SYSisEqual(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, S tol=SYS_FTOLERANCE)
Componentwise equality.
Definition: UT_Vector2.h:674
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
type
Definition: core.h:1059
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...
GLint GLint GLint GLint GLint GLint GLint GLbitfield GLenum filter
Definition: glcorearb.h:1297