CL_FreqFilter.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * NAME:        CL_FreqFilter.h ( Clip Library, C++)
 *
 * COMMENTS:
 *    Does continuous digital filtering in the frequency domain
 */


#ifndef __CL_FreqFilter__
#define __CL_FreqFilter__

#include "CL_API.h"
#include <UT/UT_Array.h>
#include <UT/UT_Complex.h>
#include <UT/UT_FFT.h>
#include <UT/UT_Interrupt.h>
#include <UT/UT_UniquePtr.h>
#include <SYS/SYS_Types.h>

#define SQRT_TWO 1.4142136

class CL_API CL_FreqFilter final
{
public:

    // Filter Type, used for Pass filters
    enum class Type
    {
        LOW_PASS,
        HIGH_PASS,
        BAND_PASS,
        BAND_STOP
    };

    // Filter Design, used for Pass filters
    enum class Design
    {
        STANDARD,
        BUTTERWORTH
    };

    // Filter shape, used for Equalization Filters
    enum class Shape
    {
        GAUSSIAN,
        BOX
    };

    CL_FreqFilter(int size, 
                  fpreal overlap = 0.1,
                  fpreal discard = 0.1,
                  bool blend_chunks = true,
                  bool filter_phase = false);

    ~CL_FreqFilter();

    void resize(int size);

    void setFilter(int size, fpreal *filter, fpreal *phase = nullptr);

    /// Ensures the filter has been set and the window size is valid
    bool isReady() const { 
        if (myFilterPhase && myPhaseFilter == nullptr)
            return false;

        return myFilter != nullptr && myWindowSize > 0;
    }

    void setOverlap(fpreal overlap) { myOverlap = overlap; }
    void setDiscard(fpreal discard) { myDiscard = discard; }
    void setBlendChunks(bool blend) { myBlendChunks = blend; }
    void setFilterPhase(bool filter_phase) { myFilterPhase = filter_phase; }

    /// Designs a pass filter, filling out the filter and phase_filter arrays
    /// Filter MUST be size floats long.
    template <class F>
    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);

    /// Variant of the designPassFilter method which does not have a parameter
    /// override method
    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)
    {
        // Override Parameters does nothing
        auto override_parameters
                = [](fpreal &, fpreal &, fpreal &, fpreal &, int, int) {};

        CL_FreqFilter::designPassFilter(
                type, design, freq1, freq2, attenuation, pass_gain, order,
                size, freqstep, 0, 0, override_parameters,
                mag_filter, phase_filter);
    }

    // Designs a eq filter, filling out the filter array
    // Filter MUST be size floats long.
    template <class F>
    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);

    // Designs a eq filter, filling out the filter array
    // Filter MUST be size floats long.
    template <class F>
    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);

    // Designs an average sideband filter, filling out the filter and phase_filter arrays
    // Filter MUST be size floats long.
    template <class F>
    static void designAverageSideband(fpreal fgain, fpreal gain, fpreal effect, int length, int passfreq,
                                   int size, const F& evaluate_data, 
                                   fpreal *filter);

    // Designs an average sideband filter, filling out the filter and phase_filter arrays
    // Filter MUST be size floats long.
    template <class F>
    static void designContinuousSideband(fpreal fgain, fpreal gain, fpreal effect, int passfreq,
                                   int size, int index, const F& evaluate_data, 
                                   fpreal *filter);

    /// Filter input data by applying the frequency domain filter
    /// @param evaluate_data lambda function which returns the data value at the
    ///                      given index
    /// @param design_filter lambda function which redesigns the filter at the
    ///                      currently processed chunk.
    template <class F, class FD>
    void filter(fpreal* dest, int dest_capacity, fpreal freqstep,
                const F& evaluate_data, const FD& design_filter);

    /// Variant of the filter method which does not need to redesign the filter depending on the chunk
    template <class F>
    void filter(
            fpreal *dest,
            int dest_capacity,
            fpreal freqstep,
            const F &evaluate_data)
    {
        // Don't redesign filter depending on the chunk
        auto design_filter = [](fpreal, int, int, fpreal *, fpreal *) {};
        filter(dest, dest_capacity, freqstep, evaluate_data, design_filter);
    }

    /// Filter method which constrains the endpoints by "bending" the smoothed
    /// curve to match the start and end point of the original curve.
    ///
    /// @param constrain_region - Number of samples over which to bend the curve
    /// at the beginning and end of the curve. Setting this equal to the sample
    /// rate, attenuated by the high frequency cutoff tends to give good 
    /// results for smoothing animation curves.
    /// Example: constrain_region = SYSmax((1.0 - high_cutoff) * sample_rate), 1)
    /// Setting this to a value <= 0 will skip the constraining step.
    /// 
    /// @param shift - If true, will apply a y-offset to the input data so the
    /// starting value is zero before smoothing. The smoothed curve will then be
    /// shifted back after smoothing. This helps avoid the step response which
    /// occurs when smoothing data with a non-zero starting y-value. This can be
    /// especially seen when applying the filter to a constant function
    template <class F>
    void constrainedFilter(
            fpreal *dest,
            int dest_capacity,
            fpreal freqstep,
            const F &evaluate_data,
            int constrain_region,
            bool shift = true);

