13 #ifndef __CL_FreqFilter__
14 #define __CL_FreqFilter__
24 #define SQRT_TWO 1.4142136
56 bool blend_chunks =
true,
57 bool filter_phase =
false);
67 if (myFilterPhase && myPhaseFilter ==
nullptr)
70 return myFilter !=
nullptr && myWindowSize > 0;
81 static void designPassFilter(
Type type, Design design,
85 const F& override_parameters,
97 auto override_parameters
101 type, design, freq1, freq2, attenuation, pass_gain, order,
102 size, freqstep, 0, 0, override_parameters,
103 mag_filter, phase_filter);
109 static void designBandEQFilter(Shape shape,
fpreal basefreq,
fpreal *bands,
int num_bands,
111 const F& override_parameters,
120 const F& override_parameters,
127 int size,
const F& evaluate_data,
133 static void designContinuousSideband(
fpreal fgain,
fpreal gain,
fpreal effect,
int passfreq,
134 int size,
int index,
const F& evaluate_data,
142 template <
class F,
class FD>
144 const F& evaluate_data,
const FD& design_filter);
152 const F &evaluate_data)
156 filter(dest, dest_capacity, freqstep, evaluate_data, design_filter);
175 void constrainedFilter(
179 const F &evaluate_data,
180 int constrain_region,
188 void computeWindow();
222 template <
class F,
class FD>
228 const F &evaluate_data,
229 const FD &design_filter)
232 int garbage =
int(myWindowSize*(myDiscard / 2.0));
233 int datachunk =
int(myWindowSize - garbage*2);
234 int blendreg =
int(datachunk * myOverlap);
239 short int intcnt = 0;
241 auto data1 = UTmakeUnique<fpreal[]>(myWindowSize);
242 auto data2 = UTmakeUnique<fpreal[]>(myWindowSize);
245 if(blendreg >= datachunk)
248 auto readdata = UTmakeUnique<fpreal[]>(blendreg);
250 if(datachunk-blendreg < dest_capacity)
252 blend = UTmakeUnique<fpreal[]>(blendreg);
255 for(j=0; j<blendreg; j++)
256 blend[j] = 0.5 + 0.5*SYSsin(sph + cph *
fpreal(j+1));
259 if (boss->
opStart(
"Filtering Data"))
261 for(index=0, chunk=0; index<dest_capacity;
262 index+=(datachunk-blendreg)*2, chunk+=2)
264 if(!intcnt-- && boss->
opInterrupt(100*index/dest_capacity))
271 data1.get() + garbage, datachunk);
273 evaluate_data(
fpreal(index+datachunk-blendreg),
274 fpreal(index-blendreg+datachunk*2-1),
275 data2.get() + garbage, datachunk);
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);
285 for(j=blendreg; j--;)
287 data1[j+datachunk-blendreg] = 0;
288 data2[j+datachunk-blendreg] = 0;
294 for(j=blendreg; j--;)
295 data1[j+garbage] = readdata[j];
298 for(j=blendreg; j--;)
299 readdata[j] = data2[garbage+datachunk-blendreg + j];
305 data1[myWindowSize-garbage+
j] =0;
306 data2[myWindowSize-garbage+
j] =0;
309 window(data1, data2, myWindow);
313 data1.get(), data2.get(), myReal.
data(), myImag.
data(),
317 design_filter(freqstep, index, chunk, myFilter, myPhaseFilter);
322 for(j=0; j<myWindowSize; j++)
324 myReal[
j] *= myFilter[
j];
325 myImag[
j] *= myFilter[
j];
331 for(j=myWindowSize; j--;)
333 myReal[
j] *= myFilter[
j];
334 myImag[
j] *= myPhaseFilter[
j];
339 myReal.
data(), myImag.
data(), data1.get(), data2.get(),
342 window(data1, data2, myUnWindow);
348 if(index+blendreg >= dest_capacity)
350 if(index >= dest_capacity)
353 n = dest_capacity -
index;
360 v1 = dest[index+
j]*(1-blend[
j]);
361 v2 = data1[j+garbage] * blend[
j];
363 dest[index+
j] = v2 +
v1;
369 dest[index+j] += data1[j+garbage];
373 for(j=blendreg; j<datachunk && index+j<dest_capacity; j++)
374 dest[index+j] = data1[j+garbage];
376 for(j=blendreg; j<datachunk && index+datachunk-blendreg+j<dest_capacity; j++)
377 dest[index+datachunk-blendreg+j] = data2[j+garbage];
380 if(index+datachunk-blendreg < dest_capacity)
383 if(index+datachunk >= dest_capacity)
384 n = dest_capacity - (index+datachunk-blendreg);
390 v1 = dest[index+datachunk-blendreg+
j]*(1-blend[
j]);
391 v2 = data2[j+garbage] * blend[
j];
393 dest[index+datachunk-blendreg+
j] = v2 +
v1;
399 dest[index+datachunk-blendreg+j]+=data2[j+garbage];
407 if(index+datachunk >= dest_capacity)
409 if(index >= dest_capacity)
412 n = dest_capacity -
index;
416 dest[index+j] = data1[j+garbage];
419 if(index+datachunk-blendreg < dest_capacity)
422 if(index+datachunk >= dest_capacity)
423 n = dest_capacity - (index+datachunk-blendreg);
429 v1 = dest[index+datachunk-blendreg+
j]*(1-blend[
j]);
430 v2 = data2[j+garbage] * blend[
j];
432 dest[index-blendreg+datachunk+
j] = v2 +
v1;
438 dest[index-blendreg+datachunk+j] += data2[j+garbage];
443 for(j=blendreg; j<datachunk && index+datachunk-blendreg+j < dest_capacity; j++)
444 dest[index+datachunk-blendreg+j] = data2[j+garbage];
458 const F &evaluate_data,
459 int constrain_region,
462 auto shifted = UTmakeUnique<fpreal[]>(dest_capacity);
463 evaluate_data(0, dest_capacity - 1, shifted.get(), dest_capacity);
465 fpreal original_start = shifted[0];
471 for (
exint i = 0; i < dest_capacity; i++)
472 shifted[i] -= original_start;
481 data[i] = (start + i) < dest_capacity ? shifted[start + i] : shifted[dest_capacity - 1];
484 filter(dest, dest_capacity, freqstep, evaluate_shifted);
486 if (constrain_region > dest_capacity / 2)
487 constrain_region = dest_capacity / 2;
489 fpreal delta_start = start - dest[0];
491 fpreal delta_end = end - dest[dest_capacity - 1];
493 for (
exint i = 0; i < dest_capacity / 2; i++)
496 fpreal a = SYSpow(2.0, -10.0 * t);
498 dest[i] += delta_start *
a;
500 dest[dest_capacity - 1 - i] += delta_end *
a;
506 for (
exint i = 0; i < dest_capacity; i++)
507 dest[i] += original_start;
517 int chunk_size,
fpreal freqstep,
int index,
int num,
518 const F& override_parameters,
523 fpreal dc_gain,gain, exp1, filter_order;
525 int half = chunk_size/2;
533 override_parameters(freq1, freq2, gain, attenuation, index, num);
549 k = SYSpow(freq2,exp1);
554 response = (w2) / (w2 + (-(s * w1_inv).
pow(exp1)).
pow(filter_order));
556 response = (w2*k)/((s+w1).pow(exp1));
563 for(i=0,f=0; i<=
half; i++,f+=freqstep)
571 response = (w2) / (w2 + (-(s * w1_inv).
pow(exp1)).
pow(filter_order));
574 response = (w2*k)/((s+w1).pow(exp1));
578 filter[i] = response.
real();
591 k = SYSpow(freq1,exp1);
593 response = (w2) / (w2 + (-(w1 / s).
pow(exp1)).
pow(filter_order));
595 response = (s.
pow(exp1)*k) / (s+w1).
pow(exp1);
603 for(i=0,f=0; i<=
half; i++,f+=freqstep)
612 response = (w2) / (w2 + (-(w1 / s).
pow(exp1)).
pow(filter_order));
615 response = (s.
pow(exp1)*k)/(s+w1).
pow(exp1);
619 filter[i] = response.
real();
643 response = (s*k) /((s+w1) * (s+w2));
651 for(i=0,f=0; i<=
half; i++,f+=freqstep)
655 response = (s*k)/((s+w1)*(s+w2));
659 filter[i] = response.
real();
670 for(i=0; i<=
half; i++)
671 filter[i] = SYSpow(filter[i],attenuation);
676 for(i=0; i<=
half; i++)
680 max = (max!=0)?(gain/max):gain;
682 for(i=0; i<=
half; i++)
690 for(i=0; i<=
half; i++)
692 max = SYSsqrt(filter[i]*filter[i] + phase_filter[i]*phase_filter[i]);
693 max = (gain-
max)/max;
695 phase_filter[i]*=
max;
700 for(i=0; i<=
half; i++)
701 filter[i] = gain - filter[i];
706 for(i=1; i<
half; i++)
708 filter[chunk_size-i] = filter[i];
710 phase_filter[chunk_size-i] = -phase_filter[i];
717 int chunk_size,
fpreal freqstep,
int index,
int num,
718 const F& override_parameters,
725 fpreal lstart, lend, lfreq, lwidth, lcenter;
727 override_parameters(basefreq, bands, index, num);
733 for(i=0; i< chunk_size; i++)
737 for(i=0; i<num_bands; i++)
739 if(endi>=chunk_size/2)
745 starti =
int((band/width)/freqstep + 0.5);
746 endi =
int((band*width)/freqstep + 1.5);
750 lstart = SYSlog(band/width);
751 lend = SYSlog(band*width);
752 lwidth = (lend-lstart)/2.0;
753 lcenter = SYSlog(band);
756 if(starti>chunk_size/2)
758 if(endi>chunk_size/2)
764 for(j=starti; j<=endi; j++)
766 lfreq = SYSlog(j*freqstep);
769 freq = 3.0*(lfreq-lcenter)/lwidth;
771 amp = SYSexp( -0.5 *freq*freq);
772 amp = (bands[i] - 1.0) * amp;
778 filter[
j] = bands[i];
786 for(i=1; i<chunk_size/2; i++)
787 filter[chunk_size-i] = filter[i];
795 int chunk_size,
fpreal freqstep,
int index,
int num,
796 const F& override_parameters,
802 fpreal lstart, lend, lfreq, lwidth, lcenter;
804 override_parameters(freq, width, fgain, pgain, shape, atten, index, num);
806 for(i=0; i<chunk_size; i++)
809 width = SYSpow(2.0,width/2.0);
811 starti =
int((freq/width)/freqstep + 0.5);
812 endi =
int((freq*width)/freqstep + 1.5);
815 if(starti>chunk_size/2)
817 if(endi>chunk_size/2)
822 lstart = SYSlog(freq/width);
823 lend = SYSlog(freq*width);
824 lcenter = SYSlog(freq);
825 lwidth = (lend - lstart)/2.0;
828 for(j=starti; j<=endi; j++, ind+=freqstep)
835 amp = (1.0-
SYSabs(lfreq-lcenter)/lwidth);
837 amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
838 amp = fgain*(0.5-shape)*2 + amp*(shape*2);
845 freq = 3.0*(lfreq-lcenter)/lwidth;
847 amp = SYSexp( -0.5 *freq*freq);
848 amp = (fgain - pgain)*attenuate(amp,atten) + pgain;
849 filter[
j] = amp*(shape-0.5)*2;
851 amp = 1.0-
SYSabs(lfreq-lcenter)/lwidth;
852 amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
853 filter[
j] += amp*(1.0-shape)*2;
857 for(i=1; i<chunk_size/2; i++)
858 filter[chunk_size-i] = filter[i];
865 int chunk_size,
const F& evaluate_data,
875 auto temp = UTmakeUnique<fpreal[]>(chunk_size);
876 auto real = UTmakeUnique<fpreal[]>(chunk_size);
877 auto imag = UTmakeUnique<fpreal[]>(chunk_size);
879 for(j=1; j<chunk_size; j++)
883 for(i=0,count=0; i<
length; i+=chunk_size,count++)
886 (
fpreal)i, (
fpreal)i + chunk_size - 1, temp.get(), chunk_size);
889 freqtransform.
toFrequency(temp.get(),real.get(),imag.get(),chunk_size);
892 for(j=1; j<chunk_size; j++)
893 temp[j] = (real[j]*real[j] + imag[j]*imag[j]);
895 for(j=1; j<chunk_size; j++)
896 filter[j] += temp[j];
900 for(j=1; j<chunk_size; j++)
906 for(j=(chunk_size>>1); j--;)
912 for(j=(chunk_size>>1); j--;)
914 start =
int(j/effect);
919 if(end > chunk_size/2)
922 for(i=start; i<=
end; i++)
928 factor = SYSexp(-1.5*factor*factor);
929 filter[i] += temp[
j] * factor;
934 filter[chunk_size-j] = filter[j];
939 for(j=0; j<chunk_size; j++)
940 if(j==0 || filter[j]>max)
944 for(j=0; j<chunk_size; j++)
948 for(j=0; j<chunk_size; j++)
953 for(j=0; j<chunk_size; j++)
954 filter[j] = (1.0+gain) - filter[
j];
958 for(j=0; j<chunk_size; j++)
966 int chunk_size,
int index,
967 const F& evaluate_data,
977 auto real = UTmakeUnique<fpreal[]>(chunk_size);
978 auto imag = UTmakeUnique<fpreal[]>(chunk_size);
980 evaluate_data((
fpreal)index, (
fpreal)index+chunk_size-1, filter, chunk_size);
983 freqtransform.
toFrequency(filter,real.get(),imag.get(),chunk_size);
986 for(j=1; j<chunk_size; j++)
987 filter[j] = (real[j]*real[j] + imag[j]*imag[j]);
992 for(j=0; j<=chunk_size/2; j++)
998 for(j=0; j<chunk_size/2; j++)
1000 start =
int(j/effect);
1001 end =
int(j*effect);
1005 if(end > chunk_size/2)
1008 for(i=start; i<=
end; i++)
1014 factor = SYSexp(-1.5*factor*factor);
1015 filter[i] += real[
j] * factor;
1019 for(j=1; j<chunk_size/2; j++)
1020 filter[chunk_size-j] = filter[j];
1026 for(j=0; j<chunk_size; j++)
1027 if(j==0 || filter[j]>max)
1031 for(j=0; j<chunk_size; j++)
1035 for(j=0; j<chunk_size; j++)
1040 for(j=0; j<chunk_size; j++)
1041 filter[j] = (gain+1.0) - filter[
j];
1045 for(j=0; j<chunk_size; j++)
bool isReady() const
Ensures the filter has been set and the window size is valid.
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
void toTime(T *sourcer, T *sourcei, T *dest, exint npnts)
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
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)
GLboolean GLboolean GLboolean GLboolean a
GLuint GLsizei GLsizei * length
void setBlendChunks(bool blend)
GLfloat GLfloat GLfloat v2
UT_ComplexT pow(T exp) const
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.
ImageBuf OIIO_API pow(const ImageBuf &A, cspan< float > B, ROI roi={}, int nthreads=0)
vint4 blend(const vint4 &a, const vint4 &b, const vbool4 &mask)
int opInterrupt(int percent=-1)
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)
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)
GLdouble GLdouble GLint GLint order
void setDiscard(fpreal discard)
SYS_API fpreal32 SYSfloor(fpreal32 val)
void toFrequency(T *source, T *destr, T *desti, exint npnts)
bool SYSequalZero(const UT_Vector3T< T > &v)
ImageBuf OIIO_API resize(const ImageBuf &src, string_view filtername="", float filterwidth=0.0f, ROI roi={}, int nthreads=0)
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
void setFilterPhase(bool filter_phase)
int opStart(const char *opname=0, int npoll=0, int immediate_dialog=0, int *opid=0)
bool SYSisEqual(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, S tol=SYS_FTOLERANCE)
Componentwise equality.
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)
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