private:
    void window(
            UT_UniquePtr<fpreal[]> &data1,
            UT_UniquePtr<fpreal[]> &data2,
            const UT_Array<fpreal> &windower);
    void computeWindow();

    static fpreal attenuate(fpreal value, fpreal power);

    int myWindowSize;

    fpreal *myFilter;
    fpreal *myPhaseFilter;

    fpreal myLastDataSample;
    fpreal myLastSlope;

    // Overlap is the region (specified as a fraction of the window size) of
    // which neightboring chunks should overlap
    fpreal myOverlap;
    // Discard is the region (specified as a fraction of the window size) of
    // data at the beginning and ending of each chunk which should be ignored
    fpreal myDiscard;
    // Whether or not the smoothed curves for neighboring chunks should be
    // blended in the overlap region
    bool myBlendChunks;
    // Whether or not the phase filter should be applied, or just the magnitude
    // filter. (phase filter must be non-null)
    bool myFilterPhase;

    UT_Array<fpreal> myWindow;
    UT_Array<fpreal> myUnWindow;

    UT_Array<fpreal> myReal;
    UT_Array<fpreal> myImag;

    UT_FFT<fpreal> myTransform;
};

template <class F, class FD>
void
CL_FreqFilter::filter(
        fpreal *dest,
        int dest_capacity,
        fpreal freqstep,
        const F &evaluate_data,
        const FD &design_filter)
{
    int j, index, n;
    int garbage = int(myWindowSize*(myDiscard / 2.0));
    int datachunk = int(myWindowSize - garbage*2);
    int blendreg = int(datachunk * myOverlap);
    int chunk;
    UT_Interrupt *boss = UTgetInterrupt();
    fpreal cph, sph;
    fpreal v1, v2;
    short int intcnt = 0;

    auto data1 = UTmakeUnique<fpreal[]>(myWindowSize);
    auto data2 = UTmakeUnique<fpreal[]>(myWindowSize);
    UT_UniquePtr<fpreal[]> blend = nullptr;

    if(blendreg >= datachunk)
        blendreg = (int)SYSfloor(datachunk*0.95);

    auto readdata = UTmakeUnique<fpreal[]>(blendreg);

    if(datachunk-blendreg < dest_capacity)
    {
        blend = UTmakeUnique<fpreal[]>(blendreg);
        cph = M_PI / fpreal(blendreg+1);
        sph = -0.5 * (fpreal)M_PI;
        for(j=0; j<blendreg; j++)
            blend[j] = 0.5 + 0.5*SYSsin(sph + cph * fpreal(j+1));
    }

    if (boss->opStart("Filtering Data"))
    {
        for(index=0, chunk=0; index<dest_capacity;
            index+=(datachunk-blendreg)*2, chunk+=2) 
        {
             if(!intcnt-- && boss->opInterrupt(100*index/dest_capacity))
                break;
         
             // get data
             if(myBlendChunks)
             {
                 evaluate_data(fpreal(index), fpreal(index+datachunk-1),
                               data1.get() + garbage, datachunk);

                 evaluate_data(fpreal(index+datachunk-blendreg), 
                               fpreal(index-blendreg+datachunk*2-1),
                               data2.get() + garbage, datachunk);
             }
             else
             {
                 evaluate_data(fpreal(index), fpreal(index+datachunk-1-blendreg),
                               data1.get() + garbage, datachunk-blendreg);
                 evaluate_data(fpreal(index+datachunk-blendreg),
                               fpreal(index-blendreg*2+datachunk*2-1),
                               data2.get() + garbage, datachunk-blendreg);

                 for(j=blendreg; j--;)
                 {
                     data1[j+datachunk-blendreg] = 0;
                     data2[j+datachunk-blendreg] = 0;
                 }
             }

             if(index>0)
             {
                for(j=blendreg; j--;)
                    data1[j+garbage] = readdata[j];
             }
         
             for(j=blendreg; j--;)
                readdata[j] = data2[garbage+datachunk-blendreg + j];

             for(j=garbage; j--;)
             {
                data1[j] =0;
                data2[j] =0;
                data1[myWindowSize-garbage+j] =0;
                data2[myWindowSize-garbage+j] =0;
             }

             window(data1, data2, myWindow);

             // do freq transform; multiply; do time transform
             myTransform.toFrequency(
                     data1.get(), data2.get(), myReal.data(), myImag.data(),
                     myWindowSize);

             // Update the filter depending on the current index & chunk
             design_filter(freqstep, index, chunk, myFilter, myPhaseFilter);

             if(!myFilterPhase)
             {
                // Magnitude scale only
                for(j=0; j<myWindowSize; j++)
                {
                    myReal[j] *= myFilter[j];
                    myImag[j] *= myFilter[j];
                }
             }
             else
             {
                // Mag & phase filter
                for(j=myWindowSize; j--;)
                {
                    myReal[j] *= myFilter[j];
                    myImag[j] *= myPhaseFilter[j];
                }
             }

             myTransform.toTime(
                     myReal.data(), myImag.data(), data1.get(), data2.get(),
                     myWindowSize);

             window(data1, data2, myUnWindow);
             
             if(index>0)
             {
                // blend data part 1
                n = blendreg;
                if(index+blendreg >= dest_capacity)
                {
                    if(index >= dest_capacity)
                        n = 0;
                    else
                        n = dest_capacity - index;
                }
                
                if(myBlendChunks)
                {
                    for(j=n; j--;)
                    {
                        v1 = dest[index+j]*(1-blend[j]);
                        v2 = data1[j+garbage] * blend[j];

                        dest[index+j] = v2 + v1;
                    }
                }
                else
                {
                    for(j=n; j--;)
                        dest[index+j] += data1[j+garbage];
                }
                
                // store data part1 & 2
                for(j=blendreg; j<datachunk && index+j<dest_capacity; j++)
                    dest[index+j] =  data1[j+garbage];
                
                for(j=blendreg; j<datachunk && index+datachunk-blendreg+j<dest_capacity; j++)
                    dest[index+datachunk-blendreg+j] = data2[j+garbage];

                // blend data part 2
                if(index+datachunk-blendreg < dest_capacity)
                {
                    n = blendreg;
                    if(index+datachunk >= dest_capacity)
                        n = dest_capacity - (index+datachunk-blendreg);

                    if(myBlendChunks)
                    {
                        for(j=n; j--;)
                        {
                            v1 = dest[index+datachunk-blendreg+j]*(1-blend[j]);
                            v2 = data2[j+garbage] * blend[j];
                            
                            dest[index+datachunk-blendreg+j] = v2 + v1;
                        }
                    }
                    else
                    {
                        for(j=n; j--;)
                            dest[index+datachunk-blendreg+j]+=data2[j+garbage];
                    }
                }
             }
             else
             {
                // store data parts 1
                n = datachunk;
                if(index+datachunk >= dest_capacity)
                {
                    if(index >= dest_capacity)
                        n = 0;
                    else
                        n = dest_capacity - index;
                }
                
                for(j=n; j--;)
                    dest[index+j] = data1[j+garbage];
                
                // blend data part 2
                if(index+datachunk-blendreg < dest_capacity)
                {
                    n = blendreg;
                    if(index+datachunk >= dest_capacity)
                        n = dest_capacity - (index+datachunk-blendreg);

                    if(myBlendChunks)
                    {
                        for(j=n; j--;)
                        {
                            v1 = dest[index+datachunk-blendreg+j]*(1-blend[j]);
                            v2 = data2[j+garbage] * blend[j];

                            dest[index-blendreg+datachunk+j] = v2 + v1;
                        }
                    }
                    else
                    {
                        for(j=n; j--;)
                            dest[index-blendreg+datachunk+j] += data2[j+garbage];
                    }
                }

                // store data part 2
                for(j=blendreg; j<datachunk && index+datachunk-blendreg+j < dest_capacity; j++)
                    dest[index+datachunk-blendreg+j] = data2[j+garbage];
                
             }
        }
    }
    boss->opEnd();
}

template <class F>
void
CL_FreqFilter::constrainedFilter(
        fpreal *dest,
        int dest_capacity,
        fpreal freqstep,
        const F &evaluate_data,
        int constrain_region,
        bool shift)
{
    auto shifted = UTmakeUnique<fpreal[]>(dest_capacity);
    evaluate_data(0, dest_capacity - 1, shifted.get(), dest_capacity);

    fpreal original_start = shifted[0];

    // shift so that start = 0
    // this avoids the initial step response, especially when applied to constant functions 
    if (shift)
    {
        for (exint i = 0; i < dest_capacity; i++)
            shifted[i] -= original_start;
    }

    fpreal start = shifted[0];
    fpreal end = shifted[dest_capacity - 1];

    auto evaluate_shifted = [&](int start, int, fpreal *data, int size)
    {
        for (exint i = 0; i < size; i++)
            data[i] = (start + i) < dest_capacity ? shifted[start + i] : shifted[dest_capacity - 1];
    };

    filter(dest, dest_capacity, freqstep, evaluate_shifted);

    if (constrain_region > dest_capacity / 2)
        constrain_region = dest_capacity / 2;

    fpreal delta_start = start - dest[0];

    fpreal delta_end = end - dest[dest_capacity - 1];

    for (exint i = 0; i < dest_capacity / 2; i++)
    {
        fpreal t = static_cast<fpreal>(i) / constrain_region;
        fpreal a = SYSpow(2.0, -10.0 * t);

        dest[i] += delta_start * a;

        dest[dest_capacity - 1 - i] += delta_end * a;
    }

    // shift back to original range
    if (shift)
    {
        for (exint i = 0; i < dest_capacity; i++)
            dest[i] += original_start;
    }
}

// Designs a filter, populating the filter and phase_filter arrays
template <class F>
void 
CL_FreqFilter::designPassFilter(Type type, Design design, 
                                fpreal freq1, fpreal freq2,
                                fpreal attenuation, fpreal pass_gain, int order,
                                int chunk_size, fpreal freqstep, int index, int num,
                                const F& override_parameters, 
                                fpreal *filter, fpreal *phase_filter)
{
    fpreal      f,k;
    UT_Complex  s,w1,w2,w3, w1_inv, response;
    fpreal      dc_gain,gain, exp1, filter_order;
    int         i;
    int         half = chunk_size/2;
    fpreal      max;

    exp1 = 2.0;
    filter_order = static_cast<fpreal>(order); // Order of the filter. Used for Butterworth filters

    // Potentially overrides the frequency, gain, and myAttenuationuation parameters
    // depending on the index and num
    override_parameters(freq1, freq2, gain, attenuation, index, num);

    if(freq1<=SYS_FTOLERANCE)
        freq1 = SYS_FTOLERANCE;
    if(freq2<=SYS_FTOLERANCE)
        freq2 = SYS_FTOLERANCE;

    switch(type)
    {
    case Type::LOW_PASS:
        w2.set(1,0);
        w1.set(freq2,0); // Cutoff frequency (high)
        w1_inv = w2 / w1;

        s.set(0,freq2/100);

        k = SYSpow(freq2,exp1);

        if (design == Design::BUTTERWORTH)
            // Frequency response of the Butterworth Filter
            // See: https://en.wikipedia.org/wiki/Butterworth_filter#Transfer_function
            response = (w2) / (w2 + (-(s * w1_inv).pow(exp1)).pow(filter_order));
        else
            response = (w2*k)/((s+w1).pow(exp1)); // default type

        dc_gain = response.magnitude();
        if(SYSequalZero(dc_gain))
            dc_gain=1.0;
        gain = 1.0/dc_gain;

        for(i=0,f=0; i<=half; i++,f+=freqstep)
        {
            s.set(0,f);

            if (design == Design::BUTTERWORTH)
            {
                if (f < SYS_FTOLERANCE)
                    s.set(0, SYS_FTOLERANCE);
                response = (w2) / (w2 + (-(s * w1_inv).pow(exp1)).pow(filter_order));
            }
            else
                response = (w2*k)/((s+w1).pow(exp1));

            if(phase_filter)
            {
                filter[i] = response.real();
                phase_filter[i] = response.imaginary();
            }
            else
                filter[i] = response.magnitude()*gain;
        }
        break;

    case Type::HIGH_PASS:
        w1.set(freq1,0);
        w2.set(1, 0);
        s.set(0,freq1*100);

        k =  SYSpow(freq1,exp1);
        if (design == Design::BUTTERWORTH)
            response = (w2) / (w2 + (-(w1 / s).pow(exp1)).pow(filter_order));
        else
            response = (s.pow(exp1)*k) / (s+w1).pow(exp1);

        dc_gain = response.magnitude();
        if(SYSequalZero(dc_gain))
            dc_gain=1.0;

        gain = 1.0/dc_gain;

        for(i=0,f=0; i<=half; i++,f+=freqstep)
        {
            s.set(0,f);

            if (design == Design::BUTTERWORTH)
            {
                if (f < SYS_FTOLERANCE)
                    s.set(0, SYS_FTOLERANCE);

                response = (w2) / (w2 + (-(w1 / s).pow(exp1)).pow(filter_order));
            }
            else
                response = (s.pow(exp1)*k)/(s+w1).pow(exp1);

            if(phase_filter)
            {
                filter[i] = response.real();
                phase_filter[i] = response.imaginary();
            }
            else
                filter[i] = response.magnitude()*gain;

        }
        break;

    case Type::BAND_STOP:
    case Type::BAND_PASS:

        if(freq1>freq2)
        {
            f=freq1;
            freq1=freq2;
            freq2=f;
        }
    
        w1.set(freq1,0);
        w2.set(freq2,0);
        s.set(0,freq1);

        k = freq1 * freq2;
        response =  (s*k) /((s+w1) * (s+w2));

        dc_gain = response.magnitude();
        if(SYSequalZero(dc_gain))
            dc_gain=1.0;
            
        gain = 1.0/dc_gain;

        for(i=0,f=0; i<=half; i++,f+=freqstep)
        {
            s.set(0,f);

            response = (s*k)/((s+w1)*(s+w2));

            if(phase_filter)
            {
                filter[i] = response.real();
                phase_filter[i] = response.imaginary();
            }
            else
                filter[i] = response.magnitude()*gain;
        }
        break;
    }

    if(!SYSisEqual(attenuation,1.0))
    {
        for(i=0; i<=half; i++)
            filter[i] = SYSpow(filter[i],attenuation);
    }

    gain = SYSpow(CONST_FPREAL(10), pass_gain/20);
    max = 0;
    for(i=0; i<=half; i++)
        if(filter[i] > max)
            max = filter[i];

    max = (max!=0)?(gain/max):gain;
    
    for(i=0; i<=half; i++)
        filter[i] *=max;

    // band stop is just 1-bandpass
    if(type == Type::BAND_STOP)
    {
        if(phase_filter)
        {
            for(i=0; i<=half; i++)
            {
                max = SYSsqrt(filter[i]*filter[i] + phase_filter[i]*phase_filter[i]);
                max = (gain-max)/max;
                filter[i]*= max;
                phase_filter[i]*=max;
            }
        }
        else
        {
            for(i=0; i<=half; i++)
                filter[i] = gain - filter[i];
        }
    }
           
    
    for(i=1; i<half; i++)
    {
        filter[chunk_size-i] = filter[i];
        if(phase_filter)
            phase_filter[chunk_size-i] = -phase_filter[i];
    }
}

template <class F>
void 
CL_FreqFilter::designBandEQFilter(Shape shape, fpreal basefreq, fpreal *bands, int num_bands,
                              int chunk_size, fpreal freqstep, int index, int num,
                              const F& override_parameters, 
                              fpreal *filter)
{
    int          i, j;
    fpreal       band;
    int          starti,endi;
    fpreal       freq,amp;
    fpreal       lstart, lend, lfreq, lwidth, lcenter;

    override_parameters(basefreq, bands, index, num);

    band = CONST_FPREAL(10) * SYSpow(CONST_FPREAL(2), basefreq);

    fpreal width = (shape == CL_FreqFilter::Shape::GAUSSIAN) ? 2.0 : SQRT_TWO;

    for(i=0; i< chunk_size; i++)
        filter[i] = 1.0;

    endi = 0;
    for(i=0; i<num_bands; i++)
    {
        if(endi>=chunk_size/2)
            continue;
        
        if(shape == Shape::BOX)
            starti = endi;
        else
            starti = int((band/width)/freqstep + 0.5);
        endi = int((band*width)/freqstep + 1.5);
        
        if(bands[i] != 1.0)
        {
            lstart = SYSlog(band/width);
            lend = SYSlog(band*width);
            lwidth = (lend-lstart)/2.0;
            lcenter = SYSlog(band);

            // filter contraints.
            if(starti>chunk_size/2)
                continue;
            if(endi>chunk_size/2)
                endi = chunk_size/2;
            if(starti<0)
                starti=0;


            for(j=starti; j<=endi; j++)
            {
                lfreq = SYSlog(j*freqstep);
                if(shape == Shape::GAUSSIAN)
                {
                    freq = 3.0*(lfreq-lcenter)/lwidth;
                
                    amp = SYSexp( -0.5 *freq*freq);
                    amp = (bands[i] - 1.0) * amp;

                    filter[j] +=amp;
                }
                else
                {
                    filter[j] = bands[i];
                }
            }
        }
            
        band *=2;
    }

    for(i=1; i<chunk_size/2; i++)
        filter[chunk_size-i] = filter[i];

}

template <class F>
void 
CL_FreqFilter::designParaEQFilter(fpreal shape, fpreal freq, fpreal fgain, fpreal pgain,
                                  fpreal width, fpreal atten,
                                  int chunk_size, fpreal freqstep, int index, int num,
                                  const F& override_parameters, 
                                  fpreal *filter)
{
    int          i, j;
    int          starti,endi;
    fpreal       amp,ind;
    fpreal       lstart, lend, lfreq, lwidth, lcenter;

    override_parameters(freq, width, fgain, pgain, shape, atten, index, num);

    for(i=0; i<chunk_size; i++)
        filter[i] = pgain;

    width = SYSpow(2.0,width/2.0);

    starti      = int((freq/width)/freqstep + 0.5);
    endi        = int((freq*width)/freqstep + 1.5);
    
    // filter contraints.
    if(starti>chunk_size/2)
        return;
    if(endi>chunk_size/2) 
        endi = chunk_size/2;
    if(starti<0)
        starti=0;

    lstart = SYSlog(freq/width);
    lend = SYSlog(freq*width);
    lcenter = SYSlog(freq);
    lwidth = (lend - lstart)/2.0;

    ind = freq/width;
    for(j=starti; j<=endi; j++, ind+=freqstep)
    {
        lfreq = SYSlog(ind);

        if(shape <= 0.5)
        {
            // blend between box (0) and triangle (0.5)
            amp = (1.0-SYSabs(lfreq-lcenter)/lwidth);
            
            amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
            amp = fgain*(0.5-shape)*2 + amp*(shape*2);
            
            filter[j] =amp;
        }
        else
        {
            // blend between triangle (0.5) and gaussian (1.0)
            freq = 3.0*(lfreq-lcenter)/lwidth;
                    
            amp = SYSexp( -0.5 *freq*freq);
            amp = (fgain - pgain)*attenuate(amp,atten) + pgain;
            filter[j] = amp*(shape-0.5)*2;

            amp = 1.0-SYSabs(lfreq-lcenter)/lwidth;
            amp = (fgain-pgain)*attenuate(amp,atten) +pgain;
            filter[j] += amp*(1.0-shape)*2;
        }
    }

    for(i=1; i<chunk_size/2; i++)
        filter[chunk_size-i] = filter[i];

}

template <class F>
void 
CL_FreqFilter::designAverageSideband(fpreal fgain, fpreal gain, fpreal effect, int length, int passfreq,
                                     int chunk_size, const F& evaluate_data, 
                                     fpreal *filter)
{
    UT_FFT<fpreal> freqtransform(chunk_size);

    int          i,j,count;
    fpreal       factor;
    int          start,end;
    double       max;

    auto temp = UTmakeUnique<fpreal[]>(chunk_size);
    auto real = UTmakeUnique<fpreal[]>(chunk_size);
    auto imag = UTmakeUnique<fpreal[]>(chunk_size);

    for(j=1; j<chunk_size; j++)
        filter[j]=0;

    
    for(i=0,count=0; i<length; i+=chunk_size,count++)
    {
        evaluate_data(
                (fpreal)i, (fpreal)i + chunk_size - 1, temp.get(), chunk_size);

        // do freq transform; multiply; do time transform
        freqtransform.toFrequency(temp.get(),real.get(),imag.get(),chunk_size);

        // calculate power spectrum
        for(j=1; j<chunk_size; j++)
            temp[j] = (real[j]*real[j] + imag[j]*imag[j]);

        for(j=1; j<chunk_size; j++)
            filter[j] += temp[j];
    }

    if(count>1)
        for(j=1; j<chunk_size; j++)
            filter[j] /= count;

    // 'effect' widens the sideband filter on a sample basis
    if(!SYSisEqual(effect,1.0))
    {
        for(j=(chunk_size>>1); j--;)
        {
            temp[j] = filter[j];
            filter[j] = 0;
        }

        for(j=(chunk_size>>1); j--;)
        {
            start = int(j/effect);
            end = int(j*effect);

            if(start < 0)
                start=0;
            if(end > chunk_size/2)
                end = chunk_size/2;

            for(i=start; i<=end; i++)
            {
                // geometric interp
                factor = SYSlog(fpreal(i)/fpreal(j))/SYSlog(effect);

                // gaussian distribution
                factor = SYSexp(-1.5*factor*factor);
                filter[i] += temp[j] * factor;
            }
        }

        for(j=1; j<(chunk_size>>1); j++)
            filter[chunk_size-j] = filter[j];
    }
    // don't filter DC
//    filter[0] = 1.0;
    
    for(j=0; j<chunk_size; j++)
        if(j==0 || filter[j]>max)
            max = filter[j];

    if(!SYSequalZero(max))
        for(j=0; j<chunk_size; j++)
            filter[j] /=max;

    if(!SYSisEqual(fgain,1.0))
        for(j=0; j<chunk_size; j++)
            filter[j] *= fgain; 

    if(passfreq)
    {
        for(j=0; j<chunk_size; j++)
            filter[j] = (1.0+gain) - filter[j];
    }
    else
    {
        for(j=0; j<chunk_size; j++)
            filter[j] += gain;
    }
}

template <class F>
void 
CL_FreqFilter::designContinuousSideband(fpreal fgain, fpreal gain, fpreal effect, int passfreq,
                                        int chunk_size, int index,
                                        const F& evaluate_data, 
                                        fpreal *filter)
{
    UT_FFT<fpreal> freqtransform(chunk_size);

    int    i,j;
    int    start,end;
    fpreal factor;
    double max;

    auto real = UTmakeUnique<fpreal[]>(chunk_size);
    auto imag = UTmakeUnique<fpreal[]>(chunk_size);

    evaluate_data((fpreal)index, (fpreal)index+chunk_size-1, filter, chunk_size);

    // do freq transform; multiply; do time transform
    freqtransform.toFrequency(filter,real.get(),imag.get(),chunk_size);

    // calculate power spectrum
    for(j=1; j<chunk_size; j++)
        filter[j] = (real[j]*real[j] + imag[j]*imag[j]);

    // 'effect' widens the sideband filter on a sample basis
    if(!SYSisEqual(effect,1.0))
    {
        for(j=0; j<=chunk_size/2; j++)
        {
            real[j] = filter[j];
            filter[j] = 0;
        }

        for(j=0; j<chunk_size/2; j++)
        {
            start = int(j/effect);
            end = int(j*effect);

            if(start < 0)
                start=0;
            if(end > chunk_size/2)
                end = chunk_size/2;

            for(i=start; i<=end; i++)
            {
                // geometric interp
                factor = SYSlog(fpreal(i)/fpreal(j))/SYSlog(effect);

                // gaussian distribution
                factor = SYSexp(-1.5*factor*factor);
                filter[i] += real[j] * factor;
            }
        }

        for(j=1; j<chunk_size/2; j++)
            filter[chunk_size-j] = filter[j];
    }
    
    // don't filter DC
//    filter[0]=1;
    
    for(j=0; j<chunk_size; j++)
        if(j==0 || filter[j]>max)
            max = filter[j];

    if(!SYSequalZero(max))
        for(j=0; j<chunk_size; j++)
            filter[j] /=max;
    
    if(!SYSisEqual(fgain,1.0))
        for(j=0; j<chunk_size; j++)
            filter[j] *= fgain; 
    
    if(passfreq)
    {
        for(j=0; j<chunk_size; j++)
            filter[j] = (gain+1.0) - filter[j];
    }
    else
    {
        for(j=0; j<chunk_size; j++)
            filter[j] += gain;
    }
}

#endif